如何统计BAM文件中的reads数 11

当完成测序的比对工作之后,我们得到了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数。具体请参考:NGS分析入门:操作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,其中即NGS分析入门:操作SAM/BAM文件文中所谓的标记(第2列)。这些标记的具体描述如下:

Flag Chr Description
0x0001 p the read is paired in sequencing
0x0002 P the read is mapped in a proper pair
0x0004 u the query sequence itself is unmapped
0x0008 U the mate is unmapped
0x0010 r strand of the query (1 for reverse)
0x0020 R strand of the mate
0x0040 1 the read is the first read in a pair
0x0080 2 the read is the second read in a pair
0x0100 s the alignment is not primary
0x0200 f the read fails platform/vendor quality checks
0x0400 d the read is either a PCR or an optical duplicate

从该表中,我们可以知道,0x0004标记表明了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 0x0004)以及mate是unmapped(flag 0x0008)。0x0004 + 0x0008 = 12. 小工具提供给大家。关于samtools的更多有趣的事情,可以阅读:A Sourceful of Secrets

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

11 thoughts on “如何统计BAM文件中的reads数

  1. Reply Puriney 12月 30,2013 1:35 下午

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

  2. Reply Dakz 5月 27,2014 12:15 下午

    http://rous.mit.edu/index.php/Manipulate_alignment_files_using_UNIX_commands

    这里有不少基于unix系统自带命令对sam处理的技巧,可以统计比如每个染色体分别比对上了多少reads等。

  3. Reply hailongcorn 7月 22,2014 1:53 上午

    请问博主,我用bedtools multicov高通量计算基因的reads数,怎么进行才能得到的不是aligment而是reads数

  4. Reply hailongcorn 7月 24,2014 10:03 下午

    上面不是说,aligment要比实际read数要多,我想计算从bam里计算某一段intron的read数,用multicov得到的不是实际的reads数吧,我现在在用htseq-count做,还请博主提供点建议

    • Reply admin 7月 25,2014 8:00 上午

      做count时,我一般也是用htseq-count。大多数人也都是推荐这个。而对于multicov,没有仔细研究过。aligments可能要比实际的reads数,如果说一个reads被align到多个位点的话,可能是aligments要多些。但我想这个时候,因为你没有做not uniquely mapped reads过滤的话,按照htseq_count上讲述的,它也不一定会正确处理,注意其中的描述:__alignment_not_unique: reads (or read pairs) with more than one reported alignment. These reads are recognized from the NH optional SAM field tag. (If the aligner does not set this field, multiply aligned reads will be counted multiple times, unless they get filtered out by due to the -a option.)

  5. Reply hailongcorn 7月 25,2014 3:25 上午

    而且我的intron是成批量的,多谢博主

  6. Reply hailongcorn 7月 28,2014 3:14 上午

    多谢博主

  7. Reply newer 3月 7,2015 11:09 下午

    您好,博主,我想请教您一下,我现在拿到了几个个体的bam文件,已经是别人处理过,现在我想知道这几个个体分别的X数,怎么能看到呢

  8. Reply Dulecha 3月 8,2015 9:08 下午

    博主,我需要统计sam/bam文件中的sequence覆盖率,需要得到覆盖的reads数和全部的reads数,怎么得到?覆盖的reads数就是文中的samtools view -c -f 1 -F 12 test.bam吗?

Leave a Reply

  

  

  

%d 博主赞过: