【直播】我的基因组23:对比对结果文件进行过滤
在此之前,我们分别在十七,十八讲中讲到了multiple mapping,PCR duplication的情况,在分析我的基因组时呢,肯定要去除掉multiple mapping、PCR duplication及低质量比对结果。那么如何简单的去掉它们呢?
用管道,直接一个命令就搞定了,是不是很简单呢!
samtools view -h -F4 -q 5 P_jmzeng.final.bam |samtools view -bS |samtools rmdup -o P_jmzeng.filter.rmdup.bam
samtools index P_jmzeng.filter.rmdup.bam
过滤和过滤后的文件大小相差了12个G。
肯定会有人问什么是管道呢?管道命令就是利用Linux所提供的管道符“|”将两个命令隔开,管道符左边命令的输出就会作为管道符右边命令的输入。这样大家是不是好理解一些了?
bam文件是二进制文件,我们需要samtools view的命令进行格式转换。这个管道的意思分开来说就是运行第一步时过滤的时候已将bam文件转成我们能看的sam格式。其中h是在输出的结果中包含头header,F和q是过滤掉没有mapped上的reads(也就是multiple mapping的情况)和低质量的reads。第二步是将上一步得出的sam文件再转成bam,第三步就是用samtools rmdup过滤掉PCR duplication的情况了,最后得到了过滤了multiple mapping、PCR duplication及低质量比对的bam文件。
最后利用samtools index对过滤后的bam文件建立索引。
当然,其实这个过滤不一定那么重要,因为很多需要bam文件作为输入文件的软件有设置参数问你是否需要不考虑那些低质量的,PCR duplication的,或者multiple mapping的reads。
但是那些软件的这些参数默认是不开启的,而大多数人跑程序又喜欢用默认参数,所以这个步骤对他来说,就异常重要了。
文:Jimmy、阿尔的太阳
图文编辑:吃瓜群众