測序產生的bam文件，有一些reads在cigar值里顯示存在softclip，有一些存在hardclip，究竟softclip和hardclip是怎么判斷出來的，還有是怎么標記duplicate的reads的，我懷着這些問題進行了探究。
測試步驟
- 編輯兩個bed文件，分別含有我們需要的read1和read2位置，這里每個文件包含兩條read1或者兩條read2，read1、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