RNA-seq差异表达分析工作流程

之前写过博文《从RNA-seq结果到差异表达》给出了一些关于RNA-seq分析的描述,这篇博文的目的是给出一个示例性质的工作流程。

需要使用到的工具:
TopHat
Cufflinks
Samtools

参考:http://vallandingham.me/RNA_seq_differential_expression.html

首先安装的tophat需要事先安装好bowtie。至于安装方面的问题,这里不至赘述。

整个pipeline非常明确:Sequences → TopHat → Manual Check → Cufflinks → Analysis

第一个问题,是否需要做duplicate removal,如果要做,什么时候做?在回答这两个问题之前,我们还是先来看看什么是duplicate。我们将deep sequence中完全相同的序列统称为duplicate。通常这种重复会有几个来源,一,测序模板中存在一模一样的片断;二,测序过程中PCR产生的重复;三,信号读取过程中读到了同一pcr产物。按照这里的讨论,对于 copy number detection, SV detection, ChIP-seq, and RNA-seq都应该做duplicate removal。去除的优点是可以大量的减少计算,降低假阳性。但是去除的话也有造成数据大量损失的风险,也就是说会降低真阳性结果。有文章对相同的library做了两次测序,一次是single end, 一次paired end。比较发现,SE的duplicate高达28%,而PE的duplicate只有8%。当把PE的结果当成SE结果来处理时,duplicate又升至28%。还有些私下的讨论认为,实际的duplicate应该只有1%左右。这里强调了去除duplicate对于数据完整性的影响。那么为什么人们在做CN/SV/ChIP-seq/RNA-seq的时候倾向于做duplicate removal呢?这主要的理论依据是在准备library的步骤中,所有模板小片段都是由超声波震断的,而相同的mRNA分子在同一地方被打断的可能性几乎为零。另一方面,当测序深度过深时,不可避免的,同一模板会被多次测序。这时候更应该去除duplicate,可以消除饱和。对于一些由酶切产生的片段,比如clip-seq, REDseq (Restriction Enzyme digestion sequence)等,就不需要做去除duplicate。在做去除duplicate之前,首先要在genome browser中观察一下mapped好的序列,看看其duplicate的存在的程度。肉眼观察这种事情,因为没有一定的标尺,所以非常不好总结。做这件事情的唯一好处就是,看得多了,就明白什么是好的测序结果。

那么duplicate removal什么时候做呢?现在的观点一般都认为是在map之后做。这样的好处是不依据序列一致就去除它,因为同一段序列可能map至不同的位点。在map之后,使用samtools rmdup或者Picard MarkDuplicates。这里需要注意的是,无论是samtools还是Picard,在duplicate removal时,所有的mapping结果对于每个read,应该只保留一个位置,而不是多个位置。2,对于PE的结果,LR的名字应该一致,否则程序可能无法识别。这些工具的出发点都是PE的,如果是SE的测序,可能需要指定参数。

 java -Xmx2g -jar /path/to/MarkDuplicates.jar INPUT=accepted_hits_sorted.bam OUTPUT=duplicated.removed.bam METRICS_FILE=picard_info.txt REMOVE_DUPLICATES=true ASSUME_SORTED=true VALIDATION_STRINGENCY=LENIENT

或者

 samtools rmdup   

第二个问题,bowtie的index文件哪里下载?在bowtie2的主站上提供了一个很有用的链接:iGenome, 这里集中了目前大部分的index文件以及相关的注释文件,可以很方便的下载。本教程就是使用从这里获取的打包文件。

tophat

tophat是针对mRNA-seq对bowtie的map结果进行了优化,它表现在两个方面,第一,从基因组中提取出mRNA junction的可能组合,对没有map结果的reads进行二次比对。第二,对于pe测序的结果进行依照mRNA测序的特点进行调整。它在参数设置方面,非常简单,只需要搞清楚几个重要参数即可。

Average Mate-Pair Inner Distance: -r/–mate-inner-dist

这个参数就是设置mate的两个测序reads之间的平均距离应该是多少。通常PE测序制库时通常片段大小的平均值为300bp左右,这个300bp包括了两端的adapter, barcodes, 以及序列本身,而-r参数需要设置的是两个测序结果之间的距离,所以它应该是总长-2*(adpter+barcode+reads)。如下图所示:
mateDist
我们需要注意的是,很多adapter可能比我们想象的要长,所以需要搞清楚具体的实验时的adapter长度。

Gene Model Annotations: -G/–GTF

通常这个参数我们并不设置。设置了之后可能会提高效率,但是也可能会产生倾向性。对于参照基因模型,通常的做法是在后面的步骤中再传入。

Threads:-p/–num-threads

线程数。现在多核处理器非常普遍,所以如果有四核的话,我们不妨设置-p 3,使用其中三个核,留下一核用于其它任务。但有一点要非常注意,如果你使用MPI的话,这个程序并不是MPI书写的程序,所以在cluster上运行时,需要申请独占模式,也就是使用-pe single参数而不是openmpi参数,并且使用-l mem_free=16G申请16G以上的内存以独占计算机,而不是多个内核。否则,多线程不安全。

Output:-o/–output-dir

输出目录。

–library-type

库类型,有三种,fr-unstranded Standard Illumina, fr-firststrand dUTP, NSR, NNSR, 以及fr-secondstrand Ligation, Standard SOLiD。默认为fr-unstranded。

下面是一个运行tophat2的例子:

export PATH=/share/bin/samtools:/share/bin/bowtie2/:/share/bin/tophat2/:/share/bin/python:$PATH
python /share/bin/tophat2/tophat --library-type fr-unstranded --mate-inner-dist 70 -p 8 \
        --output-dir tophat_output \
	/path/to/Genomes/UCSC/mm10/Sequence/Bowtie2Index/genome \
	R1_reads.fastq R2_reads.fastq

运行完tophat之后,我们需要检查一下map的质量,这时可以使用fastqc或者samtools来查看。但其实在run log中就已经记录了所有这些信息,可以从日志文件中查找我们需要的质量信息。

此进可以使用前文所述的脚本进行duplicates removal动作。

Cufflinks

对于Cufflinks,无论我们是否需要发现新的基因,最好还是走完三步。否则的话,结果可能会有一定倾向性。这三步是cufflinks, cuffmerge以及cuffdiff

对于每一个bam文件,我们都需要运行一遍cufflinks。

/share/bin/cufflinks/cufflinks -p 8 -o cufflinks_outputs accepted_hits.bam

然后我们将所有生成的transcripts.gtf文件的路径写入一个assemblies.txt文件中。我们知道,cufflinks的目的是通过重组mRNA来确定所有的transcripts,所以它会生成gtf文件来标记哪些是exon,这些exon组成了哪些transcripts。而下一步运行cuffmerge就是为了将所有发现的transcripts与已知的基因模型进行合并。
这里是一个典型的assemblies.txt的范例:

~/scratch/RNAseq/cufflinks_sample1/transcripts.gtf
~/scratch/RNAseq/cufflinks_sample2/transcripts.gtf

运行完cufflinks之后运行cuffmerge, 这一步需要把注释文件传入。

/share/bin/cufflinks/cuffmerge -g /path/to/Genomes/UCSC/mm10/Annotation/Genes/genes.gtf \
          -s /path/to/Genomes/UCSC/mm10/Sequence/WholeGenomeFasta/genome.fa \
          -o merged_asm -p 8 \
          /path/to/assemblies.txt

之后就是运行cuffdiff了。这一步将需要比较的两组的bam文件都传来。需要介绍几个参数,其中,–min-reps-for-js-test是指每组中至少多少个样品,默认为3,可以依据具体的实验调整。–labels之后跟着的是两个组标,用于标记传入cuffdiff计算的两个组的名字,注意不要把顺序标错。而之后的gtf文件就是第二步cuffmerge之后的merged文件。

/share/bin/cufflinks/cuffdiff -o cuffdiff_outputs -p 8 \
        --labels group1label,group2label --min-reps-for-js-test 3 \
	merged_asm/merged.gtf \
	gp1rep1/accepted_hits.bam,gp1rep2/accepted_hits.bam,gp1rep3/accepted_hits.bam \
	gp2rep1/accepted_hits.bam,gp2rep2/accepted_hits.bam,gp2rep3/accepted_hits.bam

结过了以上的步骤,我们就可以得到诸如gene_exp.diff, isoform之类的结果。具体的,可以在输出目录中去一一查看。

之后就可以使用R来分析过滤结果了。这里就不做介绍。

333 thoughts on “RNA-seq差异表达分析工作流程

  1. 接上条:
    3.博主有没有在玩biomaRT遇到
    Request to BioMart web service failed. Verify if you are still connected to the internet. Alternatively the BioMart web service is temporarily down.
    我的mart <- useMart("ensembl", "btaurus_gene_ensembl")
    可以执行
    但getBM(attributes=c("hgnc_symbol", "ensembl_gene_id"),filters ="ensembl_gene_id",values = EnID[,1], mart = mart))执行不了,老报上述错误

  2. 博主您好。我想知道Deseq2/edgeR可以做多类差分分析么?
    我看的是http://www.nature.com/nprot/journal/v8/n9/full/nprot.2013.099.html
    没介绍但我觉得肯定可以
    我现在三个类,每类三个样本。我想做的是C1 vs C2, C2 vs C3, C3 vs C1,但如果不能多类比较不是得跑三遍么?我想一遍跑完。

  3. 此外,关于阈值,cuffdiff是根据什么的?是FDR<0.05还是<0.1还是?谢谢!
    因为edgeR,DESeq要和Cuffdiff比较必须基于相同阈值

  4. 博主好,请问公司给的RNA-seq原始数据已经去除接头barcode,从这里开始tophat分析,请问参数–mate-inner-dist需要设置吗?不设置的话影响大不大?

    1. 最好还是设置。默认值为50,std为80.如果你运行完成之后,发现mapping比例太低,可以重跑一次。tophat1的时候,这个参数很重要。现在tophat2了,变得似乎没有以前重要了。

  5. 问题1:
    细菌的RNA-seq,我做完cuffdiff 发现转录本将几个基因合并了,例如
    原来的
    NZ_AXZL01000048.1:17578-19153
    NZ_AXZL01000048.1:19167-20946
    NZ_AXZL01000048.1:20961-21669
    生成
    NZ_AXZL01000048.1:17578-21669
    原来cufflinks做完时未合并,请问这是什么原因,(自己做的基因组拼接、预测、注释,是不是预测的不正确),能不能通过设置参数不合并。
    问题2:
    发现cufflinks cuffdiff的FPKM值标准化的默认方法不一致,请问需要设置成一致吗?不一样有没有影响?
    问题3:
    cuffdiff结果genes.read_group_tracking文件中包含的经样品间标准化的FPKM值(3个平行),是怎么计算得到gene_exp.diff中的value值的?
    备注:当时我做的cuffdiff命令如下,3个condition(F1,F2,F3),每个condition3个平行,打算做condtion间的两两比较。
    cuffdiff -o diff_out_3 -p 20 -L F1,F2,F3 -u merge_out_3/merged.gtf ./F1-1/accepted_hits.bam,./F1-2/accepted_hits.bam,./F1-3/accepted_hits.bam ./F2-1/accepted_hits.bam,./F2-2/accepted_hits.bam,./F2-3/accepted_hits.bam ./F3-1/accepted_hits.bam,./F3-2/accepted_hits.bam,./F3-3/accepted_hits.bam
    请举例说明。谢谢!

    问题4:
    不太理解 cuffdiff中 –min-reps-for-js-test 的意义,如果设置3的话,意思是是说3个平行中必须检测到3个不为0的FPKM值才能继续计算,若只检测到2个或1个,就不进行计算?像我的数据情况3个conditon,每个condition3个平行,做contion的两两比较,要设置为几合适?

    问题先问这么多,请多多指教!多谢@

    1. 您这是在考我呢?请举例说明……呵呵。

      你cuffmerge的时候有没有传入已知的基因组?如果没有,传入试试。看看它还会不会把已知的多个基因合并。

      问题2不是很明白。通常来讲,不同的比较对象,cuffdiff出来的FPKM值可能会不太一样,这个是可以接受的。

      问题3.value值的计算非常复杂,你可以去查看原始文献。近似的算法应该是3个平行的平均值。你可以计算确认一下。

      问题4 我一般设置为大于样品数的1/2。如果一组内有3个平行样品的话,我设置为2. 而且–min-reps-for-js-test 只在splicing.diff, cds.diff 和 promoters.diff 中起作用,不影响gene.diff的结果。

  6. 博主你好,一直不太cuffmerge是怎么工作的,assemblies.txt中必须要有一个以上的gtf文件吗,还有就是不同的gtf合并的时候是按照什么原则合并的呢,比如说不同的gtf有overlap,但是不完全一样,是都保留呢还是合并成一个exon?
    还有就是加了-G以后,怎么提取新发现的基因?其实我主要是想用cufflinks来找新发现的基因,而且是一个reads文件中的新基因,该怎么操作?谢谢

    1. 因为是merge,所以应该必须要一个以上的gtf文件。我还从来没有遇到过只有一个样品需要使用cufflinks比较的情况。合并原则需要你去阅读其原始文献。大约是如果有已知注释文件,那么先对已知的注释进行较正。如果是新发现的transcripts,就给个新名字。它的做法是先将所有输入的gtf文件转换成sam文件,然后运行cufflink来重新对transcripts进行拼接,然后得到新的转录本。这个过程去过滤一些假转录本。你可以试运行一下一个转录本的cuffmerge。但我不认为这有很大的生物学意义。

  7. 你好,我是没有生物背景的,机缘巧合,最近在做RNA-seq mapping。我看到网上比较常见的(paper也比较多的)pipeline是,tophat,cufflink,cuffmerge,cuffdiff,或者htseq,deseq2这样的。我也想按照这样的做,但是我老板不想用cufflinks。
    我想请教一下,如果只有tophat,就用htseq(gtf放ucsc genome网上可以下的gtf文件),和做完tophat,用cufflink(然后cuffmerge不同condition的gtf文件),再跑htseq,有什么区别么?

    1. 它的最大的区别就是你是否需要发现新的转录本。如果不需要,那么你没有必要去运行cufflinks。但是这个时候,你可能需要重新运行tophat的结果,把已知的gtf文件传给tophat,让它按照已知转录本优先的办法来重新map。htseq就是依据给定的转录组进行记数。

  8. 谢谢 博主 你的答案我正在慢慢消化,问题1中我在cuffmerge和cuffdiff中转入gtf 和细菌的genome信息,发现合并情况没有消失,只是比之前的结果中少合并了一些。我可不可以理解为这些合并的基因集是一个操纵子下的基因集。单独的基因的FPKM值就是约等于合并后的FPKM值。

    有一个疑惑的问题请教一下:
    问题5:
    FPKM值通过计算公式看出是map到基因的fragment数除以map到genome的所有fragment数与RNA的长度。
    那FPKM值与测序量有没有关系,还要不要再做样品间的标准化?

  9. 你好,关于pipeline的比较(除qPCR之外)纯从统计分析来比较有什么好的建议?
    我现在做了几个图。一个是几种pipeline预测的de genes的overlap,还有一个是前100(按P-value排列)的DE genes的spearman correlation。你觉得这可以? 我想从侧面比较方法的关联度,从而帮助生物学家在做qPCR之前决定哪个pipeline更合理

    1. 这个很难说。因为每种方法各有特点。我曾经使用过不同的办法分析同一组数据,相同的部分非常少。你可以将结果adjustPvalue 0.05 cutoff 之后按照abs(log2foldchange)逆序排列。之后,你可以提供给做bench work的人三组数据,top20各一组,overlap top20一组。然后决定哪个pipeline更合理。其实,使用多种pipeline的意义并不是很大。关键还是在于数据本身的质量。

  10. 博主 如果我有70个sample,每个sample大概30M的reads 做的就是差分比较(不同的time points,不同的conditions)。我想除了texudo还得用用edgeR和Deseq2之类的,多用几个pipeline比较。你客观的估计一下我这块大概需要多长时间(纯生物信息处理,不包括sequencing,从tophat2
    到下游de分析)。我老板要写项目本,需要问我大概时间。我估测大概要6个月不知道您觉得少了还是多了?求资深人士建议!

    1. 如果是使用cluster的话,一个月足够了。这只是计算而言。如果是单机的话,每组样品约需要1个星期。如果需要包括特殊分析,自己写代码,估计时间就不能确定了。我想计算力可以想办法解决,所以6个月时间有点多。你报一半吧。

  11. 博主您好,看你的博客学习很久了,学了不少东西!感谢!现在有一个问题请教您!我在用tophat2处理转录本时出现错误如下:
    [2014-06-03 20:42:44] Preparing reads
    left reads: min. length=36, max. length=90, 12693654 kept reads (198 discarded)
    right reads: min. length=36, max. length=90, 12693697 kept reads (155 discarded)
    [2014-06-03 20:53:28] Mapping left_kept_reads to genome ara_genome with Bowtie2
    [FAILED]
    Error running bowtie:
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    Error while flushing and closing output
    terminate called after throwing an instance of ‘int’
    (ERR): bowtie2-align died with signal 6 (ABRT) (core dumped)
    我的命令如下:
    ]$ tophat2 -p 8 -o /home/lxb/rna-seq/tophat_result/Car1_tophat_out /home/lxb/bowtie2_index/ara_genome /home/lxb/rna-seq/ara_rna_seq/clean_read/Car1_L1_1.clean.fq /home/lxb/rna-seq/ara_rna_seq/clean_read/Car1_L1_2.clean.fq

    clean reads是用Trimmomatic处理得到的,我怀疑是clean reads的问题,但不知道怎么解决,还望您指点一下!谢谢!

      1. /dev/mapper/vg_dell123456-lv_root
        50G 4.8G 42G 11% /
        tmpfs 16G 228K 16G 1% /dev/shm
        /dev/sda1 485M 38M 423M 9% /boot
        /dev/mapper/vg_dell123456-lv_home
        8.0T 112G 7.5T 2% /home

        硬盘应该没问题的~

    1. 你的问题可能是由于tophat和bowtie的版不兼容,你到官网上看看你所用的tophat应该使用哪个版本的bowtie。希望对你有帮助

  12. 大神你好。我现在是三组样本,每组三个。我已经完成了我认为生物信息学方面的工作 (生物学家想要我给他们提供差分genes) 。我已经用texudo(默认参数),counts-based(deseq2,limma,edgeR)两个pipeline得到了4个list。您说的做多pipeline也没太大意义。我就暂时做了2个pipeline4组结果。但是我老板想知道哪个pipeline最值得信任(在利用qPCR去验证这些基因)。我能做的就是通过venn图观察不同list之间的overlap,还有利用spearman correlation去计算4个list的关联度。我只能说通过这些相似性去给生物学家一点建议,缩小他们的验证范围。(最保守的当然是看几种方法都预测到的差分genes)。
    1。请问我这样做可以么?
    2。请问我在差分分析需要还做些什么?我现在就是得到了一些差分基因的list等待生物学家去 验证(当然那些软件得到的heatmap,MA plot什么的都做了)。我还需要比如改mapping的参数或者如改差分分析软件如cuffdiff,edgeR的参数再做一遍么?
    3。还是应该进行下游分析了?比如function annotation还有pathway analysis这些?
    Thank you!

    1. 以我的经验来看,texudo的验证率比较高。你这么做也没有什么错。这样至少你已经尽力了,做不出来也不能怪你。关于3,我还是非常推荐做下游分析的,因为这样做才有利于生物学家合理的安排故事,验证得再多不如一个完整的故事容易发文章。

  13. 谢谢。在Cuffdiff里 除了FDR<0.05外有没有其他删选方法?log2 fold change有什么比较经典的阈值么?谢谢!

  14. 博主我有两个问题
    1.假设一个gene有200个reads去map它,这其中有多少是有与内含子(intron)overlap的?我怎么知道???tophat可以知道么, 是否通过junction.bed?
    2. 如何知道reads里非mRNA的数量假设有核糖体RNA?
    Thank you!

    1. 1. 首先tophat会依据已知基因组或者de novo来map一次。而后依据这次的结果产生一个transcriptome。有了transcriptome,生成一个junction文件就不难了。我不明白你的问题,为什么是内含子,而不是外显子。
      2. 有专门的研究,你可以读一下SortMeRNA

  15. 谢谢楼主。你回答的很清楚。我想再问一个问题,在cuffdiff那一步,你的normalization的library的选择,是default的geometric,还是用quantile?

  16. 我的意思就是我想知道每个基因有多少mapped reads是跨内含子与外显子的。我需要统计。不知道什么软件可以? htseq-count可以么?
    谢谢你第二个网站,SortMeRNA, 我也想统计我的库里有多少rRNA,这个很有用。但我需要知道rRNA的注释文件么比如.gtf之类的?

  17. 谢谢!1.我摸索一下htseq-count, 要是搞不定再问你。
    2.关于rRNA的问题,你还没回答我。我需要其它什么东西么?比如rRNA的注释文件,还是就用你说的软件就行了?

  18. 有个幼稚的问题.
    假设我有段reads是TTAACG而reference genome是 …TTAGCG… 
    有一个突变,算不算map上?还是不算map上?tophat会如何处理这样的清况

    1. 我现在正在研究DaPars。它的程序上有一些bug,修复后可以跑起来,但是结果并不是很理想。之前试过一个叫HMMSplicer的,结果也一般化。你有了研究结果与我共享啊。

  19. Q1:你是如何判断HMMSplicer DaPars结果一般化的?通过IGV?请指教。
    Q2:关于rRNA的那个我还是看得不太懂请详解。
    谢谢你总是很及时的回复。

    1. Q1, 我们有几个阳对,已经实验证明了的,IGV验看也附合实验结果。但是这两个软件都算不出来。
      Q2,可能你的困难在于分类上。你去NCBI的分类部分去查你要看的物种,比如human,你可以从结果显示中看到每一层分类的名字(都是拉丁文,很难记的),对照着,你就可以找在arb上找到相应的物种的rRNA序列下载了。

  20. 请教您一个问题,我用htseq-count计数tophat 的结果,为什么计数不准确呢?和以前公司给处理的结果差了好多,有的基因计数为0,而实际有几百呢。还请您指教一下!

    1. 这个问题有很多条件没有限定啊,比如是否使用的相同的文件,相同的参数。htseq-count安装是否正确之类的,都没有信息。我首先会检查注释文件是否错误,其次使用igv具体查看几个例子,看看错误的发生在哪次结果中,因为这些比较好检查。

  21. 请教博主一个问题,
    有篇刊登在cell子刊的文章使用了多个来源的RNA-seq数据,并且建库方法不同(Poly+/rRNA depletion, SE/PE。并且使用RPKM 归一,没有提及其他的normalization的方法。
    而我们现在需要处理并对比多个建库来源的数据,请问文章的这个方法可靠吗?或者说有什么更好的方法?
    谢谢~

    1. 这个问题很复杂。有些时候,人们会使用不同来源的数据来相互佐证,这是可以的。有些时候,人们会将不同来源的数据进行相互比较,这个时候就要非常小心。对于测序数据,你说的不同的建库方法,以及不同的测序方法,然后使用RPKM来归一,表面上看,都是拿到了gene或者exon水平的标准化后的表达值,似乎可以相互比较。但是这个比较非常脆弱,除非是运气非常好,拿到很强的结果,否则得到的结果很可能是随机的,没有统计上的意义。但是如果只是相同实验内部比较一下,相互佐证,是完全没有问题的。如果你是需要比较不同来源的数据,这个结果只是做为初筛,是可以的。至少可以节约时间和精力,先拿到一个预结果。

  22. 我现在库里除了human gene外还有一条Luciferase gene。我还有的是Luciferase 的fasta序列。

    >Luciferase
    ATGGAAGACGCCAAAAACATAAAGAAAGGCCCGGCGCCATTCTATCCGCTGGAAGATG

    请问 如果我也想得到它的counts或者rpkm,怎么做?

  23. 你的意思用tophat不要map human的去map这个Luciferase ?可是我只有这个fasta没有gtf啊!没有bowtie2的index啊!能否说的详细点谢谢大神

    1. 自己做一个index就可以了。然后直接用bowtie来map就好了。或者你如果知道这段序列的位置,你可以使用R的Rsamtools先把bam文件读入,然后再使用GenomicRanges的findOverlaps或者countOverlaps函数来统计。

  24. 此外,接上问题,那我去拿这个库,用tophat map到human的genome上(我肯定是有human的fasta,bowtie2index和.gtf的),不会影响最终的表达结果吧(因为我要处理的库里混进一个Luciferase包括reads里可能有核糖体rna)?和一般处理步骤相同就行了么?不需要过多额外处理吧。

  25. 用tophat 对比两个fastq和生成参考基因转录组
    结果confflinks 生成的 transcripts.gtf 大小为 0b。
    请教一下大神是什么问题导致的。
    最后有这个日志。
    15774795 reads; of these:
    15774795 (100.00%) were unpaired; of these:
    15774767 (100.00%) aligned 0 times
    1 (0.00%) aligned exactly 1 time
    27 (0.00%) aligned >1 times
    0.00% overall alignment rate
    tophat产生的文件:
    -rw-r–r–. 1 root root 18810 7月 9 17:01 accepted_hits.bam
    -rw-r–r–. 1 root root 474 7月 9 17:01 align_summary.txt
    -rw-r–r–. 1 root root 233 7月 9 17:01 deletions.bed
    -rw-r–r–. 1 root root 452 7月 9 17:01 insertions.bed
    -rw-r–r–. 1 root root 52 7月 9 17:01 junctions.bed
    drwxr-xr-x. 2 root root 4096 7月 9 17:01 logs
    -rw-r–r–. 1 root root 286 7月 9 16:01 prep_reads.info
    -rw-r–r–. 1 root root 781166230 7月 9 17:03 unmapped.bam
    参考数据是genome

    1. 你map的结果太差了,只有一个uniquely mapped read。这很不正常。检查你的genome是不是正确,检查你的fastq文件是否正常,检查你的adaptor或者barcode是否去除了。

  26. 你好,您的回答我没有完全弄懂。如何利用fasta序列生成gtf文件?然后用bowtie2操作这个fasta文件生成bowtie2的index是这样么?

  27. 你好,我在运行cuffmerge的时候出现了错误
    命令:
    cuffmerge -g genes.gtf -s genome.fa -o merged_asm -p 8 ./assemblies.txt
    错误是很多下面的语句
    Error (GFaSeqGet): end coordinate (2644414) cannot be larger than sequence length 745751
    Error (GFaSeqGet): end coordinate (2659166) cannot be larger than sequence length 745751
    Error (GFaSeqGet): end coordinate (2733968) cannot be larger than sequence length 745751
    Error (GFaSeqGet): end coordinate (2732589) cannot be larger than sequence length 745751
    这个是什么错误呀?

  28. Error (GFaSeqGet): end coordinate (27841131) cannot be larger than sequence length 745751
    Error (GFaSeqGet): subsequence cannot be larger than 745751
    Error getting subseq for ENST00000378993 (28605516..29974840)!
    然后就退出了。

  29. 我的命令行:
    bowtie2 –local –very-sensitive-local -p 12 /home/LuciBowtie2Index/Luciferase -U /home/Trimed.fastq
    其中 Trimed.fastq是我的样本fastq
    而LuciBowtie2Index是我用bowtie2生成的 里面包括Luciferase的文件Luciferase.fa和bowtie2 index (e.g. Luciferase.1.bt2 )
    为啥出错?

    bowtie2-align died with signal 2 (INT)

    谢谢!

  30. 大神不知道你是否看到我的留言
    你让我直接bowtie2 去Luciferase的genome
    我的命令行:
    bowtie2 –local –very-sensitive-local -p 12 /home/LuciBowtie2Index/Luciferase -U /home/Trimed.fastq
    其中 Trimed.fastq是我的样本fastq
    而LuciBowtie2Index是我用bowtie2生成的文件夹 
    里面包括Luciferase的fasta文件Luciferase.fa和bowtie2index
    (e.g. Luciferase.1.bt2 等,你懂的)

    为啥出这个错?
    bowtie2-align died with signal 2 (INT)

    谢谢!

      1. 楼主你好,我在使用tophat+cufflinks的时候,第一次都加上了已有的注释信息,第二次则都没有加上,但cufflinks里面的FPKM结果则相差很大,比如说在基因组上同一位置拼接的转录本,未加参考注释信息的FPKM有两千多,而加了注释信息的则只有二十几,请问这是为什么呢?

        1. 按照您所描述的信息,主要原因是cufflinks对接接转录本的认定产生的差别。你需要比较相同转录本的位置是否严格一致,以及总的转录本的数量是否有显著差别。

  31. 你好 我试了 还是不行 一运行 屏幕上就是不停的滚屏 诸如
    “@HD VN:1.0 SO:unsorted
    @SQ SN:Luciferase LN:1653
    @PG ID:bowtie2 PN:bowtie2 VN:2.1.0
    SQE7Q:09324:09016 4 * 0 0 * * 0 0 TTCCACAAGCATTGCCTTCTTATTTTACTTCTTTTAGCTGTTTAACTTTGTAAGATGCAGAAACACTTCCCATTTTAGACAGCATCACAATGAGACAT 53658;;6;;;=5>97<49955423,19:46888-8:99899/828==4<::6;8879:9891878815803355)3446798886:::5:888:999 YT:Z:UU
    SQE7Q:09324:09031 4 * 0 0 * * 0 0 CAGTGCGTGCGGCTCTCGCGCTTTGAGAACGACCGCACCATCTCCTTCATCCCGCCCGACGGCGAGTTCGAGCTCATGTCCTACCGCCTCAACACGCATGTCAAGCCCT 58::878==;:5::::;:;;;=A3;;;;5:89;5;<?=5:;;;:7<7<2<6664::::;6888:;;::9;;;6<;;5;=7;;;5<;:::887786488517 YT:Z:UU

    我估计是bug 按ctrl+c最后一行就是 bowtie2-align died with signal 2 (INT)

    搞不定啊!!!!!

  32. 谢谢!
    如果我现在不是Luciferase,而是另一个基因(假设叫AN2 ),我有AN2的fasta序列(但这个基因目前是unannotated in the 目前的 human genome),也还有我实验样本reads的fastq(假设是人的,假设叫humansample1.fastq)。我想看AN2的表达水平是不是也是直接bowtie2生成AN2的index然后把整个fastq文件去map到AN2上?
    步骤:
    1.bowtie2-build /Downloads/AN2.fa AN2_bowtie2_index
    2.把这个index和fasta文件都放在一个AN2Index文件夹里
    3.bowtie2 –local –very-sensitive-local -p 12 /home/AN2Index/AN2_bowtie2_index -U /home/humansample1.fastq
    大神觉得我的做法和命令行对么?还是你另有想法?谢谢!

  33. 博主好,请问用cutdiff分析完后发现其中的很多看家基因(如50s ribosomal protein、gyrB)LOG2fold change 以及q值显示显著变化。这样的结果怎么解释?结果需不需要进一步处理!
    谢谢!

    1. 这要看具体情况。通常我遇到这种疑问时首先要做的就是:1.检查所有步骤是否正确;2.使用IGV查看这些看家基因mapping的效果,看看是否异常。然后再对其原因进行分析。

  34. 怎么查这些基因在基因组当中?没有很多基因,就一个。而且在gtf文件里是找不到的(没有注释信息).那我那样的做法是正确的么?大神

  35. 您好!我想请教一个问题,就是我在做完cufflinks之后可不可以把transcript.gtf文件中一部分基因比如说所有转录因子提取出来放到另一个文件中,然后用那个文件做cuffmerge?这样做会有什么影响么?

  36. 大神大神,还有一个问题,就是我在跑tophat的时候那个-r的参数应该是220,但是我写成了200,已经跑了比较久了,现在改的话实在很浪费时间,不知道这样有没有其他的补救的方法?或者它影响怎么样?十分感谢!

    1. 其实影响不太大。你可以这样,加上–no-coverage-search参数,快速地出一下结果。这样至少手里有东西,可以边往下做边等最开始跑的结果。

  37. 谢谢!我又有一个问题,就是那个fastq数据里面那个序列的长度length和tophat中-r参数有关系么?-r参数可以用那个length去计算么?跪谢跪谢!

    1. 我想这个问题,在文中的图例中应该回答得很清楚。如果你想准备地得到这个值,可以先使用tophat用默认的值map一遍,然后再用MISO提供的工具算出准备的值。之后再重新map一遍。

  38. 大神,我现在想做羊的mapping和差分分析。
    但是我没找到http://support.illumina.com/sequencing/sequencing_software/igenome.ilmn
    我就直接去ensembl上找 找到这个
    http://www.ensembl.org/Ovis_aries/Info/Index
    我想tophat需要annotation和dna sequence 对不对?
    可是当我点进 
    Download DNA sequence (FASTA)
    好多fasta啊(是按照染色体打包封装的) 
    请问是不是下载
    Ovis_aries.Oar_v3.1.dna.toplevel.fa.gz.?
    然后再下载ftp://ftp.ensembl.org/pub/release-76/gtf/ovis_aries/的gtf file?
    有了这两个就ok了吧?

  39. 大神大神,又有问题了,cuffdiff能同时做2个以上样本之间的差异分析么?比如我有五个样本,要比较他们之间基因表达的差异,我能不能同时放在cuffdiff里面去做呢?还是只能两两之间多次比较呢?跪谢大神,求解答。。

  40. 跪谢大神!另外还想请教一个问题,我用cuffmerge把transcripts.gtf文件合并后,跑cuffdiff之后产生的gene_exp.diff文件中基因的ID不是ENSG形式的而是XLOC编号,这个是为什么呢?要怎么把XLOC编号转换成ENSG编号呢?找了好久都没找到答案,求博主解答!

    1. xloc为前缀的是cuffdiff的分析ID。如果你的gene_exp.diff的gene栏仍然没有出现基因名,那么可能是你没有提供正确的transcripts.gtf或者你的cuffmerge的结果和已知的annotation完全对不上。你可以使用locus来重新做注释,得到gene symbol,这个可以试ChIPpeakAnno这个包。

  41. 请问有参的RNA-seq得到的差异基因绝大部分是找不到已知的locusID是怎么回事呢?不会全是newgene和ncRNA吧!测序样品有重复,tophat和cuffdiff都是默认参数。看了下18 Gb的sam文件,去掉flag>300的后剩余17 GB。不知道问题出在哪儿!谢谢

  42. 我又检查了一下tophat比对的结果,觉得这可能是tophat运行结果出了问题,我在用tophat的时候,第一次加上了已有的注释信息,第二次没有加,但两次read比对的结果却相差很大,相同的read在两次比对中比对到基因组的不同位置,其中不加注释文件的比对结果是正确的,加了注释文件的比对map的结果不对,90个碱基就对了几个,造成了同一区域内read数量差了很多,在添加了注释文件的比对中有很多是错误比对的read,不知道出现这种情况是为什么

  43. 求救大神!
    假如说我手头有三组对照A,B,C,每个对照下都有若干重复A1,A2…/B1,B2…/C1,C2…
    在cuffmerge这步我将三组全部写入assembly.txt 文件,运行后得到merged.gtf文件。
    现在问题是,cuffdiff我只想比较A和B之间,不想C参与进来,但是之前已经把C的结果merge到merged.gtf文件里了,还能用merged.gtf文件来做AB间的比较吗?C在merged.gtf中的参与会不会影响AB见的比较呢?

    1. 完全可以。是否会有影响主要取决于C的质量。如果C样品的质量高,会有正面的影响。如果质量低,可能会有一些负面的影响,但总体影响并不是很大。你注意如果单个的exon的转录本太多的话,可能是一个不好的信号。

  44. 请问一下,我在做cuffdiff 的时候,我用的GTF就是refseq的hg19的GTF,而不是用cuffmerge的结果,这样的影响是什么?
    我现有3个病人,每个病人肿瘤和周边组织各一套RNASeq数据,想找差异基因。
    谢谢!

    1. 如果你不关心新发现基因的话,影响并不大。需要注意的是,如果不关心新基因,在mapping的时候就使用注释库可能结果会更好。

  45. 我有6个样品,使用cuffdiff v1版,一起运行后,差异基因特别少,挑取其中一对样品结果看,差异基因比单独运行cuffdiff来比较这两个样品时得到的差异基因少近10倍??

  46. 博主,请问一下,我想获得转座元件的表达值,是否可以在Cuffmerge的时候自己构建一个转座元件的GTF参考基因组来使用,以求最后输出的只是转座元件的表达值?

    1. 完全可以。但是需要注意的是,如果你一开始就有这个文件,可以在map的时候就使用它,那样会得到更多的数据。

  47. 你好,我有两个问题
    1.我做RNA-Seq,map到human genome,DE gene list 结果出现了不少miRNA如MIR135B ,请问这正常么?我需要移除它们么?或者排除它们做一次mapping还是专门对miRNA做一次sequencing?
    2.我做cuffdiff发现很多DE gene的log-fold change 是inf但Q<0.05,我知道这表示什么意思,我想问下游分析这些基因该去除么?

    1. 1. 不需要。这是正常的。
      2. 两说。出现Inf,多半是因为其中一个样品该基因表达量太低。这时需要确认它们的表达值,因为p.value, FDR此时的作用并不是非常大。

发表评论

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