mappingを行ったbamファイルから更に絞り込みをかけてやる。
具体的にはSNV、INDELを持つreadを取り除いて100%matchしたreadだけにしたい。
基本はsamtools viewでbamを開いて、grepやawkで整形をすることになるだろう。
bamはバイナリなので中身が見れないが、samtools view -h sample.bam > sample.samとsamファイルに戻してやって、ファイルをテキストエディタで開いてやると・・・
@HD VN:1.3 SO:coordinate @SQ SN:AT1G51370.2 LN:1118 ・・・・・・ @SQ SN:ATMG00930.1 LN:74 @PG ID:bwa PN:bwa VN:0.7.15-r1140 CL:bwa mem -t8 /home/kuro/rnaseq/ref/arabidopsis_thaliana/bwa/TAIR10_cdna_20101214_updated ./trimmed/sample_L001_R1_001_trim_np.fastq.gz @PG ID:bwa-4217DAA5 PN:bwa VN:0.7.15-r1140 CL:bwa mem -t8 /home/kuro/rnaseq/ref/arabidopsis_thaliana/bwa/TAIR10_cdna_20101214_updated ./trimmed/sample_L002_R1_001_trim_np.fastq.gz @PG ID:bwa-15F23362 PN:bwa VN:0.7.15-r1140 CL:bwa mem -t8 /home/kuro/rnaseq/ref/arabidopsis_thaliana/bwa/TAIR10_cdna_20101214_updated ./trimmed/sample_L003_R1_001_trim_np.fastq.gz @PG ID:bwa-709CFACE PN:bwa VN:0.7.15-r1140 CL:bwa mem -t8 /home/kuro/rnaseq/ref/arabidopsis_thaliana/bwa/TAIR10_cdna_20101214_updated ./trimmed/sample_L004_R1_001_trim_np.fastq.gz NB501315:33:HGYFWBGXY:4:13505:8929:6100 16 AT1G50920.1 104 60 85M * 0 0 TATCGATAAGCGAAGATGGTGCAATACAATTTCAAGAAGATCACTGTTGTTCCCAATGGCAAGGAGTTCATTGACATCATCCTCT EEEEEEEEEEEEEEAEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE//EEEEEEEEEEEEEEEAEEEEEEEEEEEEEAAAAA NM:i:8 MD:Z:0C19T5T10G6A14G9G13T1 AS:i:52 XS:i:0 ・・・・・・
こんな感じ。なおこのmappingはリファレンスにアラビのcDNAを使っているので@SQの行が数万行は並んでいる。
@PGがmappingの実施状況を示していて、この場合４つのfastqファイルをTAIR10_cdna_20101214_updatedをリファレンス配列としてマッピングしていることがわかる。
さて、リード情報本体の意味だが
1列目 : NB501315:33:HGYFWBGXY:4:13505:8929:6100 リードの名前
2列目 : 16 flag:リードの状況
3列目：AT1G50920.1 張り付いた染色体，コンティグの名前
4列目：104 張り付いた場所
5列目：60 マッピングスコア
6列目：85M マッピング状況 この場合85ベースがマッチ(M)
7列目：* paired-endの時の相方の名前（single-endの場合*）
8列目：0 paired-endの時の相方の場所
9列目：0 paired-endの時のインサートの長さ
10列目：TATCGATAAGCGAAGATGGTGCAATACAATTTCAAGAAGATCACTGTTGTTCCCAATGGCAAGGAGTTCATTGACATCATCCTCT リードのシークエンス配列
11列目：EEEEEEEEEEEEEEAEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE//EEEEEEEEEEEEEEEAEEEEEEEEEEEEEAAAAA リードのクオリティ
12列目 : NM:i:8 MD:Z:0C19T5T10G6A14G9G13T1 AS:i:52 XS:i:0 オプションフィールド（TAG:VTYPE:VALUE）
2列目のリード状況はいろいろあるようだが、0:マップされている 16:逆相補鎖にマップされている 4:きちんとマップされていない 20:逆相補鎖にきちんとマップされていないの4つが主だったものとなるようだ。
きちんとマップされたリードだけをまずは抽出するのが1つ目。
次に12列目の部分がBWAによって付与されたタグらしく、これを次の抽出ターゲットとする。
BWA generates the following optional fields. Tags starting with ‘X’ are specific to BWA.
Tag Meaning
NM Edit distance
MD Mismatching positions/bases
AS Alignment score
BC Barcode sequence
X0 Number of best hits
X1 Number of suboptimal hits found by BWA
XN Number of ambiguous bases in the referenece
XM Number of mismatches in the alignment
XO Number of gap opens
XG Number of gap extentions
XT Type: Unique/Repeat/N/Mate-sw
XA Alternative hits; format: (chr,pos,CIGAR,NM;)*
XS Suboptimal alignment score
XF Support from forward/reverse alignment
XE Number of supporting seeds
NB501315:33:HGYFWBGXY:1:11308:11940:17581 0 AT1G71100.1 166 52 86M * 0 0 CTCTAACCTCACCCAAGAAGAACTCAAGAGAATCGCCGCCTATAAGGCTGTCGAATTCGTCGAATCTGGAATGGTCATTGGTCTCG AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE</EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEE NM:i:8 MD:Z:3A8G16A9T5A2C26T2C7 AS:i:47 XS:i:0
NB501315:33:HGYFWBGXY:1:23308:20616:13305 0 AT1G71100.1 166 60 86M * 0 0 CTCAAACCTCACGCAAGAAGAACTCAAGAAAATCGCCGCTTATAAAGCCGTCGAATTCGTCGAATCTGGAATGGTTATCGGTCTCG AAAAAEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEEEEEEEEEEEEEEEEE NM:i:0 MD:Z:86 AS:i:86 XS:i:0
下がmutationなどなにもないperfect matchで上がmismatchたくさんのread
MD:Z:3A8G16A9T5A2C26T2C7
の意味するところは数字のところは連続でそれだけmatchでATGCのところはその塩基に置換しているという感じのようだ。
NB501315:33:HGYFWBGXY:4:11608:13474:9299 16 AT1G48450.1 196 0 54M32S * 0 0 GGCTTCATGCGATTCGGCTCCGGACCTCAATTTCGTCTCAGACAAAAGTCTTCTAGACCATTGCAAAGCCGAACTAGTAATGTGAG EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEAAAAA NM:i:4 MD:Z:12A1T0A31C6 AS:i:34 XS:i:34 XA:Z:AT1G48450.2,-191,54M32S,4;AT1G48450.3,-196,54M32S,4;
例えばこれだとXAでaltanative hitsが複数ある他、86塩基中54塩基でアライメントするが32塩基は無視されている。
さて、どういう書式で処理すればミスマッチを含むリードを除くことができるだろうか？
最も簡単かつ、厳しい抽出条件としては
samtools view sample_1_sort.bam | grep "MD:Z:86" > sample_1_sort_ext.bam
という感じだろうか。これだとヘッダ情報が失われてまずいかもしれないので、一旦samに書き出して、text整形してからbamに戻してやったほうがいいかも。
ちなみに先程のmismatchてんこ盛りのファイルについて、処理前
$ samtools view sample_1_sort.bam | wc -l 42330669
なのが
$ samtools view sample_1_sort.bam | grep "MD:Z:86" | wc -l 113159
というふうに4200万リード中11万リードが残った。さてどんなものが残ったのかな。
