使用TopHat分析RNA-Seq结果 174

TopHat是一个基于Bowtie的RNA-Seq数据分析工具。它可以快速确认exon-exon剪切拼接事件。TopHat有Linux和OS X x86_64编译版本,当然也可以使用原代码编译适合自己操作系统的版本。

其上游软件是Bowtie,下游是Cufflinks。

理论上,TopHat是针对Illumina Genome Analyzer而设计的软件,它偶尔也能对其它来源的数据进行分析,但不保证成功。它针对75bp以上长度的短序进行了优化。

在使用TopHat前,必须将Bowtie的可执行文件的目录输出到PATH变量中去,例:

export PATH=$PATH:/share/sbin/bowtie

确保TopHat可以运行bowtie, bowtie-inspect以及bowtie-build。

还需要下载安装samtools

TopHat的使用范例:

tophat [options]* <ebwt_base> <reads1_1[,…,readsN_1]> [reads1_2,…readsN_2]

下载测试文件,并解压。

tar zxvf test_data.tar.gz
cd test_data
tophat -r 20 test_ref reads_1.fq reads_2.fq

如果成功的话,可以看到类似下面的输出:

[Mon May  4 11:07:23 2009] Beginning TopHat run (v1.1.1)
-----------------------------------------------
[Mon May  4 11:07:23 2009] Preparing output location ./tophat_out/
[Mon May  4 11:07:23 2009] Checking for Bowtie index files
[Mon May  4 11:07:23 2009] Checking for reference FASTA file
[Mon May  4 11:07:23 2009] Checking for Bowtie
	Bowtie version:		 0.9.9.1
[Mon May  4 11:07:23 2009] Checking reads
	seed length:	 75bp
	format:		 fastq
	quality scale:	 phred
	Splitting reads into 3 segments
[Mon May  4 11:07:23 2009] Mapping reads against test_ref with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against test_ref with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against test_ref with Bowtie
	Splitting reads into 3 segments
[Mon May  4 11:07:24 2009] Mapping reads against test_ref with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against test_ref with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against test_ref with Bowtie
[Mon May  4 11:07:24 2009] Searching for junctions via coverage islands
[Mon May  4 11:07:24 2009] Searching for junctions via mate-pair closures
[Mon May  4 11:07:24 2009] Retrieving sequences for splices
[Mon May  4 11:07:24 2009] Indexing splices
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Joining segment hits
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Mapping reads against segment_juncs with Bowtie
[Mon May  4 11:07:24 2009] Joining segment hits
[Mon May  4 11:07:24 2009] Reporting output tracks
-----------------------------------------------
Run complete [00:00:00 elapsed]

通过运行测试文件,我们知道命令当中的ebwt_base其实就是biowtie所需要的INDEXES的文件名的前缀。这些文件,我们可以从TopHat网站上下载

而之后的reads1_1[,…,readsN_1]以及[reads1_2,…readsN_2]就是需要分析的fastq文件了。如果是成对比对的话,那么文件名必须按照*_1, *_2就样成对出现。

而对于其它参数,解释其中一部分:

-o/--output-dir <string> 输出目录。默认值为 “./tophat_out”.
-r/--mate-inner-dist <int> 比对时两成对引物间的距离中值。比如说,如果你的插入片段平均有300bp,包括2个adapter各10bp, 2个barcode各6bp, 你的reads(读序长度为50bp),那么r值就应该是160=300-2*50-2*20(此为barcode在reads当中的情况),如果读序长度为100时,r值就为60=300-2*100-2*20。没有默认值,如果是末端配对比对时这个值是必须的。
--mate-std-dev <int> 末端配对时中间插入片段的长度的标准差,默认值为20bp
-a/--min-anchor-length <int> 锚定点长度”anchor length”. TopHat可以判断junction(剪切拼接)。这需要设定锚定点的最短长度,最短不能少于3,默认值为8
-m/--splice-mismatches <int> 锚定点范围内错配的个数。默认值为0
-i/--min-intron-length <int> 最短的内含子长度。默认值为70
-I/--max-intron-length <int> 最长的内含子长度。默认值为500000.
--max-insertion-length <int> 比对时插入错配最长的长度,默认值为3.
--max-deletion-length <int> 比对时缺失的最长长度,默认值为3.
--solexa-quals Use the Solexa scale for quality values in FASTQ files.
--solexa1.3-quals As of the Illumina GA pipeline version 1.3, quality scores are encoded in Phred-scaled base-64. Use this option for FASTQ files from pipeline 1.3 or later.
-Q/--quals Separate quality value files – colorspace read files (CSFASTA) come with separate qual files.
--integer-quals Quality values are space-delimited integer values, this becomes default when you specify -C/–color.
-C/--color Colorspace reads, note that it uses a colorspace bowtie index and requires Bowtie 0.12.6 or higher.
Common usage: tophat –color –quals [other options]* <colorspace_index_base> <reads1_1[,…,readsN_1]> [reads1_2,…readsN_2] <quals1_1[,…,qualsN_1]> [quals1_2,…qualsN_2]
-F/--min-isoform-fraction <0.0-1.0> TopHat filters out junctions supported by too few alignments. Suppose a junction spanning two exons, is supported by S reads. Let the average depth of coverage of exon A be D, and assume that it is higher than B. If S / D is less than the minimum isoform fraction, the junction is not reported. A value of zero disables the filter. The default is 0.15.
-p/--num-threads <int> 线程数,默认值为单线程1.
-g/--max-multihits <int> Instructs TopHat to allow up to this many alignments to the reference for a given read, and suppresses all alignments for reads with more than this many alignments. The default is 20 for read mapping (and it uses two time larger number (40) for segment mapping).
--no-closure-search Disables the mate pair closure-based search for junctions. Currently, has no effect – closure search is off by default.
--closure-search Enables the mate pair closure-based search for junctions. Closure-based search should only be used when the expected inner distance between mates is small (<= 50bp)
--no-coverage-search Disables the coverage based search for junctions.
--coverage-search Enables the coverage based search for junctions. Use when coverage search is disabled by default (such as for reads 75bp or longer), for maximum sensitivity.
--microexon-search 使用这一选项,比对时搜索microexons。这个选项只对长于50bp的短序片段起作用。
--butterfly-search 更加敏感的比对方式。如果你的测序结果是pre-mRNA样品的话,最好开启这一选项。
--library-type TopHat will treat the reads as strand specific. Every read alignment will have an XS attribute tag. Consider supplying library type options below to select the correct RNA-seq protocol.

 

174 thoughts on “使用TopHat分析RNA-Seq结果

  1. Reply Joe 3月 19,2012 6:29 上午

    Hi 翻墙过来问问题真是不易。。。
    有两个问题想请教您, 关于RNA-seq的,
    1,我下载的hg19 解压后是6个文件hg19.1-4.ebwt,hg19.rev1-2.ebwt
    如何生成hg19.fa呢?我看别人的protocol中用的都是。fa文件
    2,在cuffmerge中,显示hg19不纯在
    Error: file/directory hg19. does not exist!
    [FAILED]
    Error: could not execute cuffcompare 是不是上述的问题?
    3,gtf文件如何获得呢?

    谢谢!

    • Reply admin 3月 19,2012 8:27 上午

      @Joe, 首先,感谢您在这里留言。关于翻墙的问题,我也无奈,因为我租用服务器的时候,还没有被墙,谁知不到一年,就被墙了,而我的租约还有好几年呢……

      关于问题1,你可以先去http://tophat.cbcb.umd.edu/tutorial.html页面当中的Testing the installation当中下载测试包,打开后,你就可以看到它的测试文件,从而得知它的索引文件参数只需要一个文件前缀,可以是多个由相同前缀以及序列后缀的文件组合。这不需要转换成fa文件。

      关于问题2,的确是需要一个fa文件。而这个fa文件,你可以去http://hgdownload.cse.ucsc.edu/downloads.html下载。而问题3当中的gtf文件,也一样可以在ucsc的这个页面去找到。如果找不到,可以先阅读http://genome.ucsc.edu/FAQ/FAQformat.html页面当中的说明,然后再去前面的下载页面去找。如果实在找不到,可以在http://genome.ucsc.edu/cgi-bin/hgTables?command=start试试。如果还是不行,你可以自己生成一个。而http://cufflinks.cbcb.umd.edu/gff.html是cufflink当中自己对gtf文件的说明,你可以了解一下。

      如果还有任何问题,欢迎留言。

  2. Reply Joey 12月 22,2012 3:03 上午

    Hi, 我是tophat的新手, 我试了用http://www.broadinstitute.org/software/scripture/Walkthrough_example,发现了好几个问题,

    1. 首先我发现完全按照实例已经不可能,现在那些数据fastq已经不存在了,所以我上NCBI下载了同样的SRA的测序结果然后转换成fastq,一共是3个文件,各自大概8G左右.

    2.其次,当然tophat我用的2,bowtie也是2

    3. 运行tophat2 结果出现了accpeted_hits.bam, 但是实例里面提到的accepted_hits.sam,而且这个文件只有不到1M大, 宁外我注意到有个unmapped.bam文件大概有2G多。

    这样的结果肯定是预期的不一样的,但是从运行的过程上我没有看出有什么异样,请问大神能指点一下迷津吗?

    • Reply admin 12月 23,2012 8:19 上午

      不应该是六个fastq文件呢?成对出现?

      最好把命令行传给我,要不我无法重复你的结果。从如何得到fastq文件开始。

  3. Reply Joey 12月 23,2012 8:06 下午

    原来的6个文件的链接下载不了了,好像是NCBI改变存储的格式。所以我上http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM521650 的ftp下载到了结果的3个SRA的文件,然后有sratoolkit.2.2.2把他们转换成了fastq的格式,然后上bowtie的网站下载了mm9的index文件夹(这个是新的适用于bowtie2的bt2格式,原来的ebtw只能用bowtie1),

    然后就是执行命令了
    tophat2 -o tophat_out_SRR040001 ~/Desktop/mm9/mm9 SRR040001.fastq 我用的是tophat2

    整个执行的过程也不见的有什么问题

    2012-12-21 23:49:33] Beginning TopHat run (v2.0.6)
    ———————————————–
    [2012-12-21 23:49:33] Checking for Bowtie
    Bowtie version: 2.0.2.0
    [2012-12-21 23:49:33] Checking for Samtools
    Samtools version: 0.1.18.0
    [2012-12-21 23:49:33] Checking for Bowtie index files
    [2012-12-21 23:49:33] Checking for reference FASTA file
    Warning: Could not find FASTA ~/Desktop/mm9/mm9.fa
    [2012-12-21 23:49:33] Reconstituting reference FASTA file from Bowtie index
    Executing: ~/Desktop/RNA/base_dir/bin/bowtie2-inspect ~/Desktop/mm9/mm9 > tophat_out_SRR040001/tmp/mm9.fa
    [2012-12-21 23:51:19] Generating SAM header for ~/Desktop/mm9/mm9
    format: fastq
    quality scale: phred33 (default)
    [2012-12-21 23:51:24] Preparing reads
    left reads: min. length=152, max. length=152, 19000961 kept reads (160193 discarded)
    [2012-12-22 00:00:48] Mapping left_kept_reads to genome mm9 with Bowtie2
    [2012-12-22 04:36:38] Mapping left_kept_reads_seg1 to genome mm9 with Bowtie2 (1/6)
    [2012-12-22 06:24:05] Mapping left_kept_reads_seg2 to genome mm9 with Bowtie2 (2/6)
    [2012-12-22 08:12:58] Mapping left_kept_reads_seg3 to genome mm9 with Bowtie2 (3/6)
    [2012-12-22 10:00:17] Mapping left_kept_reads_seg4 to genome mm9 with Bowtie2 (4/6)
    [2012-12-22 11:48:10] Mapping left_kept_reads_seg5 to genome mm9 with Bowtie2 (5/6)
    [2012-12-22 13:37:14] Mapping left_kept_reads_seg6 to genome mm9 with Bowtie2 (6/6)
    [2012-12-22 15:52:24] Searching for junctions via segment mapping
    [2012-12-22 18:04:19] Retrieving sequences for splices
    [2012-12-22 18:06:02] Indexing splices
    [2012-12-22 18:09:59] Mapping left_kept_reads_seg1 to genome segment_juncs with Bowtie2 (1/6)
    [2012-12-22 18:42:57] Mapping left_kept_reads_seg2 to genome segment_juncs with Bowtie2 (2/6)
    [2012-12-22 19:16:13] Mapping left_kept_reads_seg3 to genome segment_juncs with Bowtie2 (3/6)
    [2012-12-22 19:48:40] Mapping left_kept_reads_seg4 to genome segment_juncs with Bowtie2 (4/6)
    [2012-12-22 20:21:22] Mapping left_kept_reads_seg5 to genome segment_juncs with Bowtie2 (5/6)
    [2012-12-22 20:54:02] Mapping left_kept_reads_seg6 to genome segment_juncs with Bowtie2 (6/6)
    [2012-12-22 21:35:14] Joining segment hits
    [2012-12-23 03:44:26] Reporting output tracks
    ———————————————–
    [2012-12-23 03:55:03] Run complete: 1 days 04:05:29 elapsed

    但是最后出现的文件accepted_hits.bam只有几百K,有个unmapped.bam倒是有2G左右,这显然是没有成功,但是过程又没有出错,实在找不出问题在哪儿了。

    • Reply admin 12月 26,2012 9:37 上午

      我了解的没错的话,bowtie比较适合较短的short-reads (35mer左右),而bowtie 2比较适合比较长的(50, 100, 1000mer)
      这是原文:For reads longer than about 50 bp Bowtie 2 is generally faster, more sensitive, and uses less memory than Bowtie 1. For relatively short reads (e.g. less than 50 bp) Bowtie 1 is sometimes faster and/or more sensitive.
      If your reads are shorter than 50 bp, you might want to try both Bowtie 1 and Bowtie 2 and see which gives better results in terms of speed and sensitivity. In our experiments, Bowtie 2 is generally superior to Bowtie 1 for reads longer than 50 bp. For reads shorter than 50 bp, Bowtie 1 may or may not be preferable.

      对于你的结果,似乎得不到很好的mapping。所以我建议你:
      1,将bowtie 1 和bowtie 2分别安装在两个不同的目录,然后把相应的index文件拷贝在不同的目录中。接着使用bowtie 1和bowtie 2都跑一遍,看看结果有何不同。
      2. 从你的log中看到你的short reads 有152mer这么长。这个有点奇怪。理论上应该是75X2或者说100X2(II型),所以说你的这个数据reads长度应该是100mer才比较合理。我没有下载原数据,不知道是不是正确。所以为了确认这一点,你可以使用head ~/Desktop/mm9/mm9 SRR040001.fastq来调看一下你的这个fastq文件的头几行,看看得到的fastq文件是什么样子的。
      3. 你使用了pare-end参数为75mer。我不知道这个文件是否真的包含了正反而个reads的结果。之前我说它应该是6个文件的意思就是它应该是成对出现的,比如_1, _2.fastq。从http://www.ncbi.nlm.nih.gov/sra?term=SRX019275这里看,它似乎是把正反两个文件合并在一起了,所以可能需要你将它们分开。至于如何分成两个文件,你需要自己学习一下,如果不行,再来问我。

  4. Reply Joey 12月 23,2012 8:37 下午

    [2012-12-23 17:34:04] Beginning TopHat run (v2.0.6)
    ———————————————–
    [2012-12-23 17:34:04] Checking for Bowtie
    Bowtie 2 not found, checking for older version..
    Bowtie version: 0.12.9.0
    [2012-12-23 17:34:04] Checking for Samtools
    Samtools version: 0.1.18.0
    [2012-12-23 17:34:04] Checking for Bowtie index files
    Warning: we do not recommend to have both Bowtie1 and Bowtie2 indexes in the same directory
    the genome sequence (*.fa) may not be compatible with one of them
    [2012-12-23 17:34:04] Checking for reference FASTA file
    [2012-12-23 17:34:04] Generating SAM header for test_ref
    format: fastq
    quality scale: phred33 (default)
    [2012-12-23 17:34:04] Preparing reads
    left reads: min. length=75, max. length=75, 100 kept reads (0 discarded)
    right reads: min. length=75, max. length=75, 100 kept reads (0 discarded)
    [2012-12-23 17:34:04] Mapping left_kept_reads to genome test_ref with Bowtie
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg1 to genome test_ref with Bowtie (1/3)
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg2 to genome test_ref with Bowtie (2/3)
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg3 to genome test_ref with Bowtie (3/3)
    [2012-12-23 17:34:04] Mapping right_kept_reads to genome test_ref with Bowtie
    [2012-12-23 17:34:04] Mapping right_kept_reads_seg1 to genome test_ref with Bowtie (1/3)
    [2012-12-23 17:34:04] Mapping right_kept_reads_seg2 to genome test_ref with Bowtie (2/3)
    [2012-12-23 17:34:04] Mapping right_kept_reads_seg3 to genome test_ref with Bowtie (3/3)
    [2012-12-23 17:34:04] Searching for junctions via segment mapping
    [2012-12-23 17:34:04] Retrieving sequences for splices
    [2012-12-23 17:34:04] Indexing splices
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg1 to genome segment_juncs with Bowtie (1/3)
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg2 to genome segment_juncs with Bowtie (2/3)
    [2012-12-23 17:34:04] Mapping left_kept_reads_seg3 to genome segment_juncs with Bowtie (3/3)
    [2012-12-23 17:34:04] Joining segment hits
    [2012-12-23 17:34:05] Mapping right_kept_reads_seg1 to genome segment_juncs with Bowtie (1/3)
    [2012-12-23 17:34:05] Mapping right_kept_reads_seg2 to genome segment_juncs with Bowtie (2/3)
    [2012-12-23 17:34:05] Mapping right_kept_reads_seg3 to genome segment_juncs with Bowtie (3/3)
    [2012-12-23 17:34:05] Joining segment hits
    [2012-12-23 17:34:05] Reporting output tracks
    ———————————————–
    [2012-12-23 17:34:05] Run complete: 00:00:00 elapsed

    这是用test_data执行的结果,我试过用bowtie或者bowtie2,没有什么区别。但是跟你提供的成功的输出是有区别的。

  5. Reply Joey 12月 28,2012 12:11 上午

    谢谢你的宝贵建议,我自己后来也觉得应该是两个文件合并成了一个,因为152mer正好差不多是75mer两倍,多了两个。正在看到底怎么分开。

  6. Reply xjwang 2月 26,2013 12:51 上午

    我想问一下,可变剪切分析用tophat和cuflinks应该怎么做呢?

    • Reply admin 2月 26,2013 9:05 上午

      假设你的结果是mRNA-seq,先使用tophat把所有的short reads都map到reference genome上去,然后使用cufflinks来生成BAM文件所有的transcript的gtf文件,然后再使用cuffmerge把所有需要比较的transcript的gtf文件合并起来。最后使用cuffdiff就可以得到alternative splicing的结果了。

  7. Reply xjwang 2月 26,2013 10:08 下午

    那应该看哪个文件呢?是splicing.diff这个文件吗?但是生产的这个文件最后的significant都是不显著的。
    还有一个问题就是cuffcompare也可以得到一个gtf,这个好还是cuffmerge的好呢?我邮箱wangxuejie@genetics.ac.cn太感谢了

    • Reply admin 2月 27,2013 9:34 下午

      应该是这个文件。significant不显著就需要分析具体的原因了。cuffcompare也只是比较两个gtf文件之间的不同,所以它多用于发现新基因,比如说你的merged.gtf与你提供的已知的gtf之间的区别。在做cuffdiff的时候,还是应该使用cuffmerge的结果merged.gtf。如果你不需要发现新基因,只是在已知的gtf中看alternative splicing的结果的话,你甚至可以不做cufflinks和cuffmerge这两步。

  8. Reply xjwang 2月 28,2013 1:05 上午

    你好,我们已经用tophat跑完了结果了。得到了可变剪切位点的信息即junction.bed文件,进一步应该用什么软件来找出可变剪切,以及可变剪切的分类呢?

    • Reply admin 2月 28,2013 4:42 下午

      这个junction.bed文件,记录了junction的信息,你可以用IGV结合accepted_hits.bam来查看结果。至于如何分析,我没有亲手做过这下游的事情,所以也无法回答你。如果你有了答案,请告诉我一下。

  9. Reply xjwang 2月 28,2013 9:06 下午

    多谢啦。嘿嘿,还有问一下,你做过预测新基因吗?我按照cufflinks里面的流程跑完了,但是对于结果我还是不太清楚哪个文件可以看出来新基因是哪些。我是按照这个流程跑的http://cufflinks.cbcb.umd.edu/tutorial.html#discovery

  10. Reply xjwang 2月 28,2013 9:07 下午

    还有你方便把你邮箱告诉我吗?方便交流。我的是wangxuejie@genetics.ac.cn,给我邮箱发个邮件就行,就说是糗世界

  11. Reply Chris Kuo 3月 13,2013 4:37 上午

    您好,由於看到您所敘述的設定中在使用TopHat前,必须将Bowtie的可执行文件的目录输出到PATH变量中去,例:

    export PATH=$PATH:/share/sbin/bowtie

    确保TopHat可以运行bowtie, bowtie-inspect以及bowtie-build。

    而我在安裝trinity後指令執行,卻出現一個錯誤訊息:
    which: no bowtie in (/usr/lib64/qt-3.3/bin:/usr/local/bin:/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/sbin:/opt/mpi/bin:/opt/mono/bin:/home/motion/bin:/trinityrnaseq_r2012-05-18)
    Error, cannot find path to bowtie, which is now needed as part of Chrysalis’ read scaffolding step at Trinity.pl line 537.

    猜想會不會是跟您這邊說的有相關,由於我是新手,不清楚該如何設定解決這個問題,可否請您給與意見呢?

    感謝

    • Reply admin 3月 13,2013 8:20 上午

      从你的描述来看,就是trinity找不到bowtie安装在哪里,你只需要把自己bowtie的安装路径告诉trinity就可以了。如果你还没有安装bowtie,你需要先安装bowtie。而后使用which bowtie来查看是否可以找到你安装的bowtie。如果找不到,你需要将安装目录export到PATH变量中去。例:
      export PATH=/path/to/bowtie:$PATH
      然后再运行trinity。

  12. Reply Chris Kuo 3月 13,2013 9:43 下午

    非常感謝您詳細的回答~太感激了!未來有問題在跟你請教:)

    NGS 新手感謝你 ^^

  13. Reply wujun604 4月 2,2013 1:35 上午

    你好,我想问下如果对比较短的序列(<15)做alignment的话,tophat的参数如何设置?我有个序列是14bp的,用tophat做的话出现“Error retrieving prep_reads info”错误,但是我用bowtie2的话又可以alignment, 如果我把序列长度增加到超过15就可以用tophat了。

    • Reply admin 4月 2,2013 8:02 上午

      您可否把tophat至出错前的输出传过来我看一下。同时把你的命令行,以及fastq文件头传过来。也许可以找到问题的关键。

  14. Reply wujun604 4月 2,2013 9:20 下午

    谢谢,这么迅速就回了!:)
    序列文件是.fq文件,其中前面几行是这样:
    @01:001:0001:0795
    CACTATTTCT
    +
    OLIYCEWGQF
    @01:001:0001:0279
    ACGAGTGTCT
    +
    JLLLVEGGQF
    @01:001:0001:0288
    TCTCGTTTCT
    +
    EPXISETGQE
    @01:001:0001:0309
    TAGGTTTTCT
    +
    XYWCXEFGQE
    @01:001:0001:0357
    TGTGATTTCT
    +
    CK[BOEWGQE
    @01:001:0001:0623
    TGCGCTTTCT
    +
    PPNRCEAGQF
    @01:001:0001:0770
    TGACGTGTCT
    +
    \FIQWEGGQF
    @01:001:0001:0777
    TATATTTTCT
    +
    ACQIPEPGQF
    运行出错在Preparing reads 环节,Error retrieving prep_reads info. log里面是:
    prep_reads v2.0.8 (3864)
    —————————
    613122 out of 613122 reads have been filtered out

  15. Reply wujun604 4月 2,2013 9:23 下午

    命令行是:tophat –library-type fr-firststrand -G genes.gtf -o Test_thout genome yeastSeq0N.fq
    tophat.log里面是:

    [2013-04-03 10:19:51] Beginning TopHat run (v2.0.8)
    ———————————————–
    [2013-04-03 10:19:51] Checking for Bowtie
    Bowtie version: 2.1.0.0
    [2013-04-03 10:19:51] Checking for Samtools
    Samtools version: 0.1.19.0
    [2013-04-03 10:19:51] Checking for Bowtie index files
    [2013-04-03 10:19:51] Checking for reference FASTA file
    [2013-04-03 10:19:51] Generating SAM header for genome
    Warning: found a read < 20bp in yeastSeq0N.fq
    Warning: found a read < 20bp in yeastSeq0N.fq
    Warning: found a read < 20bp in yeastSeq0N.fq
    Warning: found a read < 20bp in yeastSeq0N.fq
    format: fastq
    quality scale: phred33 (default)
    [2013-04-03 10:19:52] Reading known junctions from GTF file
    [2013-04-03 10:19:52] Preparing reads
    [FAILED]
    Error retrieving prep_reads info.
    prep_reads.log里面是:
    prep_reads v2.0.8 (3864)
    —————————
    613122 out of 613122 reads have been filtered out

    • Reply admin 4月 3,2013 8:07 上午

      看报告的意思是因为reads的长度小于20,所以全部被filter out了。你试试 –b2-L 参数,把这个参数设置成8。具体的你可以参看我的博客http://pgfe.umassmed.edu/ou/archives/2277

  16. Reply wujun604 4月 15,2013 5:32 上午

    哎呀,好不容易一不凑巧又能翻墙过来啦。大神,我后来用bowtie2做了,然后把输出格式转换了用于cufflinks读。还有个问题,我把那个R3.0编译安装了,但是画图的时后,出现“不支持半透明输出”的问题,这个是缺少哪个库么?

    • Reply admin 4月 16,2013 8:24 上午

        告诉你一个好消息,以后上我的博客不再需要翻墙了。关于R3.0,现在看好象有很多问题。很多包都无法在R3.0下编译通过。所以我们需要静等一段时间。而且我没记错的话,R本身不支持半透明输出,但是PDF支持。所以你可以使用pdf输出文件,然后使用illustrator之类的软件来编辑就可以了。

  17. Reply wujun604 4月 16,2013 8:31 上午

    不需要翻墙果然是个好消息啊!半透明输出的问题是这样的,我另外一个及其上就可以,这个机器上不行。那个机器上(2.15)是直接sudo apt-get install r-base的,所以我以为是哪个附加库的原因呢。

  18. Reply arthur 6月 11,2013 1:41 下午

    想请教您一下,我的reference是用的human_g1k_v47通过bowtie2的index得到的,然后我在整个directory底下也把原来的fasta文件放进去了,但运行tophat的时候还是会出现整个警告。
    Warning: Could not find FASTA ~/Desktop/mm9/mm9.fa
    想请教一下这是为什么?
    还有我是从ncbi下的文件,可能需要加一下header line,用tophat您知道如何加不?
    最后不知道您能不能推荐一下几个tophat的option,我目前只加了几个常用的-o之类的,
    谢谢您了啊!

    • Reply admin 6月 17,2013 8:28 上午

      可能会导致Could not find FASTA的原因有很多,你需要查看log文件来具体分析。可能的原因网上有说,1硬盘空间不足,2tophat不是最新版本,3FASTA文件格式不正确。

      关于tophat的参数,–library-type, –mate-inner-dist可能会比较有帮助。但是这种问题都没有统一的答案,具体的还是尽可能的使用正确的参数来描述你的实验。至于比对的算法方面,tophat默认参数就很好了,除非你的具体的阳对一定要看到,你可以试着调整参数,否则意义不是很大。

  19. Reply arthur 6月 19,2013 2:54 下午

    非常感谢您的解答了,但是关于-mate-inner-dist也就是-r我还有点疑问,我现在从ncbi下的seq data,并没有告诉我template是多少bp的,那我怎么设置了,我网上搜到很多人用的是manual上的-r 200,但我觉得如果不知道模板长度是不是用默认的比较好了。

    • Reply admin 6月 21,2013 1:54 下午

      完全可以使用默认值。这个值只是在alignment的时候对结果进行优化。我们常常遇到的情况是,一个read可能map到基因组的多个位置,有了这个值,可以使成对的reads尽可能的map到正确的位置。实际上,seq中的unique mapped reads应该是大多数,只有很少量的才是有多个map点的。另外,即使map至多个位点,他们成对近距离出现的可能性也很小。所以这一参数只是对算法效率上会有更多的帮助,对结果的影响其实并不大。

  20. Reply arthur 6月 27,2013 3:00 下午

    谢谢啦!!

  21. Reply AGFU 7月 2,2013 11:48 下午

    您好!请问用tophat分析之后得到的deletion.bed, junction.bed, insertion.bed有什么快速大量的处理方法来了解整个的splice情况吗?

  22. Reply zss 10月 8,2013 2:24 上午

    你好,请问用tophat 把rna-seq map到参照基因组后,对于没有map上的序列该怎么处置它们呢?

    • Reply admin 10月 8,2013 8:08 上午

      使用genome browser看一下是不是map的结果正确。然后随机抽几个blast一下,看是否是污染。如果以上都不是,可以考虑直接抛弃。

  23. Reply zss 10月 9,2013 7:17 上午

    你好!我用IGV查看的,unmap的数据在上面一点显示都没有,应该就真是没有map上的。用blast把没有比对上的序列进行两两比对吗?blast什么样的结果说明是污染呢?初学者 谢谢!

  24. Reply TY 10月 26,2013 7:48 上午

    您好,我想问一下,我使用的是pair-end数据,tophat中的-r参数的值应该如何确定?我的fastq数据长度是102,但是我不知道中间被span掉的长度有多少?

    tophat2 ../Human_Genome/bowtie_output/genome ../Reads/S_4_read_1.fastq ../Reads/S_4_read_2.fastq | tee Document.txt

  25. Reply 唐磊 12月 14,2013 9:46 上午

    安装TopHat,Bowtie,samtools这些软件是要在linux系统下吗?window系统可以吗?bowtie-inspect以及bowtie-build这些也要下吗?

  26. Reply 郑雨 1月 15,2014 2:38 上午

    您好 ,请问怎么找可变剪接的呢,直接从junction.bed 中获取么?同时怎么识别新的可变剪接的呢。有什么比较好用的注释文件的呢?

    • Reply admin 1月 15,2014 9:39 上午

      对于alternative splicing,其实我没有太多的经验。以前我自己写过一些不成熟的代码,现在也淹没在电脑硬盘里了。我的步骤大约是,获得exon level expression, get difference, annotation by exon gff from UCSC or ensembl, confirm the alternative splicing by junction.bed. 现在比较完整好用的注释文件,我首推ensembl,然后是ncbi cds database。

  27. Reply 郑雨 1月 18,2014 3:39 上午

    十分感谢博主在这么短的时间内回答我的问题。博主,我还想问一个问题,就是在cuffdiff的结果文件中会出现多个以下类似的情况:
    SNORD97 chr11:10818592-10830582 s2 s4 OK 0 8.13043 inf #NAME? 0.0005 0.0277552 yes
    CALML6 chr1:1844811-1848733 s2 s5 OK 0 1.46656 inf -nan 5e-05 0.00428151 yes
    在另一个样本中表达量很小,这样的差异明显么?
    XLOC_003358 XLOC_003358 ITGB1 chr10:33176213-33249395 s4 s5 OK 624.777 98.7652 -2.66127 -2.5435 0.0202 0.245265 no
    或者是表达差异很明显 却校正出来为no 呢?
    还有博主知道tophat是怎么去计算这个test-stat的值的呢?

    • Reply admin 1月 19,2014 7:41 上午

      你要知道,是yes还是no是基于pvalue的,你看到显示no的p.value值明显是大于0.01的。当然你可以自己调整这个cutoff value,让它也显示为yes。这个p值的计算,你可以在其帮助文件和发表的文章中找到答案。

  28. Reply zhang Jing 1月 19,2014 11:26 上午

    您好,看了您很多的博文,收获很大。想向您请教一个问题,我最近要做两个物种转座子表达的差异分析,但是由于转座子和基因不大一样,它们的序列相似度非常高,用tophat做read-mapping时,一个reads可能会比对到很多位置,影响最后的结果,请问您有什么建议呢?

    • Reply admin 1月 19,2014 4:09 下午

      在这种情况下,我建议您考虑一下BWA。这个软件在比对的过程中精确性会更好。但是它并不是针对mRNA-seq设计的,所以如果你需要做mRNA-seq的分析,会有两种选择,1,还是选用tophat2,并传递更严格的比对参数给bowtie2,具体你可以查看tophat2中的Bowtie 2 specific options部分,以及bowtie2中相关参数的说明。2,选用bwa,但是你需要选择合适的pipeline。据我所知,现在已有一些文章在使用bwa来进行mRNA-seq的比对。

  29. Reply super 1月 19,2014 3:34 下午

    博主您好,我是一个生物信息学初学者,正在学tophat,cufflink等软件。想和您多多交流和请教。老板现在给我一个数据集是有很多reads的,但是 .fastq格式的,请问如何生成您所说的“*_1, *_2”的文件呢?谢谢您!

  30. Reply 郑雨 1月 19,2014 8:23 下午

    博主您好!
    1、cuffdiff的默认FDR为0.05 判断显著性根据p-value与q-value的大小比较。如果仅考虑p-value这个值,对结果的真实性会产生影响么?
    2、博主我想再问一个问题,因为这个显著性的问题,后来我尝试了在cufflinks这一步骤中加入-G的参数 ,发现结果识别为yes的差异表达的基因数目增多了,请问这个怎么回事的呢?
    3、我按照您之前的文献参考 对RNA-seq 的tophat的生成结果进行去duplicate再进行cufflinks等一系列操作。我的数据是双端数据,我有点疑惑就是关于samtools 的rmdup中的-s 与 -S 选项的选择。

    • Reply admin 1月 21,2014 9:27 上午

      1. cuffdiff首先会用FDR来做过滤。默认值为FDR<.05。我是比较相信cuffdiff的结果的。 2. 在cufflinks时就加入-G参数时,它会首先考虑一个reads是否归属已知的transcripts。所以计数会受到影响。 3. 你的数据是双端数据,所以在使用-s参数时,会将只有单端的数据当中重复的也去除。在使用-S参数时,所有的双端数据被会被视为单端数据,去除重复。-S会产生的差异是:对于双端数据,如果一端的测序数为100,而另一端测序数为200,显然,这不是因为seq重复或者pcr重复,它应该是有差异的。但是如果都视为单端数据的话,那么两端的数据都会被去除得只余1,1.而无法体现其中一端与另一端有差异的事实。

  31. Reply super 1月 20,2014 10:35 上午

    谢谢!我的就是single end sequence,该怎么处理呢?我的序列格式是.fastq的。谢谢!

  32. Reply super 1月 21,2014 11:26 上午

    您的意思是不是就是先要用bowtie和bowtie2生成所有的文件放在一个文件夹下才可以做tophat?谢谢

  33. Reply super 2月 11,2014 1:09 下午

    斑竹你好
    我用的你说的这篇经典文章(http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html
    )和它的数据集试玩tophat,在他说的网站下载了数据,运行正常,用的tophat 和bowtie 的index。语句
    tophat -p 8 – G genes.gtf -o Output genome GSM794483_C1_R1_1.fq GSM794483_C1_R1_2.fq
    但后来今天我又换了另一台服务器上下载数据并运行,同样的语句却给我了一个错误提示,说找不到bowtie 2的index,我考虑是不是这台机器的bowtie版本不一样,于是拷贝了bowtie2文件夹里的所有index 放在运行路径下(当然没有混上bowtie1的索引),结果又来个warning: could not find FASTA file genome.fa(可明明就在路径所在文件夹里);还有一个错误,说error locating program: bowtie2-inspect。请问这是怎么回事?两台服务器的软件都不是我装的(我不是管理员没有权限)。是管理员的问题么?
    此外,新服务器的系统是XFCE UBUNTU,我了过去,bash竟然没有menu bar,也不知道咋调出来,试着用右键“显示菜单栏”但是右键却仿佛只能全选啥也做不了(可惜斑竹您不留邮箱否则我可以截图你看),所以也没法复制粘贴。真心不知道怎么回事?不知道斑竹能不能解答?
    谢谢!

  34. Reply super 2月 11,2014 1:43 下午

    你好忽略我第一个问题,我已经解决了。就是第二个问题:
    新服务器的系统是XFCE UBUNTU,我了过去,bash竟然没有menu bar,也不知道咋调出来,试着用右键“显示菜单栏”但是右键却仿佛只能全选啥也做不了(可惜斑竹您不留邮箱否则我可以截图你看),所以也没法复制粘贴。真心不知道怎么回事?不知道斑竹能不能解答?

    • Reply admin 2月 13,2014 8:43 上午

      不好意思,最近忙,没有看到您的留言。第一个问题我想多半是路径问题。还好你已经解决了。这个过程相信你学到了更多。
      bash为什么要有menu bar呢?不过你可以试一下Alt-F2 > xfce4-panel restored the panels.看看能不能出来menu bar。这好象是一个系统bug。参见:https://bugs.launchpad.net/ubuntu/+source/xfce4-panel/+bug/53897

  35. Reply list 3月 13,2014 3:56 上午

    比对时两成对引物间的距离中值。比如说,如果你的插入片段有300bp,而每个引物有50bp,那么r值就应该是200=(300+50*2)/2。没有默认值,如果是末端配对比对时这个值是必须的。

    这个你理解错了,应该300-50*2=200,插入片段已经包含reads(“引物”)本身了。

  36. Reply super 4月 9,2014 8:31 上午

    博主你好我发现您留言里提到可以不用cufflinks和cuffmerge直接拉去cuffdiff?请问能稍微详细解释一下么?
    “如果你不需要发现新基因,只是在已知的gtf中看alternative splicing的结果的话,你甚至可以不做cufflinks和cuffmerge这两步。”
    我知道用egdeR/DeSeq这些 pipeline的时候是不需要的。但texudo就要。现在概念就有点混了,谢谢

  37. Reply super 4月 10,2014 4:22 上午

    1.也就是说发现新基因这是texudo的优势之一了?不过cuffdiff的输出 比如isoform,TSS ,CDS 我真心不知道有啥用?感觉好像只有DE genes比较热门(而egerR/DESeq好像也只能处理DE Genes)。

    2.此外,我现在做的两组样本一组是被病毒感染过的动物细胞,一种是被变异后的病毒感染的。我用tophat把reads map到transcirptome上是吧?使用的是动物(老鼠)的参考基因组genome ,包括调用其bowtie2的index。对吧?
    那我需要不需要使用病毒的reference genome去做呢?有没有意义?

  38. Reply super 4月 10,2014 5:37 上午

    问题还有一个,当使用cuffdiff时会生成一系列文件。记得和您上次讨论过,比较有用的就是gene_exp.diff,然后去做下游分析(GO/pathway/network)。
    其它的呢?假设我想找相同gene or isoform的不同tss,但从tss_group_exp.diff里我也看不懂那个文件 什么TSS10000 ,TSS10004。。。让我很晕。不知道博主可否指导一下。

  39. Reply super 4月 24,2014 5:07 下午

    请问tophat也可以map到transcriptome上么?和genome上有何不同?

    • Reply admin 4月 25,2014 8:41 上午

      你要清楚,bowtie和bwa是往genome上map的工具,而tophat设计的就是为了map到transcriptome上的,所以它才用于RNA-seq的mapping工作。它会多次调用bowtie,首先map至genome,然后再将exon以外的部分以及没有map结果的reads再map至junction的位置。

  40. Reply nature 5月 18,2014 9:21 下午

    计划做RNA-edit的流程,在寻找SNV位点的过程中输出空白文件,然后在寻找原因的过程中,看到tophat mapping的tophat.log里面有这样一个问题,是在说什么事呢?能帮忙解答一下吗,谢谢!

    software/tophat-2.0.8b.Linux_x86_64/gtf_to_fasta: /lib64/libz.so.1: no version information available (required by /home/yuy/software/tophat-2.0.8b.Linux_x86_64/gtf_to_fasta)

  41. Reply FUB 5月 30,2014 7:03 上午

    您好!我用CummeRbund得到了RNA-seq的一些图,可是由于对生物信息学以及统计学基本是小白,所以不知道怎么去描述这些图片。特别是disp<-dispersionPlot(genes(cuff), dens<-csDensity(genes(cuff)), principal component analysis (PCA)所分别得到的图。他们分别的统计学和生物学的意义都是什么?dispersionPlot是表示基因表达的分布符合什么规律,然后说明统计是可信的?我的图片中各样本在dens<-csDensity(genes(cuff))中得到的是几乎都重叠在一起的曲线,比CummeRbund muanual里面更集中地重叠在一起,这又说明什么呢?最后,principal component analysis (PCA)的那些黑点点是代表所有的基因吗,处理组A和未处理组A聚在一起,处理组B和未处理组B聚在一起又有什么统计学和生物学意义呢?
    多谢赐教!

    • Reply admin 5月 30,2014 7:32 上午

      dispersionPlot图告诉你表达量与离散度的关系,这个图你应该了解录表达量多大的时候,其不同样品间的表达差异是比较可信的。当然这具体的会体现在pvalue和adjust pvalue值上。PCA告诉你哪些样品比较相近。理论上同一组的样品应该聚拢在一起。如果处理组A和未处理组A聚在一起,处理组B和未处理组B聚在一起,说明样品A与样品B之间的区别大于处理与未处理对样品总体的影响。density曲线告诉你表达值分布统计,应该符合power law法则。多个样本几乎重叠说明样品间的表达值分布接近。这对于大多数情况而言是好事。

  42. Reply FUB 5月 31,2014 1:45 下午

    您好!谢谢您的回复!关于dispersionPlot中您提到pvalue和adjust pvalue值(是最后结果里的q value?),但是在得到的图中,X轴是count,Y轴是dispersion.具体怎么来判断表达差异是否可信我还是很迷糊,可以发个图片给您看看吗?
    density曲线中的重叠表示样品间表达式分布接近,是表示FPKM归一化后,个样品比较基因的差异表达更为可信吗?

    BTW,您的新书后天可以到手了,非常期待拜读!

    • Reply admin 5月 31,2014 4:17 下午

      你以adjust pavlue值为准就好了。density曲线的问题,你的想法也是我的观点。但记住算法毕竟是算法,一切都还是以实验验证为准。当然,算法好的话,可以减少假阳性。
      BTW,所谓我的新书,我都没有见过,哈哈。

  43. Reply Ekoms 7月 28,2014 1:15 下午

    博主好,
    最近发现一个小问题:在用tophat做rna-seq的alignment,样品来源是女性,但是Y染色体上也有一点reads。。因为X和Y染色体上有一些同源基因,所以alignment的时候两边都可以,结果有些就选到Y上去了。。
    请问有没有办法把Y染色体关掉或者后期手动合并reads的方法呢?
    谢谢~

    • Reply admin 7月 28,2014 1:26 下午

      如果你对发现新基因不感兴趣的话,你可以使用非de-novo的mapping方式。将gff注释文件中所有Y染色体上的基因全部删除,然后通过-G参数传入tophat, 当然同时也要设计–transcriptome-index参数。这样你就可以得到其它染色体区域优先的mapping结果了。

  44. Reply zzbin 7月 30,2014 10:52 上午

    你好,请问我在跑
    cuffmerge -g genes.gtf -s genome.fa -p 8 assemblies.txt
    时,出现
    [00:05:31] Assembling transcripts and estimating abundances.
    Processed 11325 loci.
    [Thu Jul 31 00:05:52 2014] Comparing against reference file genes.gtf
    You are using Cufflinks v2.2.1, which is the most recent release.
    Warning: couldn’t find fasta record for ‘2LHet’!
    Warning: couldn’t find fasta record for ‘2RHet’!
    Warning: couldn’t find fasta record for ‘3LHet’!
    Warning: couldn’t find fasta record for ‘3RHet’!
    Warning: couldn’t find fasta record for ‘U’!
    Warning: couldn’t find fasta record for ‘XHet’!
    Warning: couldn’t find fasta record for ‘YHet’!
    Warning: couldn’t find fasta record for ‘dmel_mitochondrion_genome’!
    [Thu Jul 31 00:06:01 2014] Comparing against reference file genes.gtf
    You are using Cufflinks v2.2.1, which is the most recent release.
    Warning: couldn’t find fasta record for ‘2LHet’!
    Warning: couldn’t find fasta record for ‘2RHet’!
    Warning: couldn’t find fasta record for ‘3LHet’!
    Warning: couldn’t find fasta record for ‘3RHet’!
    Warning: couldn’t find fasta record for ‘U’!
    Warning: couldn’t find fasta record for ‘XHet’!
    Warning: couldn’t find fasta record for ‘YHet’!
    Warning: couldn’t find fasta record for ‘dmel_mitochondrion_genome’!
    是什么原因呢,该怎么解决?

  45. Reply zzbin 7月 30,2014 9:09 下午

    太感谢啦!!

  46. Reply zzbin 7月 31,2014 3:21 上午

    又出问题了:
    [st@localhost my_rnaseq_exp]$ cuffdiff -o diff_out -b genome.fa -p 8 -L C1,C2 -u merged_asm/merged.gtf C1_R1_thout/accepted_hits.bam,C1_R2_thout/accepted_hits.bam,C2_R1_thout/accepted_hits.bam,C2_R2_thout/accepted_hits.bam
    You are using Cufflinks v2.2.1, which is the most recent release.
    Error: cuffdiff requires at least 2 SAM files

    • Reply admin 7月 31,2014 5:50 上午

      你需要将要比较的两者分开,不要用逗号联接它们。逗号联接的是同一组内的样品。

      • Reply zhaoli 11月 4,2014 3:06 下午

        楼主您好,我运行Cuffdiff时也遇到了同样的问题
        -bash-3.2$cuffdiff -o diff_out -b genome.fa -p 8 -L C1,C2 -u merged_asm/merged.gtf\
        ./C1_R1_thout/accepted_hits.bam,./C1_R2_thout/accepted_hits.bam,./C1_R3_thout/accepted_hits.bam\
        ./C2_R1_thout/accepted_hits.bam,./C2_R3_thout/accepted_hits.bam,./C2_R2_thout/accepted_hits.bam
        -bash:you are using Cufflinks v2.2.1,which is the most recent release.
        -bash-3.2$:Error:cuffdiff requires at least 2 SAM files
        请问楼主,我应该怎么办呢,不胜感激

  47. Reply steve 7月 31,2014 9:26 下午

    你好,
    > cuff_data < – readCufflinks('diff_out')
    错误: 找不到对象'cuff_data'
    请问这个怎么解决?

  48. Reply steve 8月 3,2014 7:45 上午

    就是由cuffdiff获得不同样本的diff_out文件后,要安装cummeRbund包来可视化结果。

    我的操作步骤如下
    > source(‘http://www.bioconductor.org/biocLite.R’)
    > biocLite(‘cummeRbund’) #过程中提示有些对象被屏蔽了
    #Load the CummeRbund package into the R environment:
    > library(cummeRbund)
    #Create a CummeRbund database from the Cuffdiff output:
    > cuff_data < – readCufflinks('diff_out')
    错误: 找不到对象'cuff_data'

    我是菜鸟,刚接触生物信息,按照这篇文献做的:Differential gene and transcript expression analysis

    of RNA-seq experiments with TopHat and Cufflinks

  49. Reply Meteorick 8月 6,2014 7:27 下午

    请问一下,我的tophat2总是报错,exec format error,找不到错误啊!

  50. Reply Meteorick 8月 6,2014 9:27 下午

    博主,上面那个错误是在preparing reads之后报错啊,但是我又使用了官网上的test文件,还是一样的错。
    还有一个问题就是再运行tophat1时,checking reads之后,quality scale:–phred33-quals之后就自动退出了,输出的文件都是空的,这是为何?
    我的电脑32位fefore!希望博主可以为我解答,谢谢!

Leave a Reply

  

  

  

%d 博主赞过: