bam文件softclip ， hardclip ，markduplicate的探究

  測序產生的bam文件，有一些reads在cigar值裏顯示存在softclip，有一些存在hardclip，究竟softclip和hardclip是怎麼判斷出來的，還有是怎麼標記duplicate的reads的，我懷着這些問題進行了探究。

測試步驟

• 編輯兩個bed文件，分別含有咱們須要的read1和read2位置，這裏每一個文件包含兩條read1或者兩條read2，read一、read2一對做爲原始的reads（序列名primer_pri），另外一對做爲截取的材料(這裏取序列名爲other)
• 使用bedtools getfasta，從參考基因組得到reads的序列，將read2反向互補。將原始reads放入兩個文件，一個test_R1.fa，一個test_R2.fa

• 在test_R1.fa中添加其它修改過的原始reads，並在test_R2.fa中也添加相應的read2，不過read2不修改

  read1名稱以下

• primer_pri：原始read
• pimer_duplicate1：primer_pri的重複，如出一轍
• pimer_duplicate2：read1 primer_pri去掉5‘兩個鹼基
• pimer_duplicate3：read1 primer_pri去掉5’兩個鹼基，再去掉3'兩個鹼基
• pimer_changeR2Termi5base：read1修改了read2 5‘端的鹼基
• primer_halfother：read1截掉後面reads，用other的5‘部分reads補全
• pimer_change3Termi5base_change5Termi5base_sametwo：read1和read2同樣，而且5'端和3‘端都改變了5個鹼基

• pimer_change3Termi5base_change5Termi5base：read1 5'端和3‘端都改變了5個鹼基，可是read2保留primer_pri的read2

結果

softclip和hardclip

其中

• primer_halfother read1 82M65S，有SA tag，SA:Z:chr12,5378700,+,79S68M,60,0
• pimer_change3Termi5base_change5Termi5base_sametwo 兩條reads均爲5S137M5S
• pimer_change3Termi5base_change5Termi5base read1 5S137M5S
• primer_halfother read1 79H68M ，有SA tag，SA:Z:chr12,5378502,+,82M65S,60,0

結論（部分分析參考SAMv1.pdf文件）

• 對於map到一個位置的read，兩端map不上的叫作clip ，map到一個位置的狀況下以softclip顯示(好比 pimer_change3Termi5base_change5Termi5base_sametwo和 pimer_change3Termi5base_change5Termi5base read1)
• 對於嵌合比對的read（能夠map到多個區域，而且這比對上的區域很大部分非overlap），好比primer_halfother read1比對上兩個位置，一個比對到chr12 : 5378502，一個比對到chr12：5378700，而且兩次hit的位置的鹼基overlap少，產生的這種狀況是由於read前一部分比對到了前者，然後一部分又能夠比對到了後者，所以不管比對到任何一個位置都這條read都是部分match（這種叫作Chimeric alignment/嵌合比對）。
• 嵌合比對的read，有一條是最優的read，由於咱們map的時候設置了-M參數，所以認爲較短的split的reads判定爲優，這裏是的62 clip 的hit判定爲優。所以65個比對不上的顯示爲softclip，而另一個hit，79 clip顯示爲hard clip，序列中不顯示，而且存入0x800(supplementary alignment flag)
• 爲何82M65S對應的是79H68M呢，理論上應該是82H65M纔對，這是由於這裏兩個比對有三個鹼基的overlap，因此前面有65+3個match，後面有79+3個match（製造reads的時候碰巧截取的primer read 3'端三個鹼基和截取的other read 5‘部分read 三個鹼基同樣）
• 這種嵌合比對的reads含有SA tag

duplicate

其中mark爲duplicate 的reads 對（duplicate 是按fragment算） 有 primer duplicate1，primer_duplicate3,pimer_changeR2Termi5base，primer_halfother（82M65S，144M（未改的read2）），pimer_change3Termi5base_change5Termi5base

不屬於duplicate的有
primer_pri，pimer_duplicate2，primer_halfother（79H68M，一條），pimer_change3Termi5basechange5Termi5base_sametwo

結論

• fragment的start和end同樣(read1和read2（由於read2是測對鏈的，reads的5‘端都是fragment的末端）的5’的位置都相同)判斷爲duplicate，只取最優的不標記爲duplicate
• primer_pri的duplicate是 primer duplicate1, primer_halfother
• pimer_duplicate2的duplicate是primer_duplicate3,pimer_changeR2Termi5base，pimer_change3Termi5base_change5Termi5base
• 沒有duplicate的是primer_halfother（79H68M，一條），pimer_change3Termi5basechange5Termi5base_sametwo
• pimer_change3Termi5basechange5Termi5base_sametwo 5'有5 softclip，map的位置從M的鹼基算起(見圖），因此它沒有duplicate