使用TopHat分析RNA-Seq结果

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. 楼主你好
    我正在做一个将不同物种在一个FQ文件中的RNA-seq reads分开的处理,想用tophat将混合的reads往物种的转录组上比对,比对上的再拿出来。但是遇到一个问题,就是tophat的结果,有accepted.bam和unmapped.bam。这两个文件是不是分别就是能比对上的和不能比对上的reads的报告呢?

  2. 博主,能给我们介绍下ASTALAVISTA这个软件,如何鉴定可变剪接?还有具体详细介绍Cufflink软件的使用呢?谢谢博主

    1. ASTALAVISTA这款软件我没有用过。不好意思。关于cufflinks,我也是初通皮毛,要我具体介绍,还得等待积累。

  3. 博主,用cummeRbund做Cuffdiff结果分析,其中做的Dispersion plot,横纵坐标的值不知道是什么意思?其中的count指的是什么?dispersion又是指的是什么?图中每个点代表什么意思?谢谢博主!

  4. 博主,您好!又有问题需要麻烦您了。
    我在用tophat分析公司给我的RNA-Seq数据时,不管是用tophat预编译的版本,还是利用源码自己编译的版本,都会在Reporting output tracks这一步出错。在reports.log文件中,是loading了所有数据后,最后一行显示的是loading done。其他log文件中也没有明确的错误信息。用官网给的test数据,没有出现任何问题。
    请问这个会是什么原因造成的呢?

  5. 博主您好!我有一个初级的问题,请问cufflinks在给transcript算FPKM的时候(bam文件来自tophat),是不是只考虑完全map在该transcript的reads? 也就是说如果一个read overlap一条transcript的5’ end 是不是不算?

    1. 这并不是一下初级问题。这个问题很复杂。你可以看一眼http://www-huber.embl.de/users/anders/HTSeq/doc/count.html网页中的图。希望你能得到解答。

  6. 楼主,你好!首先我很敬佩你!!(你的才华和能力,以及乐于助人的精神)
    我的问题如下:
    我去掉了接头和低质量序列后,原本pairend测序的fastq文件*_1,*_2,出现了不匹配的,我根据reads的id提取不匹配的放在SE.fastq中,针对SE.fastq文件,我是否可以用如下命令将其map到基因组上去:tophat -G genome.gtf -r 260 genome PE_reads_1.fa PE_reads_2.fa,SE_reads.fa (也就是SE.fasta是否可以map到基因组,若不map会有什么影响?map有会有什么影响?)

    1. 您可以尝试一下。这种处理方法我没有做过。得到结果了分享给我吧。我个人觉得应该看不配对的数量的大小,如果比例相当大,多半是会有影响的。

  7. 博主,不配对的reads的文件SE.fastq是否要转换为fasta,因为tophat的官网显示的是采用fasta格式:http://ccb.jhu.edu/software/tophat/manual.shtml

    1. %我TOphat跑完出现了大量了discordant alignments(98%).我想问一下这种情况会对结果产生影响嘛?如果产生,你有什么办法处理吗?

  8. 博主你好!看了很多留言和你的指导,还是不太懂tophat的-r(引物间的距离中值)参数的问题,pair-end测序的懂了,就是把平均的总长减去2倍的length值,那么single-end测序的怎么算呢?还有Solexa single-end测序的平均总长一般是多少?我实在查不出来,求助!!谢谢!

    1. single-end不考虑这个参数,因为它没有paired mate。一般来说,插入片段的大小由实验决定。如果是按照指导步骤来的话,通常在250-300之间。

  9. 请教纠结的水稻基因组的问题!用phytozome上Osativa_204作参,tophat-cuffDiff结果中已知基因的value全是0啊!可是Tophat再HTseq发现value是有值(Reads)!

    1. 一般都有建好的,你可以先去http://support.illumina.com/sequencing/sequencing_software/igenome.html查看有没有你要的物种,没有了再自己做。

  10. 谢谢博主!请问在命令行执行cuffmerge出问题Error: could not execute gtf_to_sam怎么回事呢?导致服务器上运行cuffmerge(sh封装)时在当前目录生成core文件,sh中断!谢谢

    1. 你好,我也遇到了您之前的问题。就是关于cuffmerge,报错是
      Error: duplicate GFF ID ‘xxx’ encountered!
      [FAILED]
      Error: could not execute gtf_to_sam

      不知道最后是怎么解决的呢?

        1. 不好意思,我第一次跑这套流程,具体也不太清楚。我把报错的转录本grep出来,大致如下:(的确出现了很多转录本,这是不是之前的cufflinks出了问题?)

          1 Cufflinks transcript 25800176 25800230 1
          1 Cufflinks exon 25800176 25800230 1 +
          1 Cufflinks transcript 25800176 25800413 1
          1 Cufflinks exon 25800176 25800413 1 +
          1 Cufflinks transcript 25800176 25818224 2
          1 Cufflinks exon 25800176 25818224 2 +
          1 Cufflinks transcript 25800231 25800231 1
          1 Cufflinks exon 25800231 25800231 1 +
          1 Cufflinks transcript 25800231 25800413 1
          1 Cufflinks exon 25800231 25800413 1
          1 Cufflinks transcript 25801043 25801160 1
          1 Cufflinks exon 25801043 25801160 1 +
          1 Cufflinks transcript 25802016 25802117 1
          1 Cufflinks exon 25802016 25802117 1 +

          1. 这是那个报错的“xxx”的结果。不过我不太会弄,后面的没能复制上来。

          2. 最后还是无法解决,只能用cuffcompare代替了merge,也没有报错,cuffdiff也正常运行

  11. 博主你好!我们跑TopHat发现reads长度最小值89,最大值101,插入片段平均长度220,请问-r参数改设多少呢?谢谢,

    1. 对于长度不一致的,可以使用最短的为标准(当然,如果最短的比例非常少,可以使用中值)。你描述的的,插入片段平均长度220,如果是片段本身220的话,220 = 2 * (reads length + barcode length) + r. 如果你的barcode不在测序片段内,或者说你的reads length已经包括了barcode的就不用再算了。于是你的r值应该是220 – 2*90 = 40 左右。这个值不用太精确。

  12. 你好,今天早上run Tophat有这个错误文件。我仔细检查了GFF3文件的和sequence文件,发现染色体名字都一致。请问你遇到过这样的问题么?谢谢!
    [2015-01-05 09:48:25] Beginning TopHat run (v2.0.13)
    ———————————————–
    [2015-01-05 09:48:25] Checking for Bowtie
    Bowtie version: 2.2.4.0
    [2015-01-05 09:48:26] Checking for Bowtie index files (genome)..
    [2015-01-05 09:48:26] Checking for reference FASTA file
    [2015-01-05 09:48:26] Generating SAM header for ref/Striga
    [2015-01-05 09:48:28] Reading known junctions from GTF file
    [2015-01-05 09:48:30] Preparing reads
    left reads: min. length=30, max. length=101, 8733918 kept reads (706 discarded)
    right reads: min. length=30, max. length=101, 8732626 kept reads (1998 discarded)
    [2015-01-05 09:53:30] Building transcriptome data files str-1/tmp/StrigaC
    [2015-01-05 09:53:32] Building Bowtie index from StrigaC.fa
    [FAILED]
    Error: Couldn’t build bowtie index with err = 1

      1. 博主,您好,最近在做cufflinks与cummeRbund,但是总是在
        > cuff cuff
        Error in sqliteSendQuery(con, statement, bind.data) :
        error in statement: no such table: genes
        其他文件都没有错误的,后来发现在cuff_diff的输出文件中的cds所有的文件都是空的,不知道是在cufffdiff的哪一步出错了,cuffdiff的代码如下: cufflinks-2.2.1.Linux_x86_64/cuffdiff -p 16 -L A,B -u merged3_out/merged.gtf /share/BC_bam/3-1-A.bam,/share/BC_bam/3-2-A.bam,/share/BC_bam/3-3-A.bam,/share/BC_bam/3-4
        -A.bam /share/BC_bam/3-1-B.bam,/share/BC_bam/3-2-B.bam,/share/BC_bam/3-3-B.bam,/share/BC_bam/3-4-B.bam ,求助博主,可能是哪里出问题了呢?不知道是不是cds是空的才导致结果出错,还是本身就算有cds还是会出错呢?

        1. 我猜测你的cuffdiff没有运行结束,或者说它被终止了。cuff_diff一开始跑,就会生成一些文件占位,但最后会重写这些文件。你需要检查你的log文件,看cuffdiff在哪里出错了。从命令行来看,没有太大的总是。

  13. 博主老师好!请问cufflinks拼接基因所有的isoforms的问题!
    tophat有参基因组以genes.gtf作参Mapping,
    接下来De novo Transcript Discovery using Cufflinks
    命令是cufflinks -o OutputDirectory/ mappedReads.bam按理说cufflinks能发掘所有的isoform和新基因吧!
    然而结果是cufflinks不能拼出一个基因所有的可能剪切形式的isoforms,例如一个基因的3个isoform能拼接出来,而另一个可能的isoform就出不来!

    不知道-g 结合 -b同上述De novo Transcript Discovery所用的cufflinks -o OutputDirectory/ mappedReads.bam有何区别?
    我们的目标是拼接基因的所有的isoforms! 有朋友建议scripture,我想cufflinks若能即可!呵呵

    附:
    -g | –GTF-guide
    提供GFF文件,以此来指导转录子组装(RABT assembly)。此时,输出结果会包含ref
    erence transcripts和novel genes and isforms。
    -b | –frag-bias-correct
    提供一个fasta文件来指导Cufflinks运行新的bias detection and correct
    ion algorithm。这样能明显提高转录子丰度计算的精确性。

    1. cufflinks之后还有一步cuffmerge。如果你做的是没有参考transcriptome的,也许cufflinks并不是最好的选择.

  14. 谢谢博主老师的回复!cuffdiff是跑完了的,而且没有出错信息!我加了很多参数如下,用最宽松的方式跑cufflinks -p 4 -g genes.gtf -b genome.fa -F 0.3 -j 0.3 -a 0.01 -A 0.3 -min-frags-per-transfrag 结果cufflinks后文件大小倒是翻倍了,可是最终结果不理想啊!比如我们关注的一个基因的选择性剪切形式干脆消失了(只拼出来1个isoform)!还不如以前做的呢(3个isoform能拼接出来,而另一个可能的isoform就出不来!)!

  15. 博主您好
    我在重复Complementary Sequence-Mediated Exon Circularization 这篇文章的实验结果
    使用tophat-fusion的功能时无法跑出有效的junction.bed文件 都是只有几个Kb的结果 产生的accepted_hits.bam和unmapped.bam文件的大小总和也只有几个Mb。
    请问可能是什么原因造成的呢?

  16. tophat能不能用在没有内含子的原核生物?如果可以,内含子长度的选项如何确定

  17. 博主你好!大家好!
    我是新手,没有经验。问题比较小白。请您多帮帮忙。
    我在运行tophat的时候出现如下错误提示:
    [2015-01-26 14:58:24] Beginning TopHat run (v2.0.13)
    ———————————————–
    [2015-01-26 14:58:24] Checking for Bowtie
    Bowtie version: 2.2.4.0
    [2015-01-26 14:58:24] Checking for Bowtie index files (genome)..
    [2015-01-26 14:58:24] Checking for reference FASTA file
    Warning: Could not find FASTA file /shared/community/core/bowtie2/indexes/mm9.fa
    [2015-01-26 14:58:24] Reconstituting reference FASTA file from Bowtie index
    Executing: /shared/community/core/bowtie2/2.2.4/bin/bowtie2-inspect /mm9 > /home/student/ysu/tophat-out/tmp/mm9.fa
    [FAILED]
    Error: bowtie-inspect returned an error
    Command: bowtie-inspect –wrapper basic-0 /mm9
    请问这是路径出问题了么?谢谢。

  18. 楼主 你好 我想问一下 我在使用cuffmerge的时候 不使用参数-g -s ,我只是想单纯的把不同的Experiment下的cufflinks跑出来的transcripts.gtf文件合并到一起
    然后就出错了

    [Tue Jan 27 11:29:26 2015] Beginning transcriptome assembly merge
    ——————————————-

    [Tue Jan 27 11:29:26 2015] Preparing output location ./merged_asm/
    Warning: no reference GTF provided!
    [Tue Jan 27 11:29:43 2015] Converting GTF files to SAM
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:43] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:46] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:48] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:51] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:54] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:56] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:29:59] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:01] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:04] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:06] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:09] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:12] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:14] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:17] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:20] Loading reference annotation.
    gtf_to_sam: /lib64/libz.so.1: no version information available (required by gtf_to_sam)
    [11:30:22] Loading reference annotation.
    [Tue Jan 27 11:30:28 2015] Assembling transcripts
    cufflinks: /lib64/libz.so.1: no version information available (required by cufflinks)
    Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
    Command line:
    cufflinks -o ./merged_asm/ -F 0.05 -q –overhang-tolerance 200 –library-type=transfrags -A 0.0 –min-frags-per-transfrag 0 –no-5-extend -p 8 ./merged_asm/tmp/mergeSam_fileqDarXW
    [bam_header_read] EOF marker is absent.
    [bam_header_read] invalid BAM binary header (this is not a BAM file).
    File ./merged_asm/tmp/mergeSam_fileqDarXW doesn’t appear to be a valid BAM file, trying SAM…
    [11:30:28] Inspecting reads and determining fragment length distribution.
    Processed 36182 loci.
    > Map Properties:
    > Normalized Map Mass: 960360.00
    > Raw Map Mass: 960360.00
    > Fragment Length Distribution: Truncated Gaussian (default)
    > Default Mean: 200
    > Default Std Dev: 80
    [11:30:33] Assembling transcripts and estimating abundances.
    Processed 37927 loci.
    [Tue Jan 27 11:34:42 2015] Comparing against reference file None
    Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).
    [Tue Jan 27 11:34:49 2015] Comparing against reference file None
    Warning: Could not connect to update server to verify current version. Please check at the Cufflinks website (http://cufflinks.cbcb.umd.edu).

    1. 没看到错误的提示,全部都是警告。似乎停在了cuffcompare前后。你试着先跑一下cufflinks的示例文件,看看是否是程序问题,在确认程序没有问题之后,再检查一下磁盘空间。而后把中间文件检查一下。

  19. 博主,多谢您的回复。
    我没有找到mm9.fa的文件。请问这个应该喝index(mm9.*.bt2)在一起吗?还是应该在哪个路径下。我不知道这个文件时干什么用的。很难找到。
    多谢您的帮助。

  20. 博主,
    我才知道,我把-p设置太大了,所以被K下来。改小一点就可以跑起来。
    现在的情况时,tophat找不到mm9.fa,就自己创建了一个。然后我进去看了一下,全是NNNNN。我不知道是应该就这样呢,还是会有错。等跑完了,我再看看。可能还要您给指导。
    多谢您的帮忙。

    1. 你可以去http://support.illumina.com/sequencing/sequencing_software/igenome.html下载所需要的文件。
      fasta文件一开始和最后几千基本都是N,这个是正常的。

  21. 博主,
    我的设置果然有问题,您能帮我看看吗?谢谢!
    output文件夹下有三个文件。如下:
    drwxr-xr-x+ 2 ysu biostats 7 Jan 27 19:34 logs
    -rw-r–r–+ 1 ysu biostats 184 Jan 27 16:18 prep_reads.info
    drwxr-xr-x+ 2 ysu biostats 13 Jan 27 16:18 tmp

    tmp下的文件如下:
    -rw-r–r–+ 1 ysu biostats 1158893439 Jan 27 16:18 left_kept_reads.bam
    -rw-r–r–+ 1 ysu biostats 351800 Jan 27 16:18 left_kept_reads.bam.index
    -rw-r–r–+ 1 ysu biostats 241 Jan 27 16:19 left_kept_reads.mapped.bam
    -rw-r–r–+ 1 ysu biostats 0 Jan 27 16:18 left_kept_reads.mapped.bam.index
    -rw-r–r–+ 1 ysu biostats 241 Jan 27 16:19 left_kept_reads_unmapped.bam
    -rw-r–r–+ 1 ysu biostats 0 Jan 27 16:18 left_kept_reads_unmapped.bam.index
    -rw-r–r–+ 1 ysu biostats 2771195237 Jan 27 19:37 mm9.fa
    -rw-r–r–+ 1 ysu biostats 0 Jan 27 19:38 mm9_genome.bwt.samheader.sam
    -rw-r–r–+ 1 ysu biostats 1244987552 Jan 27 16:18 right_kept_reads.bam
    -rw-r–r–+ 1 ysu biostats 351962 Jan 27 16:18 right_kept_reads.bam.index
    -rw-r–r–+ 1 ysu biostats 0 Jan 27 19:37 temp.samheader.sam
    刚开始用,不知道怎么诊断错误。您能帮我看看可能是哪里出问题了吗?多谢

    1. 这个应该有log文件,会告诉你错在哪里。因为因素太多了,有些时候是程序问题,有些时候是硬件或者服务器限制问题,没有log文件,我也没办法。

  22. 楼主 我想问一下 cufflinks 参数设置 -G跑出的结果 所有的FPKM值都是0 是怎么回事

  23. 博主 我运行tophat遇到些问题
    Error: segment-based junction search failed with err =127
    /usr/local/bin/segment_juncs: error while loading shared libraries: libboost_thread.so.1.50.0: cannot open shared object file: No such file or directory

    此外,虽然有错误,我却得到了mapping 的结果……(比如accepted_hits.bam and the align_summart.txt,)
    请问这是怎么回事?

  24. 欧老师我现在运行tophat 用12个线程。还是得到这个错误 而且连accepted_hits.bam and the align_summart.txt都得不到了。
    Error: segment-based junction search failed with err =127
    /usr/local/bin/segment_juncs: error while loading shared libraries: libboost_thread.so.1.50.0: cannot open shared object file: No such file or directory
    请问如何解决?
    我以前一直可以用的。不知道是不是线程开多了还是哪里出了问题,,,,

    1. 你查一下你的libboost_thread.so是不是还正常。你使用ldd segment_juncs或者otool -L segment_juncs看一下链的是哪些文件。

  25. 欧老师关于illumina平台的数据,处理的参数都是类似的么?比如MiSeq,HiSeq,NextSeq,等等。是不是都是很类似的还是各有区别?谢谢

    1. MiSeq和HiSeq我是不分的。但是NextSeq的数据刚开始接触,还没有太多的经验可以分享。现在我们的想法还是先从以前的流程开始。

      1. 是的,我也是刚接触到nextseq的数据,所以就想都用默认的来试试。不知道妥否?(反正miseq我就是用默认的,此外illumina的adaptor应该是一样的,去adaptor和trim的流程差不多我想)如果有何发现希望和欧老师多多交流。

      2. 我用的是trim_galore做预处理的,细胞是人的,用的是length:20,quality:5,stringency:5,不知道妥否?

        1. 换我的话,我会把quality提高一些,至少要20,我一般设置为30。我不知道你的是不是有3’adapter,或者5’adapter。我拿到的一般都没有。所以设不设置stringency对我来说意义都不大。但是要设置的话,我想我会设置成3。当然这其实影响不大。而且我会跑两次,一次是不trim 3‘端的,一次是trim的。因为提交至GEO的时候提交不trim的会方便很多。

          1. 谢谢!请问您除了上述参数,trim的时候还有什么需要设置的?

  26. 你好,今天我用tophat2比对RNA-seq,出现了这个错误,请问这是怎么回事
    [2017-04-11 20:19:34] Beginning TopHat run (v2.1.1)
    ———————————————–
    [2017-04-11 20:19:34] Checking for Bowtie
    Bowtie version: 2.1.0.0
    [2017-04-11 20:19:34] Checking for Bowtie index files (genome)..
    [2017-04-11 20:19:34] Checking for reference FASTA file
    [2017-04-11 20:19:34] Generating SAM header for /home/share/data/pig_genome/pig102/update_20161212/updata_pig102
    [2017-04-11 20:19:36] Reading known junctions from GTF file
    [2017-04-11 20:19:42] Preparing reads
    [FAILED]
    Error running ‘prep_reads’
    Error: qual length (106) differs from seq length (150) for fastq record !

  27. 欧阳老师您好,为什么我的cuffdiff后的gene_expwen文件里所有显著的都没有loc号,然后我想着是在cufflink那一步没加-G,我加了之后,反而一个显著的基因都没有,请问一下这是为什么

    1. 你不加G,它就无法依照已知的信息来拼接transcripts,当然不会有基因名了。你加了之后,transcripts就不一样了,很可能差异表达的结果就不一样。当然,不同版本的cufflinks在pvalue的计算上不尽相同。

      1. 我在tophat 网上下载了果蝇的参考基因组,然后给了一个这个命令

        ln-s ./Drosophila_melanogaster/Ensembl/BDGP5.25/Sequence/BowtieIndex/genome

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

        因为是初学者。。。

          1. 写全了genome的路径 也是同样报错 后来用bowtie2重新构建了index 就好了。但不是很理解这个操作

nature进行回复 取消回复

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