NGS分析入门:操作SAM/BAM文件 40

当测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B取自binary)。

那么SAM文件的格式是什么样子的呢?如果你想真实地了解SAM文件,可以查看它的说明文档。SAM由头文件和map结果组成。头文件由一行行以@起始的注释构成。而map结果是类似下面的东西:

HWI-ST1001:137:C12FPACXX:7:1115:14131:66670     0       chr1    12805   1       42M4I5M *       0       0       TTGGATGCCCCTCCACACCCTCTTGATCTTCCCTGTGATGTCACCAATATG     CCCFFFFFHHGHHJJJJJHJJJJJJJJJJJJJJJJIJJJJJJJJJJJJIJJ     AS:i:-28        XN:i:0  XM:i:2  XO:i:1XG:i:4   NM:i:6  MD:Z:2C41C2     YT:Z:UU NH:i:3  CC:Z:chr15      CP:i:102518319  XS:A:+  HI:i:0
 
HWI-ST1001:137:C12FPACXX:7:2313:17391:30032     272     chr1    13494   1       51M     *       0       0       ACTGCCTGGCGCTGTGCCCTTCCTTTGCTCTGCCCGCTGGAGACAGTGTTT     CFFFFHHJJJJIJJJJIJJJJJJJJJJJJJJJJJJJJJHHHHFA+FFFC@B     AS:i:-3 XN:i:0  XM:i:1  XO:i:0  XG:i:0NM:i:1   MD:Z:44G6       YT:Z:UU XS:A:+  NH:i:3  CC:Z:chr15      CP:i:102517626  HI:i:0
 
HWI-ST1001:137:C12FPACXX:7:1109:17518:53305     16      chr1    13528   1       51M     *       0       0       CGCTGGAGCCGGTGTTTGTCATGGGCCTGGGCTGCAGGGATCCTGCTACAA     #############AB=?:*B?;A?<2+233++;A+A2+<7==@7,A<A<=>     AS:i:-5 XN:i:0  XM:i:2  XO:i:0  XG:i:0NM:i:2   MD:Z:8A21T20    YT:Z:UU XS:A:+  NH:i:4  CC:Z:chr15      CP:i:102517592  HI:i:0

看上去很类似fastq文件,它也有read名称,序列,质量等信息,但是又不完全一样。首先,每个read只占一行,只是它被tab分成了很多列,一共有12列,分别记录了:

  1. read名称
  2. SAM标记
  3. chromosome
  4. 5’端起始位置
  5. MAPQ(mapping quality,描述比对的质量,数字越大,特异性越高)
  6. CIGAR字串,记录插入,删除,错配以及splice junctions(后剪切拼接的接头)
  7. mate名称,记录mate pair信息
  8. mate的位置
  9. 模板的长度
  10. read序列
  11. read质量
  12. 程序用标记

显然,其中chromosome至CIGAR的信息都是非常重要的。但是这些对我们不重要,我们只需要了解SAM/BAM文件是什么,就可以了。重要的是如果进行下游的操作。

要操作SAM/BAM文件,首先需要安装samtools。它的安装过程和所有的linux/unix程序一样,都是经过make之后生成可执行程序,然后把它的路径告知系统,或者放在系统可以找到的位置就可以了。
比如:

tar zxvf samtools-0.1.18.tar.bz2
cd samtools-0.1.18/
make
samtoolpath=`pwd`
PATH=PATH:$samtoolpath

然后就可以按照samtools主页上介绍的工具进行各种操作了。我们最常见的几步操作比如
0. SAM,BAM转换

samtools view -h file.bam > file.sam
samtools view -b -S file.sam > file.bam

1. sorting BAM文件。大多数下游程序都要求BAM文件是被排过序的。

samtools sort file.bam outputPrefix

2. 创建BAM index。这也是被大多数下游程序所要求。

samtools index sorted.bam

3. index模板基因组。这也是被大多数下游程序所要求。

samtools faidx Homo_sapiens_assembly19.fasta

在很多时候,我们还会看到一种扩展名为BED的mapping文件。其具体格式也是几经变化,但是现在以UCSC的描述为准。从BAM文件转换成BED文件,我们需要安装BEDtools。下载安装就不多说了。示例一个如何从BAM文件转换成BED文件的命令:

bamToBed -i reads.bam > reads.bed

。更多的具体内容可以参见其说明文档

当然,还有很多种格式来记录mapping的结果,大多数都收录在UCSC的帮助文档中。比如上次有人问及的.bw是什么文件(bigWig文件)之类的,都可以在那里找到答案。

上次谈及fastq文件时,有讲过其质量评估的问题,那么在mapping之后,如何对mapping的结果进行评估呢?
最简单的,就是通过samtools来评估mapping质量了。

samtools idxstats aln.sorted.bam

注意,这一步之前需要经过sort和index。结果会显示:

chr1	195471971	6112404	0
chr10	130694993	3933316	0
chr11	122082543	6550325	0
chr12	120129022	3876527	0
chr13	120421639	5511799	0
chr14	124902244	3949332	0
chr15	104043685	3872649	0
chr16	98207768	6038669	0
chr17	94987271	13544866	0
chr18	90702639	4739331	0
chr19	61431566	2706779	0
chr2	182113224	8517357	0
chr3	160039680	5647950	0
chr4	156508116	4880584	0
chr5	151834684	6134814	0
chr6	149736546	7955095	0
chr7	145441459	5463859	0
chr8	129401213	5216734	0
chr9	124595110	7122219	0
chrM	16299	1091260	0
chrX	171031299	3248378	0
chrY	91744698	259078	0
*	0	0	0

其中第一列是染色体名称,第二列是序列长度,第三列是mapped reads数,第四列是unmapped reads数。

如果是RNAseq,我们可以使用broad institute的RNA-SeQC来得到更加完整的报告。下载到文件之后,也许需要安装BWA来获取更精准的结果,但是如果不安装的话,也可以进行分析。一般来说,这一步不需要特别精准的结果,所以我很少使用BWA选项。下载的文件如果是.zip结尾的,直接把它改写成.jar就可以运行了。

在它的主页上下载所需要的Example RNA-seq Data。下载结束之后,该解压的解压缩。接下来运行:

samtools index example/ThousandReads.bam
samtools faidx example/Homo_sapiens_assembly19.fasta
java -Xmx2048m -jar RNA-SeQC_v1.1.7.jar -n 1000 -s "TestId|example/ThousandReads.bam|TestDesc" -t example/gencode.v7.annotation_goodContig.gtf -r example/Homo_sapiens_assembly19.fasta -o ./testReport/ -start gc -gc example/gencode.v7.gc.txt

以上的参数只有一个与其说明文档不一样的地方就是使用了-Xmx2048m来指定java虚拟机的内存大小为2G。如果遇到java.lang.OutOfMemoryError,还可以指定得再大些。
当然如果是自己的文件的话,还需要多两步:
1.BAM,reference及GTF文件的基因组名称必须一致。
2.需要使用picard工具包中的CreateSequenceDictionary来构建一个dictionary文件。

40 thoughts on “NGS分析入门:操作SAM/BAM文件

  1. Pingback: 如何统计BAM文件中的reads数 ← 糗世界

  2. Reply zzy 5月 16,2014 7:27 下午

    博主,你好!
    我想请教:有没有工具或者软件可以把SAM格式文件转变成一种描述所有被mapping的基因座位的reads counts或类似这样形式的文件?

    谢谢

  3. Reply zzy 5月 26,2014 10:37 下午

    博主您好!
    还要请教,我手中有SAM文件,有gtf注释文件,已经用bedtools genomecov转换成了各个位点的reads count形式,有没有工具可以达到以下目的:
    给具有reads count的染色体位点加上基因信息注释?
    谢谢

  4. Reply 清泉一边 6月 7,2014 3:19 上午

    博主您好:
    向您请教个问题,我有用tophat处理转录本得到的bam文件,有没有什么办法得到一个excel表格或文本呢?(表格中包含每个基因覆盖的reads数目,rpkm值,p-value值)
    谢谢您!

  5. Reply chhy 8月 14,2014 6:36 上午

    楼主,我想请问一个问题:
    请问, SAM格式中的optional elds中的TAG如果是XM和MD是代表什么意思?就是你这里说到程序标记字段
    还有,其中的POS字段中的位置,是相对于整个基因组的一开始,还是相对于当前染色体的开始位置?
    我是学计算机的,这个实在搞不清楚…看了文档之后还是难以理解..

  6. Reply zhangyan 9月 2,2014 6:53 下午

    博主您好,我有三个问题想要咨询您。
    1公司给我的bam文件据称是没有比对过的,这可能吗?那我应该用什么软件来操作这个bam文件来进行比对呢?(还有同步的fastq文件)
    2几种比对软件的联合分析的大致流程和操作是什么呀?
    3ion torrent平台的测序结果在比对操作时有什么技巧和不同吗?

  7. Reply liang 9月 18,2014 12:12 下午

    博主您好,

    我想画一个图,横坐标是基因组位置,纵坐标是rpm/bp, 请问一下samtools可以得到这样的数据吗?能给我一些提示吗?我有bam文件。谢谢博主。
    我的想画的图如链接中所示的figure1.
    http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3653129/figure/F1/

    • Reply admin 9月 20,2014 6:01 下午

      你可以试着搜一下UCSC genome browser。或者IGV, IGB之类的。如果还是不行,可以试试bioconductor的rtracklayer, Gviz, 以及trackViewer。

  8. Reply Mars 10月 14,2014 7:05 下午

    楼主,你好!
    我想知道如何从tophat的结果bam/sam文件中求得map到exon,intron的reads counts?还有从ensemble下载的注释文件,是否要进行筛选,只保留下染色体的基因注释部分?

    • Reply admin 10月 20,2014 8:52 上午

      这个可以使用Rsamtools + ChIPpeakAnno的assignChromosomeRegion来解决。我手头有一个开发版的ChIPpeakAnno,做起来比较方便,可以从http://pgfe.umassmed.edu/ChIPpeakAnno/处下载试用。

    • Reply qighu 2月 8,2017 2:54 上午

      你好,我刚接触这方面知识,想问下,你的这个问题得到解决了吗?从bam文件和gff文件求得map到exon,intron的reads counts。我想问下你的具体做法。最近比较愁。。多谢

  9. Reply huyudong 10月 29,2014 3:19 上午

    老师你好,我做的是使用Illumina通过bowtie-cufflinks获得mRNA-seq结果,我想通过测序获得的数据,分析每个gene的转录起始位点区域的碱基序列从而得到基因的启动子,请问该怎么操作呢?我尝试过将cuffmerge得到的总gtf文件导入IGV浏览器里,但是发现IGV里gtf文件都只显示CDS区域,并没有之前的5‘UTR区。
    可不可以从最初的bam文件合并得到原始的5’UTR,该怎么操作呢,谢谢

    • Reply admin 10月 30,2014 2:22 下午

      这个结果比较奇怪,通常我们都到的结果都会包含5’UTR, 3’UTR区域。
      你可以从最初的bam文件中直接提取5’UTR的reads数。这个没有问题。操作时可以有很多办法。你可以使用Bedtools来完成这个步骤。

    • Reply zhihuawu 11月 27,2014 7:10 上午

      你好,我用IGV导入gtf文件,为什么不成功啊,要提示 建立索引 ,怎么建立?谢谢

  10. Reply smrs 11月 16,2014 9:12 下午

    博主,您好。

    用bwa将短序列map到参考序列时,产生的sam文件中第二列flag会有很多数字,比如0,这代表短序列和参考序列能够匹配,而flag=4,通常情况下表示不能匹配。

    但是falg=4时,还会有很多其他结果:

    1)比如短序列t0012773 AAAGCTCAGGAGGGATAGCGCC和下述参考序列
    >ath-miR390a-3p MIMAT0031902 Arabidopsis thaliana miR390a-3p
    CGCTATCCATCCTGAGTTTCA

    匹配结果如下:

    t0012773 4 ath-miR390a-3p 21 25 22M * 0 0 AAAGCTCAGGAGGGATAGCGCC * XT:A:U NM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:22

    这儿的flag=4表示什么匹配情况呢?

    2)短序列t0031842 TCCCAAATGTAGACAAAGCATT和下述参考序列
    >ath-miR158a-3p MIMAT0000176 Arabidopsis thaliana miR158a-3p
    TCCCAAATGTAGACAAAGCA

    t0031842 4 ath-miR158a-3p 1 25 22M * 0 0 TCCCAAATGTAGACAAAGCATT * XT:A:U NM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:22

    这块实际上能匹配上,那第二列flag为什么不是0,而是4呢?

    谢谢博主。

  11. Reply YidaZhang 12月 1,2014 8:10 下午

    楼主,你好!
    我想请教一下,我用tophat做完mapping之后得到一个bam文件,我想把里面unique mapping的read和multiple mapping的read区分开来,比如说分成2个文件,有什么办法能做到吗?

  12. Reply XiaoyunLei 7月 21,2016 9:12 下午

    谢谢博主。看评论,“用bwa将短序列map到参考序列时,产生的sam文件中第二列flag会有很多数字,比如0,这代表短序列和参考序列能够匹配,而flag=4,通常情况下表示不能匹配。但是falg=4时,还会有很多其他结果。” 帮助我解决了一个疑问。 fastq匹配到参考基因组后,得到的bam文件,如果 samtools -hb -F 4 in.bam out.bam 。就可能把很多其实匹配的reads丢掉了。怪不得我原来有4亿多个reads, 用了这一步后就只剩下约5百万个reads了。

  13. Reply wangdongi 8月 2,2016 8:25 下午

    博主你好,我手里的有mapping之后得到的bam文件,现在我要用Fusim做融合基因数据的仿真,该软件需要与数据对应的gtf文件,请问我该怎么做?

  14. Reply 生信新手啊 8月 5,2016 4:43 上午

    麻烦问下 最近的项目需要从bowtie2生成的bam文件中提取read号和匹配到基因组的起始位置和终止位置。我想问有没有软件能够直接从bam里面提取read号和起始位置,然后输出。我之前用pysam读取bam,然后输出到一个txt文件中,然后发现好多chr的号都是错的,比如chr44… 谢谢

    • Reply admin 8月 8,2016 8:51 上午

      你可以直接使用samtools中的view来将bam文件转换成sam文件流,直接传给处理文本文件的程序使用。我几乎没有python编程的经验,所以对pysam出现的问题无从解答。

  15. Reply 生信新手啊 8月 8,2016 10:51 下午

    谢谢回复。发邮件问了作者。大意是bam文件很多就是没有存储chr号的,也就是输出的都是0。如果是python的话最后用的是HTSeq 读取bam并成功输出了read号和chr号和起始位置,终止位置。
    这个文件本来是sam的,通过samtools转换的bam,为了方便读取。 博主这里提到sam文件流直接传给处理文本文件的程序使用,能不能稍微详细一点,我对这个不是很懂,谢谢。
    我只知道samtools view能够直接看sam文件,但是“sam文件流直接传给处理文本文件的程序”是相当于在python之类的环境下用samtools打开sam文件,然后将数据传给下游程序,这个意思吗,谢谢指点。

    • Reply admin 8月 11,2016 8:17 上午

      是这个意思。流操作只是为了减少硬盘占用,减少操作时间。你也可以写成一个sam文件,然后用自己舒服的方式来操作。

  16. Reply 狗蛋 3月 24,2017 3:04 上午

    博主你好:
    我想问下如何找出sam文件中的这些pair-reads,其中一条read匹配成功,另一条匹配失败或者说是只有部分匹配上了。

Leave a Reply

  

  

  

%d 博主赞过: