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

当测序得到的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文件。

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

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

    谢谢

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

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

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

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

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

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

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

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

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

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

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

      1. 对于gtf文件不应该啊。你先确认一下基因组是否正确,gtf文件是否异常。然后检查内存是否足够。

        1. 谢谢你的回答,我的gtf在linux下可以正常在IGV下打开,但是windows不行,是内存的原因?但我的内存也有8g啊,期望你的回答

  8. 博主,您好。

    用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呢?

    谢谢博主。

    1. 在http://bio-bwa.sourceforge.net/bwa.shtml中对bwa的flag有明确的描述。对于你这种情况,bwa的描述是这样的,如果short reads map的区域长度超过了contigs,那么会被标记为unmapped。比如说chromosome长度为20,但是你mapped的位置是16,而你的short reads的长度是6,就会有2个碱基在contigs之外,所以会标定为unmapped.

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

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

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

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

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

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

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

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

  15. 博主你好,我想找到rna-seq上每个read相对于参考基因组的SNP信息,目前是通过解析SAM文件的cigar列,但其匹配描述“M”包含了错配。也试了pysam和Samtools mpileup,但都无法准确定位到read上相对应参考碱基的某一碱基,也就无法判断单个read上的SNP。请问博主有什么建议吗,谢谢!

    1. 说的好复杂。没有太看懂。”M’有可能是错配,也可能是正确的。如果你pileup不出来的话,先用IGV看一眼是不是能看到先。

发表评论

电子邮件地址不会被公开。 必填项已用*标注