NGS分析入门:mapping reads到基因组 61

首先,我们需要明确一条,那就是是否有reference genome(参考基因组)。如果没有,那么你应该使用的工具是de novo assembly tools,比如说

另一种情况就是已经有参考基因组。这种情况最为常见。对于基因组,我们需要有一个基本的常识,那就是随着时间的发展,对于同一物种,其基因组会时常更新,比如人类基因组现在就有:2003年六月的NCBI34/hg16, 2004年五月的NCBI35/hg17,2006年三月的NCBI36/hg18,以及现在的2009年二月的NCBI37/hg19等多个版本。现在最为常见的是hg18及hg19两个版本。所以,当我们拿到一个mapping好的文件时,必须知道它是map到哪个版本的基因组上的。
下面列举一些常用的基因组版本名称:

hg19	Human genome
hg18	Human genome
hg17	Human genome
mm10	Mouse genome
mm9	Mouse genome
mm8	Mouse genome
rn5	Rat genome
rn4	Rat genome
dm3	Fly genome
ce10	Worm genome
ce6	Worm genome
danRer7	Zebrafish genome
tair10	Arabidopsis genome
sacCer3	Yeast genome
sacCer2	Yeast genome
galGal4	Chicken genome
xenTro3	Frog genome
xenTro2	Frog genome

有一种测序是用于检查基因组修饰的,比如说MethylC-Seq, BS-Seq, or RBBS-Seq,它们使用bisulfite处理了测序样品,所以它未被修饰的A会变成T,G会变成C,所以它的方法会与常规的方法有所不同,可以在ABI wiki中去了解和掌握其分析方法。
nmeth.1828-F1

这里所需要强调的是最为常规的map方式。常用的工具有bowtie/bowtie2, BWA, SOAP1/SOAP2等。这个问题又会被分成两个问题,是基因组测序(DNA-seq)还是转录组测序(mRNA-seq)。其中的区别是对于真核生物而言,mRNA序列与DNA序列并不完全相同,在经历了后剪切之后,成熟的mRNA可能是原基因的一部分,甚至顺序及个别碱基会产生变化。如果是mRNA测序,那map工作就会在DNA测序map的基础上再多一步,map到转录组上去。所以最为流行的做法是,使用bowtie来map DNA测序,使用tophat来map RNA测序。实际上,tophat是通过调用bowtie来完成工作的。bowtie有1和2的差别。bowtie1出现的早,所以对于测序长度在50bp以下的序列效果不错,而bowtie2主要针对的是长度在50bp以上的测序的。而tophat1和tophat2的差别最主要的就是调用了bowtie1还是bowtie2。当然如果你只安装了bowtie1的话,tophat2也可以调用它。

使用bowtie的基本步骤是:

  1. 生成index文件。常用的index文件可以去bowtie网站上下载。把它拷贝到需要的目录下就可以了。
  2. 比对序列。比对最常用的命令行为:
    /path-to-bowtie-programs/bowtie2 -p <# cpu> -x <genome index prefix> <fastq file>  > <output filename>
     
    i.e.
    /programs/bowtie2 -p 8 -x hg19 Experiment1.fastq > Experiment1.sam
  3. 其中,-p参数用于指定cpu数,也就是分成几个线程并行运算。-x参数用于指定index文件的前缀(第一步下载或者生成的文件)。输出格式默认为SAM。可以指定为多种其它格式,比如说BAM(SAM的二进制文件)等等。

使用tophat的基本步骤是:

  1. 理解bowtie
  2. 使用tophat来mapping reads。其命令常见的形式为:
    /path-to-bowtie-programs/tophat -o <output directory> -p <# cpu> </path-to-genome-index/prefix> <fastq file>    
     
    For example:
    /programs/tophat -o TophatOutput/ -p 8 /programs/indexes/hg19 Experiment1.fastq
     
    Paired-end Example:
    /programs/tophat -o TophatOutputPE/ -p 8 /programs/indexes/hg19 Experiment1.r1.fastq Experiment1.r2.fastq

    可以发现,其很多参数是同bowtie是一样的。但是它有几个重要参数需要了解:

    --library-type <fr-unstranded | fr-firststrand | fr-secondstrand>
      用于生成RNA-seq的library。最常见的是使用fr-unstranded,两条链都考虑。
    -G <GTF file>
      用于加注transcriptome信息。
  3. GFF文件可以去UCSC Table Browser下载。
    ucsctablesGTF

61 thoughts on “NGS分析入门:mapping reads到基因组

  1. Reply Bo 4月 13,2013 10:37 上午

    您好。我通过tophat 2已经做完了mapping。请问怎么样从output的文件里看出所有accepted reads里有多少是exon reads,有多少是junction reads (就是中间夹了intron)呢? 谢谢

    • Reply admin 4月 14,2013 9:03 下午

      关于如何查看accepted reads里有多少exon reads,可以使用samtools idxstats工具。对于有多少junction reads,可以查看tophat2的junctions.bed文件,

  2. Reply Bo 4月 15,2013 12:19 下午

    谢谢回答!关于junction reads,我不是很清楚junctions.bed中score的含义,是不是就是每一个junction 被cover的次数?所以只要把所有的score加起来就是这个sample所有junction read数的总和了?

  3. Reply Bai 6月 26,2013 11:14 下午

    您好,我已经用tophat mapping基因组生成了结果,请问我如何能够计算出在gtf,在基因间,在内涵子中的read数?

    • Reply admin 6月 27,2013 10:18 上午

      你可以先学习一下如何使用samtools和BEDTools,而后再学习一下bioconductor/DEXSeq, 或者cufflinks。这些工具可以完美地解决你的问题。

  4. Reply zhl 7月 9,2013 10:33 下午

    您好。我通过bowtie2已经做完了mapping。请问怎么样从sam文件里看出每条染色体的unique

    • Reply admin 7月 10,2013 8:22 上午

      如果你使用的是最新版本的bowtie2的话,你应该使用–very-sensitive参数以及–non-deterministic参数。这样可以保证你的比对结果尽量的准确。但是对于uniquely mapped reads, 只能通过后期过滤的办法来解决。在bowtie2的说明中是这样写的:XS:i:
      Alignment score for second-best alignment. Can be negative. Can be greater than 0 in –local mode (but not in –end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read.
      所以,你只需要把所有带有XS:i:
      标签的reads都过滤掉,留下的就是uniquely mapped reads了。

      还有一个办法就是使用-k 2参数。这样,对于有多个maping的reads,它就会出现两个结果,你把这些有两个结果的都过滤掉,留下的就是uniquely mapped reads了。

      第三个办法是不使用-k, -a等参数,这样结果中就会出现MAPQ值。MAPQ其实就是map quality的意思,它与uniqueness正相关。也就是说,如果它的得分越高,它的结果就越unique。你可以使用某个MAPQ值来过滤,比如20。但是,这一方法最不推荐。

      就两个方法而言,还是第一个方法最为准确。

  5. Reply yolanda 12月 25,2013 4:41 上午

    O(∩_∩)O圣诞节快乐~
    您好,对于tophat的结果我一般都是用IGV看的,我现在想将具体某个染色体或某个区段上的mapping的具体情况生成文件,以及每个基因是否被mapping到全长,不知道如何能够比较好的实现,可以给我指点一下吗?谢谢!

    • Reply admin 12月 26,2013 9:52 上午

      我想你可以首先使用samtools的view工具提取你所需要的某个区段上mapping的结果,然后对gff文件也进行相应的处理,再计算就好了。

  6. Reply handong 2月 25,2014 7:46 上午

    你好,用-k 2参数之后 得到有多个maping的reads,它就会出现两个结果,要怎么把这些有两个结果的都过滤掉?

    • Reply admin 2月 25,2014 8:17 上午

      对于bowtie2需要后期自己过滤。但是对于bowtie,可以通过-m 1的设置来过滤。至于bowtie2的过滤方法,我也研究了很久,没得到很好的办法。现在最可行的办法就是使用XS:i:标签来过滤。其说明文档写道:
      XS:i:

      Alignment score for the best-scoring alignment found other than the alignment reported. Can be negative. Can be greater than 0 in –local mode (but not in –end-to-end mode). Only present if the SAM record is for an aligned read and more than one alignment was found for the read. Note that, when the read is part of a concordantly-aligned pair, this score could be greater than AS:i.

      还有很多情况,请参见:http://pgfe.umassmed.edu/ou/archives/3321

  7. Reply super 2月 26,2014 6:17 上午

    博主你好,看到你这个流程我有几个问题问一下。
    我也阅读了很多文献,为什么文献基本上没谈rna-seq没有其它NGS技术之前的去接头,剪切等预处理?而是一般都是上来就mapping,请问是这样么?
    如果不是请问需要做什么预处理? 如果需要我有几个问题
    1。怎么知道需要不需要(我目前知道的就是可以用FASTQC看reads质量)?
    2。都有哪些步骤?每个步骤有什么软件值得推荐一下?(我知道去接头好像是cupadapt)
    非常感谢博主。

    • Reply admin 2月 26,2014 9:37 上午

      通常拿到测序数据有两种,一种是multiplex的数据(就是一条lane里有多个样品,后期通过识别barcode来区别不同的样品),一种是一条lane只有一个样品的。后者不需要去接头,所以直接就可以拿来mapping。而前者,如果你提供好barcode的话,测序软件会自动demultiplex出来数据,这些软件是买测序仪时就提供的,所以通常拿到手的就是demultiplex的数据。于是也可以直接拿来mapping。但是有些时候因为实验设计比较复杂,比如说两侧都有barcode这种,就需要自己写代码来demultiplex了。这样的软件,我是使用自己写的代码,因为我还是希望保留barcode中有一个错配碱基的(如果这个碱基的评分不高的话)。所以没有其它软件可以推荐了。之前需要fastqc来看质量,应该是比较流行的做法。

  8. Reply super 2月 26,2014 6:43 上午

    补充一下,上次那个linux bash操作的问题我自己解决了。按鼠标中间的滑轮可以复制粘贴。抱歉是bionformatics的新手给您添麻烦了。

  9. Reply super 2月 26,2014 3:03 下午

    那么FASTQC查看质量不好的(比如小于20)怎么剔除?
    另外需要不需要剔除很短的序列?或者N多的?有什么软件么谢谢。
    因为发现文献大部分都是讨论DE比较的,估计觉得预处理太naive了吧
    可我是一窍不通,马上要拿到动物病毒的数据了,proton和illunima的都有。

  10. Reply super 3月 2,2014 3:27 下午

    博主你好,上面一个问题您还没有回答我“那么FASTQC查看质量不好的(比如小于20)怎么剔除?另外需要不需要剔除很短的序列?或者N多的?有什么软件么谢谢。”。
    是不是我问的太多了疏漏了没看到?辛苦!
    因为我觉得预处理肯定也是RNA-SEQ里很重要的一块但基本上论文和论坛都不说它(或许觉得它太naive了)

  11. Reply zluna 7月 9,2014 9:18 下午

    本人新手……
    1.我从ncbi下载了一些chip-seq数据,sra转成fastq了,这一步有什么特殊要求,需要在转换的时候加分割成spot参数吗?
    2.对于有些fastq文件的reads是35bp,有的比较长,是否要考虑去掉接头等信息?还有就是trim掉低质量的reads除了fastx还有更快的吗?
    3.用bowtie2还是bowtie来mapping呢,官方说50bp一下的可以用bowtie来做,请问您什么建议吗?
    4.mapping完的文件怎么获得只有unique reads的结果呢?chip-seq一般都要求unique reads的
    谢谢

    • Reply admin 7月 9,2014 9:47 下午

      1. 需要看具体情况。首先要明白什么是spot,其实就是测序的一个簇。分不分割与具体上传的文件相关。如果一行包括了两条(比如pair-end),就需要。
      2. 一定要去除接头。如果你觉得fastx不够快,你可以自己写程序。
      3. 我建议使用bowtie,因为它可以设置-m 1参数,这样就回答了你的问题4.

  12. Reply zluna 7月 9,2014 9:59 下午

    感谢及时回复
    1.对于测序文件,都是single read。还有我比较疑惑地方是,每个样品说的是测了两个lane,这个样品有两个sra文件,转出来后是两个fastq,需要用cat合并他们吗?
    2.去接头这个我真不会,目前已知的是文件的测序平台等信息,如:
    http://www.ncbi.nlm.nih.gov/sra/SRX015119
    目前没弄明白应该怎么去

  13. Reply Phate 7月 12,2014 12:52 下午

    博主你好
    想問一下,使用bwa和bowtie2所產生的sam檔,兩者之間有什麼差異?

    • Reply admin 7月 14,2014 2:13 下午

      差异有两个方面,一个是文件的tags不同(这是次要的),一个是比对算法不同而导致的比对结果不同(这是主要的)。关于后者,请查看bwa作者的分析网页:http://lh3lh3.users.sourceforge.net/alnROC.shtml。关于前者:bowtie2特有的标签有:YS, YF, YT, bwa特有的标签有:BC, X0, X1, XT, XA, XF, XE。其中,XM, NM, XS和XT比较重要。

  14. Reply 清泉一边 7月 16,2014 6:49 上午

    请教博主,有没有做过甲基化分析呢?可否给我介绍一个流程及其运用的一些软件呢?谢谢博主!!!

  15. Reply smrs 11月 22,2014 4:16 上午

    博主,你在文章提到“所以最为流行的做法是,使用bowtie来map DNA测序,使用tophat来map RNA测序”,那bwa算法呢?如果现在有一边miRNA的solexa测序数据,想把短序列map到成熟体miRNA序列上去,这个用什么算法比较好?谢谢博主。

    • Reply admin 11月 24,2014 11:43 上午

      miRNA这个可以使用bwa.但是bwa主要还是应用于DNA测序。而且我一般都使用mirdeep来处理miRNA-seq的结果。我提汲的bowtie及tophat它们是一个系列的。而bwa是另外一下系统。最初Li heng写bwa的时候是想着高速准确地map short reads至genome。他本人一开始也反对使用bwa来map RNA-seq的结果。但后来他有些松动,具体做了什么改变,我不太清楚,但我习惯了使用bowtie,就没有试bwa。只是在chip-seq, iCLIP-seq的时候会用bwa。

  16. Reply smrs 11月 22,2014 4:34 上午

    博主,您好,看了您的很多博客,收益非浅。想继续请教您一个问题。
    miRNA的solexa测序数据也是属于RNA-seq不?那这个需要使用tophat工具不?谢谢。

  17. Reply smrs 12月 1,2014 10:43 下午

    非常感谢博主的回复,博主我还有一耳光问题想请教您,用bwa 将短序列 >miR-466g 匹配 到参考序列 >ath-miR8176 时,结果的第二列是20,这种结果是什么匹配情况呢?

    短序列:
    >miR-466g
    ATAAATGACACATCATCACACA
    参考序列:
    >ath-miR8176 MIMAT0032775 Arabidopsis thaliana miR8176
    GGCCGGTGGTCGCGAGAGGGA

    匹配结果:
    miR-466g 20 ath-miR8176 21 25 22M * 0 0 TGTGTGATGATGTGTCATTTAT * XT:A:U NM:i:1 X0:i:1 X1:i:0 XM:i:1 XO:i:0 XG:i:0 MD:Z:0A21

    这块的20应该是16+4的结果,因为结果中给出的是短序列的反向互补序列,但是即使是互补序列,它和参考序列也没有成功匹配呀。

    • Reply admin 12月 3,2014 2:34 下午

      依据li heng的说法,16只代表它map到了reverse strand上,它并不代表它是否map到了参考序列上。http://seqanswers.com/forums/showthread.php?t=5347

  18. Reply strangephone 5月 27,2015 8:38 上午

    老师您好,看您的博客收益颇多,目前也在进行NGS相关解析,想请教两个问题:
    1.在进行RNA-seq mapping的时候,想对比一下用bowtie2和tophat2分别mapping之后会有多少read mapped上的数量上的差别,于是做了两遍mapping,想请问是否有更便利的方式进行这两个软件的对比。
    2.在下载reference genome文件的时候,对应bowtie2软件的参考序列后缀为.bt2, 而对应tophat2软件的参考序列后缀为.ewbt,想请问同一物种的这两种参考序列是否一致,这两种后缀有什么区别。

    谢谢您的阅读

    • Reply admin 5月 28,2015 2:17 下午

      谢谢你对我博客的支持。
      关于问题1, 你需要阅读一下tophat的原始文献。简单地讲,tophat在后台调用了bowtie。bowtie是针对genomic DNA设计的,而tophat做的就是把transcriptome拿出来做为模板,让bowtie把它当成DNA模板来比对,如果比对不上的,再直接map到基因组中。
      关于问题2,我想你可以试着自己做一下bowtie的index文件,你就了解了什么时候会出现.ewbt,bt2之类的。

  19. Reply super 7月 10,2015 6:01 上午

    请教欧老师 如果做de novo assembly的rna-seq 有很多工具 比如trinity等 如何评价它结果的可信性 就是说我怎么知道它组装的正确率等?我不记得之前是不是问过这个问题了

  20. Reply Liping 7月 29,2015 5:51 下午

    Hi,博主好
    我用bowtie分析small RNA seq的数据。设置-m 1000
    bowtie 的log文件为
    # reads processed: 16669984
    # reads with at least one reported alignment: 7921340 (47.52%)
    # reads that failed to align: 8748644 (52.48%)
    Reported 7921340 alignments to 1 output stream(s)
    # reads processed: 8748644
    # reads with at least one reported alignment: 3686286 (42.14%)
    # reads that failed to align: 5062333 (57.86%)
    # reads with alignments suppressed due to -m: 25 (0.00%)
    Reported 10264957 alignments to 1 output stream(s)

    请问
    (1)# reads processed: 8748644是表示第一次的unaligned reads 8748644 (52.48%)又被重新align了一遍吗?
    如果是,第二次align与第一次align有什么区别呢?unaligned reads到底占多少比例呢?
    如果不是,那# reads processed: 8748644表示什么意思呢?
    (2)# reads with alignments suppressed due to -m是什么意思呢?
    (3)Reported 10264957 alignments to 1 output stream(s) 的来源是什么啊? 7921340 与 3686286之和大于 10264957?

    谢谢博主

  21. Reply Liping 7月 30,2015 12:20 下午

    博主请自动忽略我的上一个问题,我犯傻了。
    我做了两次bowtie,第一次过滤掉了tRNA, rRNA等。后来用过滤掉的reads又与genome做了align。

    • Reply admin 7月 30,2015 12:23 下午

      我也说奇怪呢,我还从来没有遇到过你前文提到的情况。我还在想你是不是把多个任务的log文件写在一起了。呵呵。

  22. Reply twelvesummer 4月 14,2016 11:04 下午

    老师你好,我想问一下 GRCh38 和 GRCh37都是hg19版本吗?我拿到一个参考基因组的fasta文件,通过哪里可以查看它是什么版本呢?

    • Reply admin 4月 15,2016 11:36 上午

      GRCh38是hg38。GRCh37是hg19。如果你拿到的是一个fasta文件,你可以去http://hgdownload.cse.ucsc.edu/goldenPath/hg19/liftOver/下载hg19ToHg38.over.chain.gz文件,了解其中的不同,然后再通过查找想应的序列来确认版本。或者你可以通过已知在hg19及gh38上位置不同的基因去map至这个fa文件,来确认其版本。

  23. Reply Jack 5月 21,2016 5:56 下午

    您好,对于rna-seq,很多文章关注与mapped reads, 那对于unmapped reads呢?如果我想知道这些unmapped reads是否有有意思的结果,是不是一定要想de novo 一样拼接然后在map?然后在用些软件来检验是否是新基因,是否有表达差异?因为我现在做的是一个作物(有参考基因组和注释)的野生材料,所以想是否有必要关注一下unmapped reads,看看是否有新东西?可能很大一部分assembled contigs是因为当初map基因组时候筛选标准严格所以过滤掉了。 请问,您怎么看?如果可行的话,推荐几个软件,谢谢!

  24. Reply 谭普文 5月 30,2016 9:17 下午

    博主您好!看了您的博客获益匪浅,我想请教一下比对率低的问题,我用RSEM(调用bowtie2)比对双端测序的数据,用ensembl的人类作为参考基因组,发现overall比对率只有50%,这是一个很低的比对率,其中有10%比对了刚好一次,37.6%左右比对了>1次,还有大量的reads没有比对上,请问一下这有可能是哪些原因造成的呢?

    • Reply admin 5月 31,2016 8:50 上午

      原因有很多。我遇到的有,1. 参考基因组设置犯错误,比如应该是小鼠的,我设置成了人的。2. reads 5’端的碱基存在质量问题,3’端的碱基存在质量问题。可以先有fastqc来检查一下5’和3’端的质量,GC比,然后依照quality来trim一下reads再比对。3. barcode及adapter去除不净。使用cutadaptor及相关软件去除barcode。4. 测序质量太差。解决办法,允许更多的错配。这一步不推荐,因为结果不可靠,后期会浪费很多时间。5. 污染。有些实验人员同时在做分子克隆,污染了载体的RNA。有些是除tRNA不净,但这不应该影响比对的数量。

  25. Reply xuertowong 8月 14,2016 8:15 下午

    博主感谢您的分享!
    目前我遇到一个问题就是想通过mapping寻找到3个以上错配的相似序列。但是bowtie最多只能找到3个错配的。请问您还有别的什么软件可以实现吗?

    • Reply admin 8月 15,2016 3:47 下午

      你是指的bowtie2吗?因为bowtie可以直接通过-v来控制。在bowtie2中,你可以通过打开-N 1来允许在seed中有错配,然后再通过XM标签来过滤结果。

      • Reply xuertowong 8月 16,2016 2:20 上午

        感谢博主的回复!我用的是bowtie1 -v这个参数的范围是0-3 所以只能找出3个错配的序列。因为序列只有15bp所以没有使用bowtie2。现在CHIP抓不到片段所以希望通过信息学来预测一下,再用EMSA来鉴定。所以想找到一个程序来预测。不知道博主在这方面有没有相关的信息

      • Reply xuertowong 8月 16,2016 9:35 上午

        其实也不算是mapping吧。 就是现在有footprinting获得的序列,想通过这个序列来寻找其他相似的序列 在做EMSA验证

        • Reply admin 8月 16,2016 9:59 上午

          你可以试一下先用bowtie mapping,没map上的再用GSNAP map http://research-pub.gene.com/gmap/. 但是我没有明白你的为什么会有这么多的错配,你用紫外照过了?你是研究蛋白质和DNA/RNA相互作用的?看你的描述应该是专门研究protein-DNA相互作用的。你是想得到motif还是想得到targeting对象?

          • Reply xuertowong 8月 16,2016 11:09 下午

            博主应该是我没表达清楚。是这样,开始我想通过CHIP-sep获得一个motif。但是CHIP的数据总是出问题,因此我想预测一下在基因组上哪些序列能够和我这个蛋白结合。目前我有一个footprinting获得的DNA序列。希望通过这个序列能够找到一些相似的序列在进行验证他们是否结合。我在一片文献看到有通过信息学在全基因组范围搜索相似片段,但并没有介绍用什么软件。

          • Reply admin 8月 18,2016 1:09 下午

            ChIP-seq你都试着哪些办法,为啥认为它不work呢?是不是call不出peak?你还可以搜索是否有同源的已知的motif。如果有的话,可以使用MEME中的FIMO之类的工具搜一下。如果你说的是你只有一条或者几条DNA序列,你想在全基因组内找出相似性高的序列,你可以直接进行全基因组比对,可以使用UCSC的blat工具。如果这几条序列并不是很长,是很小的片段(只有十几个或者更短的碱基),你可以试着做一个简并的序列,然后使用FIMO搜一下。

          • Reply xuertowong 8月 29,2016 1:21 上午

            非常感谢老师的帮助

  26. Reply zhangqingyang 6月 28,2017 8:21 下午

    老师,您好,我用trinity对13套转录组数据进行拼接,但是由于服务器内存有限,老师跑了两天之后java就显示错误了,原因是线程不够吗,请问一下您有什么建议么

  27. Reply YY 7月 9,2017 10:44 下午

    博主您好,请问一下,我用一个自己binning出来的细菌contigs fasta文件作为reference,然后用bowtie2将一堆reads map到reference上去生成了sam文件。然后,我用samtool将sam转成bam,并且设置了参数将匹配上的reads过滤出来再进行sort生成新的bam文件。
    现在我想要取出这些匹配上的reads,重新进行组装。请问我该怎么做呢?我如何能将align上的reads从bam文件中拿出来呢?谢谢老师~

Leave a Reply

  

  

  

%d 博主赞过: