ARTで250bpでカバレッジ100のシングルfastqを生成。マッピングソフトによるカバレッジの差を調べる。
まずはfastqのジェネレート。
art_illumina -ss MSv3 -sam -i input.fasta -p -l 250 -f 100 -s 10 -o single-read
bwa alnは感度の面から使えないように見える。ただし、リアルデータだと少し話が変わってくることに気づいた。
次はARTのリファレンスにしたゲノムを実際にMiseqで読んだ時のデータ。表にはゲノム全体の平均カバレッジを記載した。
(カバレッジ = アップされたリード数 x 平均リード長 / ゲノムサイズ)
以前読んだシーケンスデータも3つ調べてみた。
bwa alnのマッピング率がARTのデータで極端に落ち込んだ原因は、ARTの fastqが5'側にlow qualityの部位をつくりやすいため。これは MiSeqのデータで5'側のクオリティが落ち込むことを反映しているためらしい。
ARTの人口データでも、bra alnのシード領域のミスマッチを2まで許容させるとマッピング率は7.1から63.5まで増加した。
変異部位
一番落ち込んでいるのはbwa alnの標準マッピング条件(中央が変異)。エラー許容条件ではindelでもカバレッジが落ち込まない。
最後にmismatch permitting valuの大きさの分析、数値が大きければミスマッチがより大量に蓄積していることを意味する
打ったコマンドは以下の通り。
---------------------------------------------------------------------------------------------------------------
デフォルト条件
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
--------------------------------------------------------------------------------