macでインフォマティクス

macでインフォマティクス

HTS (NGS) 関連のインフォマティクス情報についてまとめています。

InDel_Hunterのマッピングソフト検討

ARTで250bpでカバレッジ100のシングルfastqを生成。マッピングソフトによるカバレッジの差を調べる。

 

まずはfastqのジェネレート。

art_illumina -ss MSv3 -sam -i input.fasta -p -l 250 -f 100 -s 10 -o single-read

 

マッピングソフトデフォルト条件での平均カバレッジ

f:id:kazumaxneo:20170324172706j:plain

 

エラー許容マッピング条件での平均カバレッジ

f:id:kazumaxneo:20170324172736j:plain

 

bwa alnは感度の面から使えないように見える。ただし、リアルデータだと少し話が変わってくることに気づいた。

 

次はARTのリファレンスにしたゲノムを実際にMiseqで読んだ時のデータ。表にはゲノム全体の平均カバレッジを記載した。

(カバレッジ = アップされたリード数 x 平均リード長 / ゲノムサイズ)

 

 

以前読んだシーケンスデータも3つ調べてみた。

f:id:kazumaxneo:20170324172844j:plain

 f:id:kazumaxneo:20170324172850j:plain

f:id:kazumaxneo:20170324172857j:plain

 bwa alnのマッピング率がARTのデータで極端に落ち込んだ原因は、ARTの fastqが5'側にlow qualityの部位をつくりやすいため。これは MiSeqのデータで5'側のクオリティが落ち込むことを反映しているためらしい。

 

 

 ARTの人口データでも、bra alnのシード領域のミスマッチを2まで許容させるとマッピング率は7.1から63.5まで増加した。

変異部位

f:id:kazumaxneo:20170324173027j:plain

一番落ち込んでいるのはbwa alnの標準マッピング条件(中央が変異)。エラー許容条件ではindelでもカバレッジが落ち込まない。

 

 

最後にmismatch permitting valuの大きさの分析、数値が大きければミスマッチがより大量に蓄積していることを意味する

f:id:kazumaxneo:20170324173017j:plain

 

 

 

f:id:kazumaxneo:20170324173122j:plain

 

 

 

打ったコマンドは以下の通り。

---------------------------------------------------------------------------------------------------------------

デフォルト条件

 

1、bwa aln

bwa index -p $f -a is $f

bwa aln -t 20 $f $p > R1.sai

bwa samse -r "@RG\tID:HSS\tSM:WT\tPL:illumina" $f R1.sai $p > R1.sam

 

2、bwa mem

bwa index -p $f -a is $f

bwa mem -R "@RG\tID:X\tLB:Y\tSM:Z\tPL:ILLUMINA" -t 22 $f $p $q > bwa_mem_R1.sam

 

3、smalt

smalt index -k 14 -s 8 k14s8 $f

smalt map -n 22 -o R1.sam k14s8 $p

 

4、bowtie2

bowtie2-build $f R1_index

bowtie2 -x R1_index -U $p -S R1.sam -k 100 -p 20

 

 

---------------------------------------------------------------------------------------------------------------

エラー許容条件

 

1、bwa aln

bwa index -p $f -a is $f

bwa aln -l 25 -k 0 -n 225 -i 225 -t 20 $f $p > R1.sai

bwa samse -r "@RG\tID:HSS\tSM:WT\tPL:illumina" $f R1.sai $p > R1.sam

 

2、bowtie2

bowtie2-build $f R1_index

bowtie2 -x R1_index -U $p -L 25 --gbar 225 --mp 2 --end-to-en --np 0 -N 0 -S R1.sam -k 100 -p 24

 

 

--------------------------------------------------------------------------------