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