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. 博主你好,现在关于RNA-seq有很多软件。
    又是Cuffdiff的(用的好像是FPKM作为表达值),又是DESeq的,又是edgeR的(用的好像是map的基因数做表达值,我没说错吧?),五花八门。请问哪一种是比较好的?
    我现在的pipeline 就是你说的这个
    tophat+bowtie — cufflinks –cuffmerge — cuffdiff 然后用cummeRbund去看结果。

    1. tophat+bowtie — cufflinks –cuffmerge — cuffdiff 然后用cummeRbund去看结果应该是主流。但是它有缺点,比如它不能做复杂的比较。比如你有不平衡的样品设计,有多组之间的,有时间序列的。在这些时候,还是使用edgeR或者DESeq比较合适。不过,因为它们的输入是reads matrix,所以之前需要做htseq-count来对mapping的结果进行计数。当然现在R出了一个summarizeOverlaps,但就我个人使用的感觉来说,不如htseq-count的效率高,不过它在不断改进中。

  2. 你好,什么是您所说的复杂的比较,能否说的具体一点。多组之间就不好了么?比如三个condition的话,cuffdiff 就不适合了?不平衡又是什么意思?辛苦了!

    1. 比如说你有两组样品,每组样品分别是在时间点0,1,2,3取了样,每次取样取了两个重复,这时候,你想比较的可能性有很多组合,在这种情况下,cuffdiff就力不从心了,它只能两两比较,比如可以在时间点0的两组样品间比较,或者时间点1的A组样品和时间0的A组样品比较,但是不能同时考虑其它时间点。如果你了解limma的话,使用edgeR就可以进行更复杂的比较了。当时除了时间系列,还有比如因子,权重之类的。
      你讲的三个条件,简单的讲,这三个条件中,两两可以互相比较,但是你无法判定A和B,C相同比较的情况。
      样品数不平衡是指一组有两个样品,但是另一组有三个样品,这种比较就属于不平衡的比较。

  3. 多谢!那么tuxedo系列相比起R语言base的edgeR/DESeq等优势在哪?为什么成为主流?是不是差异表达基因分析的更准确还是更利于下游分析?不会因为它是全程shell command line是优势吧?

  4. 补充上一条:是不是运行时间上也有优势?

    此外关于下游分析,请问博主都有哪些?我现在知道的是找出重要的差异表达的基因进行GO注释,通路分析等。还有没有?谢谢!

    1. 下游分析依赖实验需要。对于RNA-SEQ通常也就是GO,KEGG分析。其它的我一时也想不起来。不过我做过很多结合其它seq结果一并分析的。

  5. 还有,DESeq和edgeR各自有何优缺点不知道博主是否有所研究?
    合作者告诉我他可能有4-5种animal species,每种物种可能会重复run2-3次(可能在不同的condition,比如treated 和 untreated cells,我这是病毒),所以至少要run Ion proton 20次。我这种情况应该是肯定要用edgeR/DESeq了吧?

  6. 此外对于您说的“比如说你有两组样品,每组样品分别是在时间点0,1,2,3取了样,每次取样取了两个重复,这时候,你想比较的可能性有很多组合,在这种情况下,cuffdiff就力不从心了”我不是太明白。
    比如http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.016.html#an3
    (这篇最经典的介绍tophat-cufflinks-cuffdiff的论文),它的样本是
    C1_R1_1.fq ,C1_R2_1.fq, C1_R3_1.fq (我理解为condition 1 三个时间点或者condition1 的一个时间点取了三个重复不知道对不对)
    C2_R1_1.fq,C2_R2_1.fq,C3_R3_1.fq (我理解为condition 2 三个时间点condition2 的一个时间点取了三个重复不知道对不对)
    您的意思是不是如果是两个时间点两个重复,就没法用这个pipeline了。是这样么?抱歉回复了您很多条,让您看的费劲了阿

    1. 看来我还是没有讲清楚。这样讲吧,无论一个实验条件下有多少个重复,对于cufflinks而言都可以接受。如果有两个时间点两个或者多个重复,只要你不需要考虑复杂的组合,只需要进行时间点1的两组进行比较,或者时间点2的两组进行比较,都可以使用cufflinks。但是如果你要考虑时间点1和时间点2会不一样,而你又需要比较两组(包括两个时间点)时,cufflinks就不如deseq或者edgeR。

      对于deseq或者edgeR之间的比较,这是一个比较复杂的问题。因为两者基本上都使用了类似的模型。我也没有看到很好的比较文章。

  7. 谢谢,回答很有用。我按照流程试着玩了一下egdeR感觉也还蛮不错 ,纯R分析应该比较牛,有很多数学方面统计,但比起来不足之处好像只能输出重要差异表达的基因,是这样么?还是其它功能我还没玩到?因为cuffdiff好像除了给出最重要的差异表达基因,还有transcirpts,CDS,TSS等,不知道这些R的工具包也能做到么?
    但为什么不是基因symbol例如MKNK2,TCF3这样而是FBgn0000003这种类型的ID(HTseq的输出),请问这是什么ID?

    1. egdeR和deseq暂时还没有相关TSS的分析。对于transcripts, CDS,exon等你可以通过设定不同的gff文件来实现。但这并不是一个很好的选择。
      至于基因symbol表示成trackID的,很可能是你cuffmerge或者cuffdiff时输入的gff文件不一致或者说你发现了新的transcripts。

  8. 请教一个问题:我们的RNA-seq数据中有一个基因,在cufflink分析后,两个样品中FPKM值分别是45/0,但是在使用了cuffdiff以后则显示为78/134(尽管是“两个样品”,其实它们起始的RNA库是一样的,只是有处理和不处理的区别),这样的结果是否存在异常?

  9. 您好,我没有用cuffmerge,我直接cufflinks完就用htseq-count得到的counts矩阵(按照指南上一步步操作的),然后输出就是FBgn0000003,我也不知道这是什么基因ID看起来满怪的。此外我在“ NGS分析入门” 板块也给您留言了,还是关于预处理的问题。“那么FASTQC查看质量不好的(比如小于20)怎么剔除?RNA-SEQ需要不需要剔除很短的序列?或者N多的?或者重复的reads?有什么软件么谢谢。”。

    1. FBgn0000003看起来好象是flybase的id。你从这个ID很容易找到它对应的symbol。之所以你得到的是这样的ID,多半是因为你在cufflinks过程中提供的gtf文件中相应的栏是使用的是flybaseID。从flybaseID至symbols,你可以使用bioconductor的相关注释工具,或者使用在线的转换工具,比如由flybase提供的转换工具:http://flybase.org/static_pages/downloads/IDConv.html.
      至于如何过滤fastq文件,一般我是通过bioperl来编写自己的小程序。其它的工具,可能fastx是比较常见的工具:http://hannonlab.cshl.edu/fastx_toolkit/

  10. 谢谢。关于mapping我的问题是,我现在用了一个果蝇的raw data在玩pipeline,我发现tophat2的mapping rate只有不到60%。这是不是太低了。在RNA-seq中这么低的rate会不会影响最后的结果,因为表达值不就是FPKM(cuffdiff)和mapping到gene的reads数量(edgeR)么?

    1. 对于果蝇RNA-seq来说,如果mapping后的reads有超过10M就可以了,无论mapping rate是多少。mapping rate低,只能说明实验过程中有污染或者空载体以及小片段过多。

  11. 此外您说fastx_toolkit我也听说过,就是不知道该剔除什么?做什么?比如剔除FASTQC里得分不好的?还有什么?请指教。

  12. 如果不是果蝇呢?我只是拿这个dataset 玩玩,真要实际数据应该是家畜或者人的。所以特别想知道拿到别人的raw data ,在mapping 前的 pro-processing要做些什么?能否说详细点呢谢谢。
    另外,我看有人推荐说可以用tophat+bowtie(第一次tophat unmapped的reads)提高rates,请问这样对后续处理有用么?谢谢!

    1. 如果是人的话,最好能达到300M reads。宽松一点的话,150M reads也不错。计算方法很简单,使用genome size/reads length约等于10就好了。比如果蝇基因组大小为168M,如果你的reads length为100的话,那么10倍深度就是大约17M.但实际的transcriptome不可能有那么大,通过许多别人的实验数据表明,10M左右对于果蝇就够用了。其它家畜你按照这个大约估计即可。关于是否需要使用bowtie做二次mapping的问题,个人认为完全没必要。首先tophat本身就是调用的bowtie或者bowtie2来进行mapping。其次,tophat抑制输出的部分多为intron或者基因间的部分。提高rates多是本末倒置,把本来应该去除的噪音又找回来了。

  13. 1。 谢谢,您的建议mapping这个步骤就是还是只用tophat就可以了是么?
    2。 那么之前的预处理有哪些呢?您说的“如何过滤fastq文件,一般我是通过bioperl来编写自己的小程序。其它的工具,可能fastx是比较常见的工具”这是处理什么的?重复的reads?过短的reads?还是Quality Control得分过低的reads?我主要想知道如何预处理?还是您说的其实不需要因为tophat可以做,rawdata拿来直接mapping就可以了?谢谢!

    1. 1, 是的,只用tophat就可以了。2,主要是看你的reads的质量.我们知道bowtie在mapping的时候是不考虑quality的。如果你的reads的质量大部分还是不错的话,就没必要过滤。如果质量真的是不行,可以考虑过滤。但是具体的过滤方法,并没有一个非常统一的流程。你可以试着去galaxy去找找看,那里有很多关于trim的。对于这方面,我没有太多使用过的经验。

  14. 谢谢。
    还有个问题,为什么不直接用bowtie/bowtie2去做mapping 而要用tophat呢?
    另外如果有duplicate是不是就不管它还是说等mapping之后去除呢?
    您这篇博客的意思是不是Single Ended的就要去除,Pair Ended的就不需要了?

    1. 你需要去读一下tophat的文章就明白了。简单的讲,bowtie是map到genome,而tophat是map到transcriptome.tophat就是通过多次调用bowtie来完成map至exons和junctions这个过程。关于uniquely mapped reads问题,你可以在我的博客里搜一下,应该有讲过。

  15. I found a problem that gene ID output from Cuffdiff and egdeR/DESeq are different so that I cannot find the common DE genes.
    你好我又有问题
    两个pipeline
    Tophat-Cufflinks-Cuffmerge-Cuffdiff
    Tophat-edgeR/DESeq
    我用的是最经典的那个“Tophat-Cufflinks-Cuffmerge-Cuffdiff”文章的果蝇数据集。
    但是
    The edgeR output 的gene ID 是(e.g. FBgn0000370,FBgn0000500,…)
    而 Cuffdiff output gene ID 是 (XLOC_000028,XLOC_000038,…) 以及还有 gene symbol 比如,KH1,RpLP1,…).
    如何能让它们统一?因为我想比较overlap的DE genes。

    PS:我的file
    genes.gtf的第一行
    2L protein_coding exon 7529 8116 . + . exon_number “1”; gene_id “FBgn0031208”; gene_name “CG11023”; p_id “P9062”; transcript_id “FBtr0300689”; transcript_name “CG11023-RB”; tss_id “TSS8382”;
    merged.gtf的第一行
    2L Cufflinks exon 7529 8116 . + . gene_id “XLOC_000001”; transcript_id “TCONS_00000002”; exon_number “1”; gene_name “CG11023”; oId “FBtr0300690”; nearest_ref “FBtr0300690”; class_code “=”; tss_id “TSS1”; p_id “P1”;
    可不可以用merged.gtf 去生成edgeR/DESeq的counts table呢?

    1. 可以使用merged.gtf去生成edgeR/DESeq的counts table。
      如果想统一gene symbol,你需要将注释统一即可。有很多工具,我推荐由我维护的ChIPpeakAnno软件包。

  16. 谢谢。小问题又来了。
    为什么在cufflinks 的示范教学论文上,
    tophat 是这样使用的
    tophat -p 8 -G genes.gtf -o C1_R1_thout genome C1_R1_1.fq C1_R1_2.fq
    tophat -p 8 -G genes.gtf -o C1_R2_thout genome C1_R2_1.fq C1_R2_2.fq
    tophat -p 8 -G genes.gtf -o C2_R1_thout genome C2_R1_1.fq C1_R1_2.fq
    tophat -p 8 -G genes.gtf -o C2_R2_thout genome C2_R2_1.fq C1_R2_2.fq
    tophat -p 8 -G genes.gtf -o C2_R3_thout genome C2_R3_1.fq C1_R3_2.fq
    tophat -p 8 -G genes.gtf -o C1_R3_thout genome C2_R3_1.fq C1_R3_2.fq
    也就是每个condition每次重复每次去mapping,然后再组合,然后再DE分析
    而在edgeR/DESeq的示范论文上,
    tophat是这样使用的
    tophat2 -G Drosophila_melanogaster.BDGP5.70.gtf -p 5 -o Untreated Dme1_BDGP5_70 \ SRR031714_1.fastq,SRR031715_1.fastq SRR031714_2.fastq,SRR031715_2.fastq
    tophat2 -G Drosophila_melanogaster.BDGP5.70.gtf -p 5 -o treated Dme1_BDGP5_70 \ SRR031724_1.fastq,SRR031725_1.fastq SRR031724_2.fastq,SRR031725_2.fastq

    后者是把每个样本所有的重复放在一起mapping然后去生成counts table。请问有啥区别?而前者为什么那么麻烦呢?
    谢谢!

    1. 我一般以tophat上的为准。edgeR/DESeq上的放在一起的这个,我不推荐。还有,你的这个示例是不是有点小问题啊,文件名似乎对应不上。

  17. 谢谢。如果我的library不止mRNA,还有miRNA,rRNA这些,那么需要有什么预处理么?还是说和一般的RNA-Seq预处理相同(我目前还没遇到这种情况只是感兴趣问问)

  18. 另外RNA-seq的步骤是否相同?(如果我的library不止mRNA,还有miRNA,rRNA这些)谢谢!

  19. 您说步骤相同,但是您上次说RNA-Seq直接tophat 也就是说miRNA,rRNA也是直接tophat就行么?谢谢!
    问题又来了,
    我看到很多文章去比较edgeR/DESeq/Cuffdiff/baySeq的性能(我相信您肯定没少看 )。
    当然它们多用的是simulation data, 因为我想如果是real data,怎么评判你的假阳性率就是个问题了,除非有实验验证证明它们的输出确实是差分表达的基因。
    请问那些 simulation data是怎么做的阿我都没搞清楚。或许是我没读太清楚。

    1. 对于simulation data,很多都有倾向性。你如果做了解的话,其实可以试着自己做一些,你需要首先了解一些基本的概率分布,然后就可以自己生成一些数据,再加上一些自己刻意准备的特殊情况。对于评判假阳性,你可以计算一下FDR。依据我做RNA-seq的经验,cuffdiff的结果还是比较准确的。如果有复杂设计,我才会考虑其它的方法。

  20. 博主你好,在使用Tophat 时候,我已经把 -p选项改为1,为什么还是有bowtie2-align fix_map_orderin bam2fastx,三个线程出现,最近服务器使用 有限制,我想每个fastq文件只用一个线程,可以多跑几个fastq,谢谢

  21. 您好,问一个幼稚的问题。请问关于RNA 库,里的您引文中的“fr-unstranded,fr-firststrand和fr-secondstrand Ligation, Standard SOLiD。默认为fr-unstranded。” 请问这是不是把reads map到transcirpts 上的方向阿。
    这个和single-ended 和pair-ended序列是一个概念么?我的意思是不是SE 就一定是fr-unstranded的。谢谢!

  22. 同上一个问题再加问一个naive的问题,这个RNA库的fr-un,fr-first等我怎么知道呢?是要问测序的人还是说可以从FASTQC里看?

  23. 晕,今天测序的人告诉我他是用的stranded,我却用的TOPHAT默认去mapping的,看来得重新作mapping了…..

  24. 好像不仅是mapping, htseq-count 和 cufflinks 都要设置是不是stranded的是吧?晕!默认都是unstraned。哎!悲催!

  25. 博主您好,我这边另一个测序员(我有几个数据)说他也不清楚是自己的library type是stranded 还是unstranded的。只说保留了极性并且序列是有方向的。我不太懂这个意思是不是就是stranded?谢谢!

    1. 保留了极性和序列有方向不能提供任何有效的信息。你去问问样品制备人员吧,他们应该清楚自己的样品是不是strand specific的。

  26. 他自己就是置备样本的吧。他说I’m not sure. One thing with the library prep I use though is that it retains polarity, so there is directionality in the sequences.
    但是我看tophat 有个选项在cufflinks里
    option: fr-secondstrand
    example: Directional Illumina, Standard SOLiD
    是不是要用这个呢?我糊涂了?
    这个strand到底是啥意思啊?我如果没有选正确的选项(比如用的TOPHAT/Cufflinks的unstrand是不是很不好啊 )

    1. 你查看一下Comprehensive comparative analysis of strand-specific RNA sequencing methods这篇文章,估计你就明白strand是啥意思了。从他的原话来看,他这个并非strand specific的样品。

  27. 谢谢。我正在看,可是我生物背景不足很多地方看不懂。如果您能简要的说一下那就更好了。反正我的序列是从proton仪器出来的Single-ended的序列。我在网上看到说Directional RNA就是strand specific的,但您的意思是不是说他这意思其实不是directional的 ?因为他是测序的,我问他自己也不清楚。我现在是需要设置参数的软件(tophat/cufflinks/htseq-count)的相关初步分析(mapping , making counts table等)都做了两遍(strand/unstrand),只是不知道该给他哪个结果。

    1. 对这个是关键。其实测序都是有方向的(从5’至3’)。在建库的时候区别了sense和antisense,基本上就是stranded的了。但实际上,通常stranded的测序,也只能达到80%左右的区分。不过这已经足以区分重叠部分的基因了。

  28. 博主我问题又来了。我发现我用tophat(fr-firststrand)得到的结果,和我用tophat(默认是unstranded)得到的结果进行比较。我目前没有replicates,每个condition只有一个样本我不好做统计比较,我只好先做两个样本间mapping reads的fold-change比较。 我发现,fold-change 最高的50个基因(我认为也就是所谓的TOP50 up-regualted的吧)和最低的50个基因(我认为所谓的TOP50 down-regualted的)在两种情况下(fr-firststrand vs unstranded)竟然几乎一样,只有两三个不同。请问这是什么情况?是不是表明库怎么构建无所谓?还是说我分析的步骤有问题?请详解。

  29. 补上条,mapping的reads也差不多。strand参数的甚至比unstrand参数map上的还少一点点(但mapping rate一样),请问这是咋回事?(测序员说他在建库的时候区别了sense和antisense)我预想的话参数不同应该变化蛮大吧。

    1. 我认为这是正常的。使用strand会比unstrand结果少一些。因为有些map到antisense上的reads会被mask掉。但是它们并不会占太大的比率。但其实我并没有比较过使用不同参数map的结果。只是想象的。如果对方都说了是区分了strand specific的,你就按照对方的要求做了就好了。

  30. 博主几个rna-seq构建库的问题问问,谢谢!
    1 区分正义子反义子有意义么?有何好处(据说stranded 比unstranded贵,想知道贵在哪里?)
    2 read length是不是越长越好(也是听说构建库时越长越贵?) 如果不是有没有什么比较合理的最大值?如果是,也是问问为啥贵?
    3 发现价格好像测2 泳道的比1泳道贵,看seqanswer上说好像reads的数目比reads长度更重要(同等条件下),比如one lane of 50bp SE vs. two lanes of 25bp SE, 最好选后者。这是为什么?
    4 除了reads 长度,在构建库时候还有什么需要注意的?我看好像有人说insert-size?(这个好像是针对pair ended 的是吧)?
    Thank you!

    1. 1. 区分正反义子,我想最大的好处就是可以区分重叠的外显子了。其它的好处比如更干净的数据啊之类,其实并不明显。
      2. read length并不是越长越好,要看你做什么。其实长度为50或者100的reads是比较常见的。太长的reads除了增加计算的负担,或者重测序,并不是一般获得表达图谱的理想选择。
      3. 对于数目和长度,我的观点都是够用就行。有个前提是数据干净。如果背景高,再多的reads也是无用的。如果reads数目达到深度测序的要求(比如说基因表达谱4~10X),reads长度超过30就够用了。所谓的one lane of 50bp SE vs. two lanes of 25bp SE这无法比较。应该说是100M 50bp SE reads vs 200M 25bp SE reads,如果是这两者的话,不好说,因为长度太短的reads很可能map到多个位置。我们假设如果长度超过21的reads与更长的reads map到多个位置的可能性是一样的,那么当然是reads数高的好些。
      4. 在构建库的时候,最关键是library的质量,少些背景比什么都强。有人说insert-size这个,可能是他做库的时候做了环化吧。如果是那样的话,insert-size是很重要的。

  31. 接上问,请问replicates选多少(我们知道从统计学角度more repilcates应该是越好是吧但很多时候做实验的人缺钱只能做2-5个)为好?可以从cuffdiff/edgeR/DESeq的结果中看出么?谢谢!

  32. 1.我想最大的好处就是可以区分重叠的外显子了
    这对最终的基因表达水平有帮助吗?
    2.什么叫数据干净?就是背景噪音小?
    3.关于replicates选多少(我们知道从统计学角度more repilcates应该是越好是吧但很多时候做实验的人缺钱只能做2-5个)为好?可以从cuffdiff/edgeR/DESeq的结果中看出么?您说的3个好像是前人的经验我想说有没有什么判定方法还是就是要靠经验?

    1. 1.没有太大帮助。
      2.对。
      3.是经验。主要还是因为条件限制,生物实验成本高,做得多负担不起。3个的结果基本上可能接受。

    1. 你使用boxplot做为关键词,搜索一下图像,你就比较清楚什么是outlier了。通常我们认为是差异显著,但也分情况,要看相比较的两者是不是至少有一方的表达值真正达到了可确认的范围。如果两者都很低,但是比值很高,没有太大意义。

  33. 楼主您好!
    我这边也做了差异表达的RNA测序,做了2个重复,公司帮我分析了数据,我发现2个重复的结果并不一致,很多基因,在2个重复中的表达截然相反,一个上调,一个下调。请问,这是为什么?这2个重复我该如何使用?直接舍弃一个重复,只用其中一个的结果进行后续分析可以吗?
    谢谢!

    1. 我建议你再做一组重复。如果你现在只有2个重复的情况下还要放弃一组,我想最终你的文章质量会受到比较严重的影响。

  34. 主要是我是初学者啊。博主老师的内容储备很强大;但美中不足的是这个博客留言回复不像Seqanwser那样有email提醒你回复通知而且可以显示博主是回复的谁的留言。我在想要是博主能留下自己的email就更好了,哈哈。

  35. 博主我看到一篇文章说RNA-seq比起dna微阵列的好处在于 一个是噪音小 一个是基因表达值的dynamic range更大了,请问这是个好处么?我怎么觉得应该是越小越好啊?能否给点解释?

    1. 可能是你对dynamic range的理解出现了一点偏差。所谓dynamic range就是底物浓度与荧光信号成一次线性关系的区间。对于基因芯片,当样品中要侦测的生物分子的浓度小于一定值,比如说,小于2pm时,所有的值都会变成与背景噪音相近,无法分辨,而小于10pm时,要侦测的底物浓度与荧光信号不是正比例的关系,而是某种二次曲线。当达到饱和浓度时,也是一样。在这些区域,如果你依然按照一次线性关系去推算底物浓度,显然与真实值有很大的偏差。所以,dynamic range就是用于描述可以使用一次线性关系表述底物浓度的区间。

  36. 我看一篇论文介绍RNA-SEQ除了可以发现差异表达之外还有其它应用:比如
    Alternative expression analysis
    Transcript discovery and annotation
    Fusion detection
    能简要的和我说说这三个么?有啥用?我记得您上次说主要(或许80%以上都是对DE genes感兴趣)但我也想知道这其它三个

    1. 我说过80%吗?这个数据应该不准确。不过就是大约的意思。alternative splicing分析就是后剪切分析。一个基因会有多种不同的转录本(transcripts),这些转录本由不同的exon组成。后剪切分析就是分析这个的。transcript discovery就是分析新的转录本。通常我们只是筛选已知基因。但有些时候,也会对未知基因,或者已知基因的未知转录本感兴趣。最后一个fusion detection我没有做过,这应该是进化方向会感兴趣的。我不做进化,所以无法解释了。

  37. 我现在做homo_sapiens的RNA-seq分析,准备用tophat去mapping,但是对软件参数的数值选择倍感困惑,因为我想着测序数据相对实验数据重复比较难,如果参数没选好,对后面的分析是不是有很大影响?比如–read-mismatches , –read-gap-length,–read-edit-dist,这些都是默认值为2的options,那么增加1或者减少1,区别在哪里?已知–mate-inner-dist的值,要如何计算–mate-std-dev的值呢?或者都可以直接使用默认值,因为manual里面说这个值可以不用那么精确?我要怎么去考虑这些参数的改变对数据分析的影响呢?tophat的manual我有仔细读过,但还是不懂该如何下手设置各个参数?或者一般情况下,可以直接参考博主文章里面的提到的tophat设置?

    1. –read-mismatches , –read-gap-length,–read-edit-dist这些参数使用默认值即可,如果mapping的效果特别差才做调整。你要了解这些参数的意义,先读一下序列比对那点事儿。对于mate-std-dev的值,使用默认吧。除非你自己做的实验,清楚的知道是多少。对于tophat,你先试着做一次,主要注意-mate-inner-dist, –library-type即可。其余的在有特殊要求时具体调整。

  38. 博主,我现在有一个基因list,是牛的,有ensembl id和gene name,想转化成人的ensembel id 和gene name。我现在可以用biomart的在线做,但每次都要粘贴(或者upload txt文档),然后鼠标点阿点的。这样如果有10个list我要手动操作十次,很麻烦。请问有没有R的package?好像有个 biomaRt可以做到这个么?

  39. 接上条,我打开biomaRt的mannual,好像没找到相关指令。博主应该很擅长R/bioconductor,不知可否帮忙?
    我想做的事情可以用在线完成
    http://www.ensembl.org/biomart/martview
    我想做的你懂的,主要有些基因牛的基因名大家都不熟但换成人的就都熟了

    1. biomaRt是R/Bioconductor中非常基础的一个组成部分。步骤有三:useMart来指定使用的mart数据库;listFilters and listAttributes来了解可能注释的内容;getBM来得到注释。
      但是你的工作不应该直接使用biomaRt。我想更专业的应该使用inparanoid数据库

  40. 首先我的工作可以用这个网页完成吧?我不确定我这样做对不对请您指教。
    http://www.ensembl.org/biomart/martview
    一个dataset选bos,再输入牛的基因名(ensembl ID),然后另一个dataset选homo,然后点击result。输出基因的bos的ensemble ID ,还有对应的human的ensemble ID和gene name。
    对么?

    其次我想用最好能用R的package做或者什么bash下的软件包,就是自动的不用每次都这样GUI点来点去的。这样我一旦有一个基因的list,比如txt文件就可以在服务器上操作,几条command解决问题。您说的这个inparanoid好像也是得界面操作的吧(而且只能一个基因一个基因的查询)?

  41. 您上次说如果我不想发现新的genes/transcripts,在texudo里可以省去cufflinks和cuffmerge直接用cuffdiff是这样吧?我发现经过cufflinks和cuffmerge后好像gtf的注释会发生变化。

  42. 有点困惑。
    比如我现在在用texudo做牛(Bovine)的RNA-SEQ。Tophat后就开始组装和差分分析了。我如果不用cufflinks+cuffmerge(我还没试没用,过不用cufflinks-cuffmerge是啥样,但我觉得应该会和DEseq/edgeR输出的基因名/id差不多),我估计Cuffdiff的输出应该是Ensembl ID(因为我用的是Ensembl数据库的genes.gtf注释文件,我查看了一下里面有Ensembl Id)。
    但我已经试了用cufflinks+cuffmerge后再cuffdiff后的输出结果了。出人意料的是,由于cuffmerge后的merged.gtf把Ensembl ID 改成gene name了,所以输出的significant=yes的基因就很多是基因名了(只能是gene name而没有ensembl ID)。按说是好事,可是有的重要基因没有名(这个我能理解,我想可能是有的Ensembl ID没有对应的gene name),而有的基因名并不是标准的名字,比如(“Q865A2_BOVIN”),这让我很困惑。因为这个既不是Ensembl id也不是official gene name(好像是Uniprot name)。可不可以有方法统一一下?比如都是ensembl ID 我还可以用biomart去转换?或者都是official的gene name也行。

  43. 请问具体怎么做?有没有在线的工具或者包?还是或者我换个Gtf(不用Ensembl的库而是用比如UCSC的?)
    我现在就是cuffidff的结果总有Q865A2_BOVIN这种不知道该怎么转换?

  44. 接上条: 此外我如果有四种condition,我的merged.gtf可以把四种情况下的transcripts都aseemble到一起然后用这个四种情况下的merged.gtf而不用cuffdiff比较的两两情况merge是吧?
    比如C1vsC2,C2vsC3都用的是merged.gtf(cuffmerge C1,C2,C3的transcripts.gtf生成的),对吧?

  45. 博主,我发现牛reference genome的genes.gtf注释文件很有问题(我选的是ENSEMBL库的bos4.0库)。
    它的“gene name” 有的是gene name,有的则是uniprot ID。很讨厌。

    这就导致如果用texudo会出现这样的结果
    test_id gene_id gene locus sample_1 sample_2 status value_1 value_2 log2(fold_change) test_stat p_value q_value significant
    XLOC_022293 XLOC_022293 THIOM_BOVIN 5:80512046-80522058 C1 C2 OK 59.4749 52.7758 -0.172407 -1.10325 0.09935 0.222418 no
    XLOC_008347 XLOC_008347 FOXF1 18:11584384-11588769 C1 C2 OK 1.52995 2.89378 0.919468 3.87032 5e-05 0.000369275 yes

    看到了么?第三列并不统一.
    我个人认为解决方法就是全部输出Ensembl ID,然后统一转换为gene name。这样不仅方便也可以和tophat-htseq-edgeR/DESeq的结果直接进行比较。因为由于htseq用的是genes.gtf,所以输出也是Ensembl ID,这样和texudo计算DE genes的overlap也方便点。

    可现在问题是,您说过可以略去cufflinks和cuffmerge。我发现如果跳过cufflinks和cuffmerge,就可以生成Ensembl ID,但结果却和不跳过的不太一样。

    请对比刚才的
    test_id gene_id gene locus sample_1 sample_2 status value_1 value_2 log2(fold_change) test_stat p_value q_value significant
    ENSBTAG00000000014 ENSBTAG00000000014 THIOM_BOVIN 5:80512046-80522058 C1 C2 OK 78.8215 68.6068 -0.200238 -1.32214 0.0649 0.140941 no
    ENSBTAG00000000009 ENSBTAG00000000009 FOXF1 18:11584384-11587706 C1 C2 OK 2.38301 4.11073 0.786607 4.44249 0.00015 0.000820208 yes

    请注意 主要是p-value,q-value和fold change 不一样。但是是否significant的结果是一样的。
    请问这要紧么?
    如果要紧有没有其他方法?(比如不略去cufflinks+cuffmerge但依然生成Ensembl ID的文件?)

    1. 你的cuffmerge的结果中的gene_name一栏中是什么?如果你希望它是Ensembl ID的话,那一栏中应该是Ensembl ID才好。或者就是按我说的,使用cummeRbund+biomaRt来对结果再注释。

  46. 您好,我的cuffmerge的gene_name一栏就是gene.gtf里的gene name, 也就是上一条展示给您的有的是gene name,有的则是uniprot ID。很讨厌。(这个文件是Ensembl 最新的Bos4.0 genome file)
    我现在几个问题
    1.为什么有cufflinks+cuffmerge的和没有的结果不一样?(虽然yes和no一样但p-value 和 fold change不一样)这正常么?
    2. 如何用cummRbund+biomaRt?我记得cummRbund只是一个可视化程序包啊。无非就是能输出一个你想输出的列(可以是gene name 也可以是 gene id)。您的意思是用cummRbund输出一个列然后输入biomaRT在转换?还是cummRbund自己有什么功能能转换?
    3.除了biomaRT,您以前推荐我的一个org.Bt.eg.db可以做ID conversion么?

    1. 1.因为加入计数的reads总数不同,基因的长度也可能不同。
      2.我的意思是你输出genomic location,然后再注释得到ensembl ID
      3.org.Bt.eg.db也可以,但从genomic location至ID的转换可能会比较困难一点。
      4.如果使用biomaRt遇到上述提示,考虑两点,你的网络不稳定,你没有联接至网络。

发表评论

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