转载自嘻皮悠的博客 http://blog.sina.com.cn/u/2054060511
————————————————————————————————————————————————
reads map到基因组上通常会产生bam/sam文件,但bam/sam文件中通常不直接给出reads mapping的正负链(+/-)信息。那么,如何获取所有mapping到基因组正链或者负链的reads呢?
从sam文件的说明文档中(http://samtools.github.io/hts-specs/SAMv1.pdf),我们得知,FLAG 0x10(十六进制的10,相当于十进制的16)表示reads是否显示为反向互补。因此,只需要对该位进行过滤,即可知道read到底是map到了正链还是负链上。
第一种方法:(借鉴https://www.biostars.org/p/59388/)
(1)获取所有mapping的reads
samtools view -F 4 reads.bam >mapped_reads.sam
(2)正链mapping的reads
gawk '(and(16, 2))' mapped_reads.sam > forward_mapped_reads.sam
第二种方法:(来自https://www.biostars.org/p/14378/)
samtools view -F 20 ... : forward strand
samtools view -f 16 ... : reverse strand
samtools view -f 4 ... : unmapped
FLAG信息解读网站:https://broadinstitute.github.io/picard/explain-flags.html