统计BAM文件中的reads数

评论22,277

当完成测序的比对工作之后,我们得到了bam/sam文件。那么,如何得到reads的统计数据呢?
这有很多途径:
1.读取日志文件。对于bowtie的日志,其中会包括如下的描述:

31991083 reads; of these:
31991083 (100.00%) were unpaired; of these:
6844445 (21.39%) aligned 0 times
18391269 (57.49%) aligned exactly 1 time
6755369 (21.12%) aligned >1 times
78.61% overall alignment rate

从中,我们就可以了解到测序文件的总reads数为31991083,其中100%是unpaired。其中,21.39%没有比对成功,有57.49%确定为一个比对位点,有21.12%有多余一处位点。

2.使用samtools查看。
这其中有许多办法。如果已经使用了samtools sort以及index命令的话,我们可以使用samtools idxstats来获取reads数。具体请参考:操作SAM/BAM文件处理

第二种办法是通过samtools view工具加上-c 选项来实现。例如:

samtools view -c test.bam
31991083

这里需要强调的是:

“Instead of printing the alignments, only count them and print the total number.” samtools view -c 算的不是reads,而是alignments数。当read有很多位置可以align上同时又都输出了,用samtools view -c 会比实际reads树木要多~~~
by: 严云

如果我们想知道有多少是mapped,有多少是没有mapped,我们可以通过加上-F 4选项或者-f 4选项。这其实就是使用了samtools view的过滤器。在SAM文件中有很多FLAG,其中即操作SAM/BAM文件处理文中所谓的标记(第2列)。这些标记的具体描述如下:

FlagChrDescription
0×0001pthe read is paired in sequencing
0×0002Pthe read is mapped in a proper pair
0×0004uthe query sequence itself is unmapped
0×0008Uthe mate is unmapped
0×0010rstrand of the query (1 for reverse)
0×0020Rstrand of the mate
0×00401the read is the first read in a pair
0×00802the read is the second read in a pair
0×0100sthe alignment is not primary
0×0200fthe read fails platform/vendor quality checks
0×0400dthe read is either a PCR or an optical duplicate

从该表中,我们可以知道,0×0004标记表明了read比对不成功。而samtools中的view工具使用-f INT来依照INT与标记致保留reads,使用-F INT来跳过reads。当然我们还可使用它们来过滤其它的信息。比如我们想知道有多少paired end reads有mate并且都有map时,可以使用-f 1 -F 12来过滤。

samtools view -c -f 1 -F 12 test.bam

其中-f 1指定只包含那些paired end reads,-F 12是不包含那些unmapped(flag 0×0004)以及mate是unmapped(flag 0×0008)。0×0004 + 0×0008 = 12. 小工具提供给大家。关于samtools的更多有趣的事情,可以阅读:A Sourceful of Secrets

有了这个记数,我们就可以计算cuffdiff中的FPKM值了。

原文来自:http://pgfe.umassmed.edu/ou/archives/3419

发表评论

匿名网友