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

之前写过博文《从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来分析过滤结果了。这里就不做介绍。

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

  1. Reply yichenyiyi 10月 9,2014 3:14 上午

    你好,我用tophat做转录组数据与基因组序列比对时,发现生成文件accepted_hits.bam中并不是unique的,有的甚至还比对到了两条不同的染色体上,是需要设置什么参数吗?非常感谢。

  2. Reply yichenyiyi 10月 10,2014 9:51 下午

    我刚刚又重新看了tophat的详细参数,其中-g应该可以设置为1,这样得出的accept.bam文件中就是unique的结果。
    再问一个问题,我在做cuffdiff这一步时,如果加入-b参数(基因组fasta),结果文件中全是NOTEST,没有OK的基因,但是如果不加入-b参数,结果文件中就大多数都是OK的,这个怎么解释?谢谢大神。

    • Reply admin 10月 14,2014 8:06 上午

      你的理解和我的理解有一点点不一样。tophat中的-g参数只报告按得分降序排列的n个位置。默认值为20,是指如果有相同得分的话,随机从中取得20个位置。如果不足20个的话,有多少报告多少。但是要注意,即使把它设置也为了1,你得到的并不一定就是uniquely mapped reads,只是每个mapped reads的最高得分中的一个。这其间是有差别的。

      关于第二个问题,请仔细检查你传入的fasta文件是否与你mapping的reference完全一致。

  3. Reply Mars 10月 20,2014 2:50 上午

    我在按照流程做的时候,有几个问题:
    1. cufflinks的结果文件transcript.gtf中有些transcript_id:CUFF2.3或者gene_id:CUFF2等这些是不是代表新的基因或转录本呢?
    2. 我现在有癌症患者样本1,2,3和他们的癌旁组织样本,我该如何做呢?
    3. 不知道大神,知道如何鉴定转录组中的某个基因的AS(alternative splicing)?

  4. Reply pengli 11月 3,2014 3:58 上午

    我想问一下,我用htseq-count算counts数的时候,gff文件是应该用新组装的gtf文件呢,还是老的已公布的基因组注释文件?

  5. Reply jingming 11月 3,2014 5:02 下午

    博主您好,
    您知道怎么计算差异性表达的P alue值吗? 我这边只有Fold change.不知道如何计算P alue用什么工具?
    谢谢您!

  6. Reply jingming 11月 8,2014 4:52 上午

    欧博,你知道要有什么数据才可以计算P value 吗?我这有RPKM,以及RATIO。 麻烦欧博了。你和高博的那本红宝书没说道RNA_SEQ怎么处理这些啊

  7. Reply zzbin 11月 16,2014 9:37 上午

    博主,之后的R软件包cummerbund的使用是怎么样的呢?我在R环境安装好cummeRbund后,构建cummerbund数据库时,经常出现下面的错误提示:
    > cuff_data<-readCufflinks('diff_out')
    Error: could not find function "readCufflinks"

    • Reply admin 11月 17,2014 9:47 上午

      你是否调用了cummeRbund包?使用library(cummeRbund)之后才能找到相关的function.如果不是这种原因比较复杂的话,再留言讨论。

  8. Reply super 11月 18,2014 6:40 上午

    博主您好,
    1.我记得上次您和我说如果不要发现新的基因,可以直接cuffdiff省去cufflinks,我想再确认一下可以么?这样我也懒得用cufflinks和cuffmerge了。
    2.第二个问题,关于reference genome,有没有说哪种好?我发现很多人都是人用UCSC,动物用的Ensembl,这是不是有讲究啊?还是说习惯?
    谢谢!

  9. Reply zzbin 11月 18,2014 8:48 上午

    博主,你好!library(cummeRbund)最后会出现:
    Error : objects ‘metadata’, ‘isTRUEorFALSE’, ‘isSingleString’ are not exported by ‘namespace:IRanges’
    Error: package ‘Gviz’ could not be loaded

    继续构建cummerbund数据库时,仍然出现下面的错误提示:
    > cuff_data<-readCufflinks('diff_out')
    Error: could not find function "readCufflinks"

  10. Reply tianyi 11月 26,2014 10:50 下午

    求大神讲解一些为什么cuffdiff与edger找出的差异基因不同,原理是什么,求解释,谢谢!

    • Reply admin 11月 30,2014 4:06 下午

      这个我其实理解的并不好。也不是我这个水平可以一两句话说得清楚的。我试试讲几句吧,不正确的地方请指正。cuffdiff及edgeR的模型不同。正如你所说,cuffdiff使用的是FPKM值,而edgeR以及DESeq却并不提倡使用normalization之后的值,这是因为后两者都是使用的负二项式分布模型,这一模型是基于reads计数的。这是最主要的差别。其余的,比如噪音的控制,也就是background的去除等等细节,每一种方法都非常不一样。所以即使edgeR和DESeq的出发模型是类似的,最后的结果也会非常不同。

  11. Reply tianyi 11月 26,2014 10:51 下午

    是不是FPKM的影响?求大神讲解一些为什么cuffdiff与edger找出的差异基因不同,原理是什么,求解释,谢谢!

  12. Reply tianyi 11月 30,2014 9:06 下午

    谢谢你的回复,有些明白了,但老板总说这两个的差异到底在哪里?说从read看差异,我就没怎么看出来!
    还有个疑问就是比如使用cuffdiff中的–time-series,怎么对A组的A10_1,A10_2,A20_1,A20_2和B组的B10_1,B10_2,B20_1,B20_2进行.bam的时序比较,不是很明白,谢谢!

  13. Reply Asu 12月 11,2014 1:55 上午

    顶!楼主真无私!

  14. Reply shigongyao 12月 25,2014 1:25 上午

    你好!我有个问题想请教您。我有棉花两个亚种A2和D5的基因组信息,知道它们之间转录本的对应关系,现在有它们分别10天的纤维转录组数据。怎么比较它们两者间的差异。这个好像不同于一个参考基因组多个样本的实验,而是两个参考基因组对应两个样本。怎么操作啊?cuffmerge和cuffdiff怎么做?
    谢谢!

    • Reply admin 12月 26,2014 10:28 上午

      两个基因组之间的比较啊,我没做过。我的想法是,先把基因组比对一下,然后重建一个全新的基因组,之后将已经做好的mapping结果liftover到新的基因组上去(这一步可能最为困难,因为需要对bam/sam文件进行操作),同时做好新的annotation。之后就可以使用cuffdiff来进行差分分析了。

  15. Reply tiaotiao 12月 31,2014 12:24 上午

    想请教下,cuffdiff得到的结果中isoforms.diff与splice.diff的区别。感觉这两个都是和选择性剪切有关,但二个文件的具体区别不了解。希望指正

  16. Reply Tracy 1月 7,2015 8:39 下午

    博主好!
    用tophat2.0.11 将双端测序结果map到已知基因组,从align_summary.txt 文件可看到mapped 的read占输入序列的96.5%。如你所说,后续再做一步samtool rmdup -s accepted_hits.bam accepted_hits_rmdup.bam。问题1:-s参数是否正确,我用-S参数,去重复后的文件大小只有输入的10%。问题2:如果用参数-s,从文件大小和文件行数看不出输入输出结果的差异,请问如何评估做去重复前后结果差异。如果结果正常去除量百分比是多少。我做的是细菌RNA-seq。谢谢!

    • Reply admin 1月 8,2015 9:08 上午

      你是双端测序,依照说明,默认的就行,不需要加入任何-Ss参数,-s是为单端而设定的。所以是不正确的。我建议你试着使用一下MarkDuplicates。有些研究者认为,rmdup可能会有潜在的风险。

  17. Reply Chxu 1月 18,2015 9:25 下午

    博主的东西学习后收获很大,我做的Rna-Seq有两个生物学重复,如何检测两个生物学重复之间的相关系数呢?

    • Reply admin 1月 20,2015 11:10 上午

      这的确是一个好问题。现在的做法有两种,一个是基因水平,一个是外显子水平,先分别对它们进行不同水平的计数,然后对两个重复不同水平的计数做pearson correlation coefficient(PCC),你就会得到它们的相关系数了。计数可以使用htseq-count或者cufflinks。

  18. Reply super 1月 23,2015 11:46 上午

    博主您好,请问我2组数据
    组1 牛类样本1 样本2 样本3
    组2 马类样本1 样本2 样本3
    已经用tophat做好mapping了
    可以做差分DE分析么? 好像Cuffdiff不能做不同物种的差分比较吧?
    还是先把基因的ensembl ID统一再用 edgeR、Deseq ?
    谢谢!

  19. Reply super 1月 26,2015 11:28 上午

    什么叫各自定量?什么叫内参?能否详解,谢谢!

  20. Reply cuffdiff 2月 8,2015 9:09 下午

    博主您好,我在做了cufflinks之后做了cuffdiff,发现gene_exp.diff里的样本FPKM值和cufflinks结果里对应的单个样本的genes.fpkm_tracking里的FPKM值不一致,这个正常吗,是否和cuffmerge有关?谢谢!

  21. Reply super 3月 6,2015 6:37 上午

    老师,我运行tuxedo出现了点问题:
    当我运行tophat-cuffdiff时,出现警告
    Warning: couldn’t find fasta record for ‘GJ058256.1’!
    This contig will not be bias corrected.
    Warning: couldn’t find fasta record for ‘GJ058424.1’!
    ……
    (1)这什么意思?怎么造成的?是说我基因组不全么?请问影响后续差分分析么?(我用的是最新的Ensembl库)
    如果我用 Tophat2-Cufflinks/Cuffmerge-Cuffdiff2 ((默认参数设置)), Cuffdiff报错:

    Warning: couldn’t find fasta record for ‘GJ060129.1’!
    This contig will not be bias corrected.

    ……

    Error (GFaSeqGet): subsequence cannot be larger than 16338
    Error getting subseq for TCONS_00062149 (1..16448)!

    (2)请问怎么回事?是我哪里做错了?
    谢谢!

    • Reply admin 3月 6,2015 9:05 上午

      你使用的的tophat的genome和cuffdiff的是一样版本吗?这些报错我的理解是cuffdiff使用的注释在基因组中找不到或者位置信息对不上。

  22. Reply super 3月 6,2015 9:17 上午

    欧老师 我不太懂 我用的都是最新的版本 难道每一种版本还有版本之间的兼容问题?

  23. Reply super 3月 6,2015 9:22 上午

    接上条 我查了一下 我用mapping和cufflinks/cuffdiff的genome是同一版本的! 求欧老师详解

    • Reply admin 3月 11,2015 8:21 上午

      如果是同一版本,为什么会出现gff中的record不在genome当中呢?要不这样,你把gff文件中的chromosome和mapping的genome保持一致吧。

      • Reply super 5月 12,2015 9:30 上午

        欧老师您好,现在依然报错,我做的是牛细胞。用的是UMD3.1(最新版的)当参考基因组。tophat和cufflinks用的genome是同一个(版本也一样)。请问您说的保持一致是什么意思?怎么保持?

        这是UMD3.1版本的genes.gtf的第一列 我觉得应该是染色体信息 (我uniq单一显示了)
        1
        10
        11
        12
        13
        14
        15
        16
        17
        18
        19
        2
        20
        21
        22
        23
        24
        25
        26
        27
        28
        29
        3
        4
        5
        6
        7
        8
        9
        GJ058256.1
        GJ058424.1
        GJ058425.1
        GJ058430.1
        GJ058433.1
        GJ058437.1
        GJ058729.1
        GJ059463.1
        GJ059486.1
        GJ059509.1
        GJ059556.1
        GJ059670.1
        GJ060027.1
        GJ060032.1
        GJ060118.1
        GJ060120.1
        GJ060129.1
        MT
        X

      • Reply super 5月 12,2015 11:06 上午

        ps:欧老师我发现好像运行cuffmerge时就报错了
        Error (GFaSeqGet): end coordinate (83281041) cannot be larger than sequence length 81345643
        Error (GFaSeqGet): end coordinate (83282866) cannot be larger than sequence length 81345643
        ……

      • Reply super 5月 12,2015 11:45 上午

        请忽略上一条

        ps:欧老师我发现好像运行cuffmerge时就报错了
        Error (GFaSeqGet): end coordinate (83281041) cannot be larger than sequence length 81345643
        Error (GFaSeqGet): end coordinate (83282866) cannot be larger than sequence length 81345643

  24. Reply note 3月 11,2015 5:46 上午

    楼主,请问一般用cufflinks的参数-g重构转录本的时候运行时间是多久?我的数据是100M对的pair end测序结果,人的,都运行70个小时了,out文件夹下依旧没有任何东西。换成只有1M的数据量做测试,不久out文件夹下游文件生成,但是也运行了很久都没有结束。内存20G,-p 30。

    • Reply admin 3月 11,2015 8:24 上午

      你有几点需要改变一下,30个线程,没必要,我一般依据cpu的内核数来确定线程。如果只有4个内核,你只需要设置成4就好了。多了反而慢。你可以试着运行一下sample数据。或者查看log文件,看停在哪一步了。我一般4个小时之内就会出结果(数据为300M PE).

  25. Reply ws 3月 11,2015 8:47 下午

    博主,你好,请问你是否做过没有生物学重复的RNA-seq比较,比如,我有3个样本,但每个样本均只有一套RNA-seq数据,没有生物重复或技术重复,这时改怎么做?

    • Reply admin 3月 12,2015 10:21 上午

      没有重复的我做过挺多,但是发不了高分的杂志。你需要改是只是在cuffdiff时只输入要比较的样品就可以了。注意,在tophat那一步,我示例的是paired end的,那个R1和R2并不是指重复。

  26. Reply super 5月 25,2015 8:09 下午

    欧老师,如果我做reference genome的mapping,如何从map后的bam或者sam文件里求出每个特定基因的splice variance?或者说如何得到每个基因的consensus sequence ?谢谢!

    • Reply admin 5月 28,2015 1:40 下午

      你的这个问题很奇怪啊。你是需要做alternative splicing event detection呢还是想得到所有的transcripts?怎么会有一个consensus?这是个特有名词啊?

      • Reply super 5月 29,2015 5:31 上午

        可能说的有点雾,比如对于每个差异基因,如何从bam/sam文件里得到它们的transcripts的contigs?

      • Reply super 5月 29,2015 5:37 上午

        I would like a tool that generates a transcript sequence derived from a contig of the transcript sequences obtained in the RNA-seq mapped sequence data.
        欧老师,直接用英文表达或许歧义少点,因为翻译会可能把我老板的意思曲解了。
        谢谢!

        • Reply admin 6月 1,2015 10:21 上午

          你是想从RNA-seq里直接拿到transcripts的sequence吧?这个我的思路会是cufflinks拿到gtf文件,然后再提取

          • Reply super 6月 1,2015 11:54 上午

            我的理解也应该是
            gene1
            >transcript1
            AGC…..
            >transcript2
            ACG….
            gene2
            >transcript1
            ACC…..
            >transcript2
            AGG….
            这样的。

            请问您有什么方法和步骤么?

          • Reply super 6月 1,2015 11:56 上午

            补上一条,我cufflinks里的transcripts.gtf文件大概是这样


            10 Cufflinks transcript 555301 555870 1000 + . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; FPKM "8.5261709316"; frac "1.000000"; conf_lo "6.958111"; conf_hi "10.094231"; cov "39.076693";
            10 Cufflinks exon 555301 555870 1000 + . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; exon_number "1"; FPKM "8.5261709316"; frac "1.000000"; conf_lo "6.958111"; conf_hi "10.094231"; cov "39.076693";
            10 Cufflinks transcript 877473 880429 1000 + . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; FPKM "34.4084521901"; frac "1.000000"; conf_lo "22.954865"; conf_hi "45.862040"; cov "163.956835";
            10 Cufflinks exon 877473 877598 1000 + . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; exon_number "1"; FPKM "34.4084521901"; frac "1.000000"; conf_lo "22.954865"; conf_hi "45.862040"; cov "163.956835";
            10 Cufflinks exon 880363 880429 1000 + . gene_id "CUFF.2"; transcript_id "CUFF.2.1"; exon_number "2"; FPKM "34.4084521901"; frac "1.000000"; conf_lo "22.954865"; conf_hi "45.862040"; cov "163.956835";

          • Reply admin 6月 1,2015 4:02 下午

            你的cufflinks没有传原本的注释gtf文件啊?怎么都没有基因名?不过无所谓了,你可以后期再自己注释。你拿到这个transcripts.gtf文件之后,你可以使用R的GenomicFeatures来import它成为一个TxDb对象,然后使用Gviz的ranges(GeneRegionTrack(txdb))或者使用txdb的exonsBy方法来提取所有的exons的位置,然后你在使用getSeq方法从BSgenome中拿到序列,最后把它们拼起来就好了。

          • Reply super 6月 7,2015 1:23 下午

            欧老师 可不可以这样 我用htseq-count ,利用 – O 选项输出.sam文件。
            G219V:00087:05203 0 MT 6274 42 25M1D14M1D7M * 0 0 ACTCTCGCTCCCTGTATTAGCAGCCGCATCACAATGCTATAACAGA :<877:99345+336781876667+//5586641/..0+,+)+//2 AS:i:76 XN:i:0 XM:i:0 XO:i:2 XG:i:2 NM:i:2 MD:Z:25^G14^T7 YT:Z:UU XF:Z:ENSBTAG00000043561
            里面的最后一列不就是map到的gene然后第十列就是transcript么?

          • Reply admin 6月 8,2015 12:17 下午

            那你怎么拼成transcript呢?

          • Reply super 6月 9,2015 5:07 上午

            好像是的啊,我拿到reads后我不知道该怎么组装了。
            好吧欧老师我还是乖乖用您教的方法吧。我先试试,有问题我再来找您。
            有个问题,您在上一帖说我的transcripts有点怪,请问怪在哪里?
            chr1 Cufflinks transcript 139361 139843 1000 - . gene_id "CUFF.1"; transcript_id "CUFF.1.1"; FPKM "3.0069682237"; frac "1.000000"; conf_lo "1.524715"; conf_hi "4.489221"; cov "7.587669"
            是什么原因导致的?
            我用的都是默认设置啊
            cufflinks -p 8 -o Human_clout accepted_hits.bam
            cd Human_clout
            head head transcripts.gtf

          • Reply admin 6月 10,2015 2:01 下午

            cufflinks的时候不把你的参考gtf传进去吗?如果是无参的,那就不奇怪了。

          • Reply super 6月 10,2015 4:13 下午

            我有参考基因组啊 我是reference based的,用的牛的,人的都用,你是不是觉得gene_id 是”CUFF.1″很奇怪?我用的都是默认设置啊
            cufflinks -p 8 -o Human_clout accepted_hits.bam

          • Reply admin 6月 12,2015 7:53 上午

            你试一下在cufflinks的时候传一个-g/–GTF-guide 参数。

  27. Reply super 6月 9,2015 8:55 上午

    欧老师,我上面问了一个问题(我的transcripts.gtf有什么问题),此外我补问一个问题,有人在QQ群里说可以用gffread的gtf_to_fasta从GTF得到转录组层次的contigs。 还有人说可以用bedtools getfasta得到每个基因的转录本。
    请问这也可行么?

  28. Reply populus 7月 8,2015 12:31 下午

    博主,你好。请问如果我用DEseq2或edgeR进行RNA-seq差异基因的分析的话,如何获得counts值呢?是否可以用cuffdiff genes.read_group_tracking里面的raw_frags?谢谢

    • Reply admin 7月 8,2015 12:58 下午

      实际上对于DESeq2/edgeR都不建议使用cuffdiff中的raw counts,你可以通过查看它们是否都是整数来确认这一点。之前我们推荐的是htseq-count,现在流行的是RSEM。其R中也有一些包,但速度比较慢。

      • Reply populus 7月 8,2015 7:29 下午

        感谢回复。此外,还想问下我用JGI的基因组和gff文件去做的mapping,找出差异基因后,如何从JGI的gene ids转换为Entrez gene ids?谢谢

        • Reply admin 7月 9,2015 7:38 上午

          你可以使用JGI提供的转换工具试试,或者使用基因的坐标重新注释。因为我没有使用过JGI所以也没有太多的东西可以分享。

      • Reply super 7月 9,2015 5:27 上午

        那请问htseq-count和rsem各有什么优劣?为啥现在流行rsem?

        • Reply admin 7月 9,2015 7:34 上午

          其实我一直用htseq-count,效果也很不错。听说rsem更为准确。至于具体原理,我没有细纠。如果感兴趣,你可以自己看看。

  29. Reply super 7月 9,2015 11:49 上午

    欧老师我关于edgeR的分类问题:
    我做猪基因组rna-seq 每类3个样本 一共2类
    首先我通过htseq得到counts 然后分别用edgeR做 edgeR的基因count阈值我设为默认的1
    也就是这条命令: keep 1) >= 3
    但是分类效果很差 MDS图上看的六个样本都没完全分开 于是我把1改成了100
    MDS图上六个样本终于分开,三个在上面,三个在下面。
    三个问题
    (1).这种方法合理么?我指的是为了追求明显的差异基因结果,滤去了count table里的一些基因
    (2).如果合理,这是唯一的滤除方法么?有没有其它可以推荐?
    (3). MDS图怎么看,上三红点下三蓝点算分类合理了么?因为我之前一直是左三红点右三蓝点 (你知道的PCA图左右是PC1,上下是PC2),每次都能找到1000多个差异基因,但是如果是上三个下三个,差异基因的数目只有100个(都是默认设置)。抱歉您这个博客评论没法上传图片,所以我只好描述给你听了。

    • Reply super 7月 9,2015 11:50 上午

      刚才那个命令没打出来
      是 keep 1) >=3
      我给改成了 keep 100) >=3

      • Reply super 7月 9,2015 11:52 上午

        rowSums(cpm(counttable)>1) >= 3
        rowSums(cpm(counttable)>100) >= 3

        • Reply super 7月 10,2015 11:33 上午

          接上问: 欧老师,我这种方法滤去了低counts的基因,但是也会导致很多有价值的基因被滤去是吧? 我记得您说过如果edgeR的BCV图里离散度common那条线在0.2-0.4之内就很好,但如果大于0.5就有问题。表明样本差异量很大。现在我是3v3的样本,MDS图根本没分开(也没有明显的outlier样本),BCV图的离散度也大于0.5;去掉一对样本后,MDS图倒是分开了(2vs2),但是BCV图依然是common>0.5,这就是你说的不正常吧?那么有什么的好的解决方法?(我觉得我的Mapping和htseq应该没啥问题)难道要重新做library或者sequencing么?

          • Reply admin 7月 14,2015 12:13 下午

            其实BCV更重要一些。MDS分不分开无所谓。如果BCV太大,那么最终的结果中假阳性的比例会升高。最终还是要以实验验证的结果为准。

    • Reply admin 7月 14,2015 12:08 下午

      1). 除非你认为不同组之间的表达差异总体上会很大,否则没有必要,因为不同组之间大部分基因的表达都不会产生太大的变化。很多时候,软件包的设计者为了突出它的包的功能,会使用特定的样品来做示范,这个时候,你会看到非常明显的组间差异。
      2). 你过滤的时候可以把值调高一点点,比如所有基因上的counts的下1/4。太高了可能不好。
      3). MDS图就是为了强调组间差异的。如果组间大多数基因没有差异,或者这些差异是随机的,那么MDS图是看不出来的。至于左右,上下,其实无所谓,橫纵坐标和PCA图类似,只是两个最大差异的维度。

  30. Reply super 9月 10,2015 11:36 上午

    请问欧老师,做完tophat后如何求出每个基因的RPKM或者FPKM(如果我只有一个BAM file没法用cuffdiff计算)?如果直接用cufflinks,得到的genes.fpkm_tracking文件倒是有每个基因的FPKM,但是基因名都是CUFF.1 不像cuffdiff可以有基因名。
    可以用R做么?edgeR好像有rpkm()函数,但是需要知道基因的长度吧。。。
    所以我糊涂了。。。
    求欧老师解答

  31. Reply super 1月 20,2016 6:30 上午

    欧老师按照您的说法,-r/–mate-inner-dist 是tophat里比较重要的一个参数(这个是针对PE reads的),那么Cuffdiff里的–frag-len-mean 是不是也很重要(说明说是unpaired reads only,是不是也就是针对SE reads的)?

    • Reply admin 1月 20,2016 10:14 上午

      我想你对cuffdiff的–frag-len-mean的理解有些困惑。对它的理解应该是和tophat中的该参数是一样的。但是在最新的版本当中,cuffdiff会自己学习片段长度,该参数已经不推荐了。

  32. Reply super 1月 21,2016 10:09 上午

    欧老师关于GTF的问题我想请教。我发现在hg19里,发现很多基因在不同的染色体上都有片段
    比如IL9R在X和Y染色体上都有片段

    chrX unknown exon 155231073 155231173 . + . gene_id “IL9R”; gene_name “IL9R”; p_id “P20593”; transcript_id “NM_176786_1”; tss_id “TSS10431”;
    chrY unknown exon 59340194 59340278 . + . gene_id “IL9R”; gene_name “IL9R”; p_id “P21045”; transcript_id “NM_002186_1”; tss_id “TSS14565”;
    chrY unknown exon 59342487 59343488 . + . gene_id “IL9R”; gene_name “IL9R”; p_id “P4037”; transcript_id “NM_176786”; tss_id “TSS14565”;

    (1)但为何在Cuffdiff2的输出上只有Y染色体的?是不是Cuffdiff2默认就是找基因组最后一条染色体?

    $ grep “IL9R” gene_exp.diff
    IL9R IL9R IL9R chrY:59330251-59343488 C1 C2 NOTEST 0 0 0 0 1 1 no

    (2)此外是不是那些跨染色体的基因都是NOTEST的呢?

    • Reply super 1月 21,2016 10:28 上午

      我cuffdiff的 –multi-read-correct (use ‘rescue method’ for multi-reads )选项用的是默认的,FALSE。

    • Reply admin 1月 27,2016 10:41 上午

      我认为(1)这些应该属于非编码基因,它们有相同的序列,但是在基因组上有多个拷贝。Cuffdiff2应该不会做出选择。做出选择的应该是mapping的那一步。
      (2)并不是如此。依照定义NOTEST是指数据量不足以支撑test。作者并没有描述关于跨染色体基因的情况。

  33. Reply super 1月 27,2016 10:41 上午

    欧老师您用过Kallisto/Sleuth么?这个好像是观察DE transcripts。等我run出结果后,
    您也知道一条gene对应很多transcripts。像这种pipeline可以用来测DE gene吗?谢谢!
    target_id Q-value ext_gene
    ENST00000257570 0 OASL
    ENST00000339275 1.08823325997074e-138 OASL
    ENST00000543677 4.50865532164411e-20 OASL

  34. Reply super 4月 19,2016 9:36 上午

    欧老师您好,我的样本是人细胞感染了病毒样本。两组一组健康一组有病毒,各有三个fastq样本。我看两组的差分gene,我将样本map到人的基因组上,然后用tuxedo做了。现在我还想看是否有病毒的转录本transcripts的差异表达?怎么做?我现在已有的是病毒的转录本fasta文件(大概有二十个transcripts)。
    我的思路直接用Trinity的下游计算,比如RSEM计算出每个fastq样本的transcripts的counts,再用edgeR做,这个思路是对的么?谢谢

  35. Reply super 4月 27,2016 2:00 上午

    欧老师两个问题:1.如果我cufflinks 用-g 选项和-G 选项差别是什么?我看mannual觉得他讲得不是很清楚。2.此外以前记得您说过 ,如果我不想发现新的基因或者转录本只想用现在的基因组,cufflinks+cuffmerge这两步可以省略。您帖子里说的“可能会有倾向性”是什么回事?什么是“倾向性?”

  36. Reply super 4月 27,2016 5:51 下午

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

  37. Reply super 5月 3,2016 11:33 上午

    谢谢欧老师。我有一个问题请教
    我有六个文件,做差分分析。做tophat+cuffdiff,通过。做tophat+htseq+edgeR时,出现了问题。
    我竟然得到了一个错误
    “Error occured when processing SAM input (line 4240949 of file aligned.sam):
    (“Malformed SAM line: RNAME == ‘*’ although flag bit &0x0004 cleared”, ‘line 4240949 of file aligned.sam’)”
    [Exception type: ValueError, raised in _HTSeq.pyx:1294]
    比对好的SAM文件的第三列也就是RNAME竟然是 * 。

    我的命令:
    tophat2 -o outdir –library-type fr-firststrand –keep-fasta-order –GTF genes.gtf genome.fa Trimmed.fastq
    samtools sort -o aligned.bam -n -@ 8 accepted_hits.bam
    samtools view -o aligned.sam aligned.bam
    htseq-count -s yes aligned.sam genes.gtf>subset

    我感觉我的命令行没问题啊。。
    现在应该怎么办?
    是不是说只能用cuffdiff去得到结果了?
    我可以在SAM文件去掉这一行第三列是 * 的序列么?

  38. Reply JJ 5月 3,2016 10:40 下午

    博主你好,我想用cufflinks做转录组的拼接,用的tophat做的mapping,我的问题是这样的:
    1)做拼接时cufflinks -p -o -L -g 结果输出的gtf文件里 geneid有以-L 参数命名的 还有以gff注释文件里的ID为geneid的 为什么会这样呢?
    2)另外关于library type我也不是很清楚,我下载的是公共数据,我如何知道它在测序时是区分了正义还是反义呢?
    3)cufflinks里的参数library type 和tophat里面的–library type是一样的吗?因为我下载的数据质量值是sanger,所以我在tophat的时候我是用的缺省参数。

    • Reply admin 5月 17,2016 6:27 上午

      1)是正常的。因为有一些新发现的基因。
      2)这需要联系作者来明确。一般来说,他们有可能在描述上提到样品制作的方法,也许可以提供一些线索。
      3)是一样的。对于使用缺省参数,如果参数设置不正确的话,可能有少数reads的map不正确。但大多数应该是正确的。不过,这只是我个人的经验。最好还是使用正确的参数,以确保结果真实可靠。

  39. Reply lu 5月 4,2016 4:22 上午

    老师您好,我想知道cuffdiff有没有参数可以得到区分isoform的结果,cuffdiff得到的.diff文件将同一个基因的isoform合并,位置信息也只显示一个,我试着改gtf文件,有的成功的得到isoform的名字不同但是位置信息还是一样的,有的gtf试了好几种修改都不能运行程序,gtf格式看着是与原来一样的,是cuffdiff运行的时候还会识别gene_id,transcript_id的格式吗?非常感谢!!!

    • Reply admin 5月 17,2016 6:07 上午

      这也许算是它的一个bug。这只是因为它在最后注释的时候选择了基因水平的注释。所以你可以自己写一段代码来重新注释,或者使用cummeRbund来注释。

  40. Reply super 5月 5,2016 10:54 上午

    欧老师问一个问题,
    如果我有三个class ,想做差分分析。也就是想做 C1vsC2, C2vsC3,C1vsC3。
    我想问一下,做一次cuffdiff(C1-C2-C3)和我分别做三次cuffdiff(C1-C2),cuffdiff(C1-C3),cuffdiff(C2-C3),理论上来说是不是一样?

  41. Reply Desensitize 6月 5,2016 10:05 下午

    老师您好!最近我使用cufflinks遇到了一个奇怪的问题,在我选用-G参数的时候,得到的结果里所有的基因的fpkm都是0。但是我使用-g参数的时候,发现是有表达的(我随意选了一个gene分别在两个文件中grep),而且这种情况只发生在了部分的样本结果里,请问老师知道这种-G得到全部为0的问题该怎么解决吗?

    • Reply admin 6月 13,2016 9:21 上午

      你的cufflinks的版本号是多少?-G参数是从哪里来的?

      • Reply Desensitize 6月 13,2016 9:33 上午

        我的问题已经解决了,我最开始使用的是cufflinks 2.2.1,-G参数是以前的–GTF的缩写,但是我后来发现cufflinks官网已经宣布取消了对这个参数的支持。我觉得很奇怪的是不知道为什么仍然能够运行而且不报错。同时我某些样本的测序结果在最新版的cufflinks里使用这个参数就没有问题。。。

  42. Reply yang 7月 19,2016 7:08 上午

    您好!
    我想用cuffdiff的结果做火山图,可是这个结果的p值默认是小于等于10E-5的都等于10E-5,出来很难看,请问有办法获得准确的p值吗?

  43. Reply super 8月 3,2016 8:31 上午

    请问欧老师有接触过single-Cell RNA-Seq比如Cel-Seq么?我想知道和普通RNA-Seq的pipeline有没有什么差异。好像在测序中加了UMI,不知道这个该如何处理?需要去掉么?

    • Reply admin 8月 4,2016 2:05 下午

      我做的很少。只做过一组吧,还是和别人合作的,我做的下游分析。UMI是用来unique reads用的,在消除pcr bias之后就应该去掉吧,可以看这篇http://www.nature.com/nmeth/journal/v11/n2/full/nmeth.2772.html。
      有很多single cell seq的pipeline了,现在。你可以在网上搜搜,但是据我所知效果都不是特别好,其原因是很多人做single cell seq 去解决并非single cell的问题,因为他也不知道哪些是阳对,哪些是阴对,我们也无从比较不同分析结果的正确性。
      我在处理的时候,除了去barcode,UMI之后,基本上也试了些single cell DE analysis的工具,最后还是改成直接上DESeq或者edgeR或者voom了。不知道是因为数据不好,还是其它的问题。还有就是,normalization和quality control是最为头痛的问题。你可以在分析数据的时候仔细体会一下。

      • Reply super 8月 8,2016 5:32 上午

        谢谢欧老师,果然是大牛。之前我也谷歌过感觉做的人不太多。 (1)请问UMI如何去掉?是测序时机器会自动去掉,还是用软件去掉?我看UMI是好像是在reads的header上? (2)reads的fastq文件预处理后就直接做map么?用Tophat还是Bowtie?我看了一个CEL-Seq的pipeline竟然用的是bowtie? (3)能不能具体说说您觉得最为头疼的两大问题即归一化与质量控制。比如,Normalization的话不是edgeR或者DESeq会自动代我们去做么? 谢谢!

        • Reply admin 8月 8,2016 9:25 上午

          这个测序机器应该不会自动去掉。我的是自己写的perl去除的。你可以打开fastq文件看一下文件头,可以找到UMI是否已经被处理掉的线索。我建议你看个一千行,有没有UMI。网上有很多处理UMI的程序可以下载,你可看一下UMI-tools。有时候如果UMI太短了,可能会有少量的测序错误,这样可能会引起一些UMI的合并,有人会将map好的reads找回UMI再行判断一下,但是我没有做过这一步。

          对于mapping的话,用bowtie的话,可能不是最合适的,但也不是不可以,一般会再使用GSNAP再多找回些reads,或者直接map到transcriptome上。但是我没有看到任何理由不使用RNA的mapping公具的。

          关于Normalization和质量控制,主要问题是噪声太大。edgeR默认的是TMM,但是不一定会合适。现在估计有七八十种single cell analysis 的工具可以用,但是观点太多,无法统一。你在拿到mapping结果的第一步,应该是质量控制。有些测序会有内参,好整一些。如果没有的话,我的想法是直接使用mapped reads总数过滤一下。也有使用一些gene的expression pattern来过滤的。但是这需要实验支持。我当时就是最简单的处理,因为其实和实验设计的严密性有很大关系。测序的细胞数如果少于96的话,我估计几乎很难做出什么结果来,因为质量控制这一步,就有可能去掉不少。我能想起来的就这些了。

  44. Reply xueliang 8月 7,2016 3:39 下午

    楼主好,为什么我用bowtie map reads之后明明可以看到有reads map到某个基因位点且reads量还很多,为什么用cuffdiff计算该基因位点FPKM时数值为0?(在另外一个阶段,该基因位点的FPKM却不为0)

    • Reply admin 8月 11,2016 8:40 上午

      原因有很多。你先确认一下基因的大小,另外一个阶段的count数为多大。你可以通过cuffdiff的输出的track文件去查一下到底是哪一步出现了问题。

  45. Reply super 8月 8,2016 11:39 上午

    欧老师(1)你说的是这个么?这个感觉好像只能做SE而不是PE啊。。。https://github.com/brwnj/umitools (2)tophat/star可以么?您推荐什么呢?(3)“我的想法是直接使用mapped reads总数过滤一下。” 这主要是什么意思?

    • Reply admin 8月 11,2016 8:16 上午

      1) 我想这个问题以前有,现在应该是可以PE也可以处理了。你先试试吧。
      2)我觉着tophat是可以的。star我在其它场合试过,似乎mapping的效果不是很好,很多错配。我听其它人说mapping的效果挺好的。但是我的经历却并不愉快。
      3)我只是主观地设置了一条reads总数的期待值,少于这个数的样品就不要了。你可以按照自己的需求来进行质量控制。我这方面文献读的少。你多从文献中获得经验,我个人的精力总是有限的。

  46. Reply xiaohui 8月 26,2016 10:56 下午

    老师问一下,假如我想用edger来分析基因表达差异,我要进行哪些步骤,那个数据来源是从tophat结果还是从cufflinks结果来提取,还有怎么提取?

    • Reply admin 8月 29,2016 12:56 下午

      数据源一般是tophat map好之后,通过htseq-count来获得。也可以通过cufflinks的cuffnorm来得到,或者可以使用RSEM来取得。

  47. Reply 小洋 11月 15,2016 9:16 下午

    您好~我想问一下,我做到cuffdiff这一步的时候,为什么它会反馈说不能打开accepted_hits.bam这个错误信息呢?麻烦了

    • Reply admin 11月 16,2016 9:27 上午

      这个问题我不太明白。我猜你的accepted_hits.bam是tophat输出的结果。你是否应该确认一下你的cuffdiff的输入是否指到了正确的文件路径?

  48. Reply cheng 2月 13,2017 6:28 上午

    你好,求教一下 cuffdiff只能 两两间比较,有没有方法 可以看出,基因在 多个(比如6个)组织间 是否显著差异

  49. Pingback: RNA-seq差异表达分析工作流程-2 ← 糗世界

  50. Reply YAO 5月 31,2017 10:11 上午

    博主你好,我在做tophat–cufflinks–cuffmerge–cuffdiff的时候遇到这样一个问题。cuffdiff输出的结果里边gene id都是XLOC的编号;想问一下怎样更改为参考基因注释里边的gene id呢?
    另外还有一个疑问,就是cuffdiff输出的expression dif文件里边会出现一个XLOC值带领的一行对应有两个甚至两个以上的gene name,这样的话,对于具体是哪一个基因产生了变化是怎么判断?或者这种包括了两个或者两个以上基因的结果就舍去不要么?

    谢谢!

    • Reply admin 5月 31,2017 10:24 上午

      如果你在cuffmerge的时候提供了基因组的gtf文件的话,你就应该会得到相应的注释ID。看你的第二个问题,似乎你自己已经解决了第一个问题。

      如果一个位置对应多的基因名,有时候是因为它们重叠在一起了,无法区分,只能把两个相近的基因放在一起了。你可以参考一下:http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/index.html当中的class code再做判断。

Leave a Reply

  

  

  

%d 博主赞过: