bam文件softclip , hardclip ,markduplicate的探究
来源:互联网 发布:贵金属行情分析软件 编辑:程序博客网 时间:2024/06/06 13:16
测序产生的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
阅读全文
0 0
- bam文件softclip , hardclip ,markduplicate的探究
- bam文件softclip , hardclip ,markduplicate的探究
- 多线程拆分bam文件
- WSO2-BAM的FAQ
- sam to bam文件出错
- Linux文件锁的探究
- NGUI UIPanel在softClip下不起作用的问题
- Oracle BAM Sensor 的enable batching属性
- Oracle BAM的Drill Down实现
- Oracle BAM Export Rule的方式
- bam/sam 数据格式的介绍 (一)
- bam/sam 数据格式的介绍 (二)
- BAM部署失败 - 未能加载”AdomdClient”或它的某一个依赖项。系统找不到指定的文件
- 对python文件方法open的探究
- 对python文件方法open的探究
- ODI 11g和Oracle BAM的集成
- Oracle BAM发送报警信息超长的解决方法
- NGUI裁剪粒子的实现,完全仿照原有的NGUI SoftClip逻辑思路,可多层裁剪
- Confluence集成实践:3:RestAPI集成方式
- Java常用工具类封装——String操作工具类
- 总结Java中的next()和nextLine()方法的区别
- UNIX环境编程学习笔记(13)——文件I/O之标准I/O流
- RobotFramwork(一)安装篇
- bam文件softclip , hardclip ,markduplicate的探究
- fabirc1.0商业正式版本源码解析11---peer的Admin和Endorser服务
- 常用SQL查询语句
- AtCoder Beginner Contest 072 ABCD C++&&Python3
- LinkedHashMap分析
- Dubbo服务和EDAS服务互调
- jstl时间转换
- Mac系统常用命令
- JAVA 并发编程随笔【五】Thread线程创建及运行线程任务