[转] cufflinks介绍(2018-05-29)

原文:http://blog.sina.com.cn/s/blog_751bd9440102v72b.html

一. 简介

Cufflinks下主要包含cufflinks,cuffmerge,cuffcompare和cuffdiff等几支主要的程序。主要用于基因表达量的计算和差异表达基因的寻找。

二. 安装

Cufflinks下载网页

1. 为了安装Cufflinks,必须有Boost C++

libraries。下载Boost并安装。默认安装在/usr/local。

$ tar jxvf boost_1_53_0.tar.bz2

$ cd boost_1_53_0

$ ./bootstrap.sh

$ sudo ./b2 install

2.安装SAM tools。

下载SAM tools。$ tar jxvf samtools-0.1.18.tar.bz2$ cd samtools-0.1.18$ make$ sudo su # mkdir /usr/local/include/bam# cp libbam.a /usr/local/lib# cp *.h /usr/local/include/bam/# cp samtools /usr/bin/

3. 安装 Eigen libraries。

下载Eigen$ tar jxvf 3.1.2.tar.bz2$ cd eigen-eigen-5097c01bcdc4$ sudo cp -r Eigen/ /usr/local/include/

4. 安装Cufflinks。

$ tar zxvf cufflinks-2.0.2.tar.gz

$ cd cufflinks-2.0.2

$ ./configure --prefix=/path/to/cufflinks/install --with-boost=/usr/local/ --with-eigen=/usr/local/include//Eigen/

$ make

$ make install

5. 可以直接下载Linux x86_64 binary。不需要上述繁琐步骤,解压后的程序直接可用。(推荐)

三.Cufflinks的使用

1. Cufflinks简介

Cufflinks程序主要根据Tophat的比对结果,依托或不依托于参考基因组的GTF注释文件,计算出(各个gene的)isoform的FPKM值,并给出trascripts.gtf注释结果(组装出转录组)。

注意:

1.

fragment的长度的估测,若为pair-end测序,则cufflinks自己会有一套算法,算出结果。若为single-end测序,则cufflinks默认的是高斯分布,或者你自己提供相关的参数设置。

2. cufflinks计算multi-mapped reads,一般a read map到10个位置,则每个位置记为10%。a

read mapping to 10 positions will count as 10% of a read at each

position.

3. 一般不推荐用cufflinks拼接细菌的转录组,推荐Glimmer。但是,若有注释文件,可以用cufflinks和cuffdiff来检测基因的表达和差异性。

4. cufflinks/cuffdiff不能计算出exon或splicing event的FPKM

5.cuffdiff处理时间序列data:采用参数-t

6.当你使用cufflinks时,在最后出现了99%,然后一直不动。因为cuffdiff需要更多的CPU来处理一些匹配很多reads的loci。而这些位点一般要等其他位点全部解决了后,才由cuffdiff来处理。可以用参数-M来提供相关的文件,过滤掉rRNA或者线粒体RNA。

7. 当使用cufflinks或cuffdiff出现了“crash with a ‘bad_alloc'

error”,cuffdiff和cufflinks运行了很长时间才结束————这表明计算机拼接一个高表达的基因或定量分析一个高表达的基因,运行的内存使用玩尽了!解决方法:修改选项“-max-bundle-frags”,可以先尝试500000,若错误依旧在,可以继续下调!

8. cuffdiff报道的结果里面所有的基因和转录本的FPKM=0,这表明GTF中的染色体名字和BAM里的名字不匹配。

9.cuffdiff和cufflinks的缺点:存在一定的假基因和转录本(原因:测序深度,测序质量,测序样本的测序次数,以及注释的错误)

10. large fold

change表达量不代表数据的明显性(这些基因的isform多或这些基因测序测到的少,整体较低的表达)。cuffdiff中明显表达倍数改变的基因,存在不确定性。

11.通过cufflinks产生的结果中transcript.gtf文件中cuff标识的转录本就是新的转录本。相应的,其他模块输出中CUFF标识代表着新的转录本。

12. 若出现了如下错误:

YouareusingCufflinksv2.2.1,whichisthemostrecentrelease.

open:Nosuchfileordirectory

File30doesn'tappeartobeavalidBAMfile,tryingSAM...

Error:cannotopenalignmentfile30forreading

这表明,你的参数有问题。例如“--min-intron-length”,你设置为了:“-min-intron-length”

2. 使用方法

$ cufflinks [options]*

一个常用的例子:

$ cufflinks -p 8 -G transcript.gtf --library-type fr-unstranded -o cufflinks_output tophat_out/accepted_hits.bam

3. 普通参数

-h | --help

-o | --output-dir  default: ./设置输出的文件夹名称

-p | --num-threads  default: 1用于比对reads的CPU线程数

-G | --GTF 提供一个GFF文件,以此来计算isoform的表达。此时,将不会组装新的transcripts,程序会忽略和reference transcript不兼容的比对结果

-g | --GTF-guide 提供GFF文件,以此来指导转录子组装(RABT assembly)。此时,输出结果会包含reference transcripts和novel genes and isforms。

-M | --mask-file 提供GFF文件。Cufflinks将忽略比对到该GTF文件的transcripts中的reads。该文件中常常是rRNA的注释,也可以包含线立体和其它希望忽略的transcripts的注释。将这些不需要的RNA去除后,对计算mRNA的表达量是有利的。

-b | --frag-bias-correct 提供一个fasta文件来指导Cufflinks运行新的bias detection and correction algorithm。这样能明显提高转录子丰度计算的精确性。

-u | --multi-read-correct让Cufflinks来做initial estimation步骤,从而更精确衡量比对到genome多个位点的reads。

--library-type  default:fr-unstranded处理的reads具有链特异性。比对结果中将会有个XS标签。一般Illumina数据的library-type为 fr-unstranded。

--library-norm-method具体参考官网,三种方式:classic-fpkm默认的方式。geometric针对DESeq。quartile计算时,fragments和总的map的count取75%

4. 丰度评估参数

-m | --frag-len-mean default: 200插入片段的平均长度。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主设置此值。

-s | --frag-len-std-dev default: 80插入片段长度的标准差。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主设置此值。

-N | --upper-quartile-form使用75%分为数的值来代替总的值(比对到单一位点的fragments的数值),作normalize。这样有利于在低丰度基因和转录子中寻找差异基因。

--total-hits-norm default: TRUECufflinks在计算FPKM时,算入所有的fragments和比对上的reads。和下一个参数对立。默认激活该参数。

--compatible-hits-norm Cufflinks在计算FPKM时,只针对和reference transcripts兼容的fragments以及比对上的reads。该参数默认不激活,只能在有 --GTF 参数下有效,并且作 RABT或 ab initio 的时候无效。

--max-mle-iterations进行极大似然法时选择的迭代次数,默认为:5000

--max-bundle-frags一个skipped locus/loci在别skipped前可以拥有的最大的fragment片段。默认为1000000

--no-effective-length-correctionCufflinks will not employ its "effective" length normalization to transcript FPKM.Cufflinks将不会使用它的“effective” 长度标准化去计算转录的FPKM

--no-length-correctionCufflinks将根本不会使用转录本的长度去标准化fragment的数目。当fragment的数目和the features being quantified的size是独立的,可以使用(例如for small RNA libraries, where no fragmentation takes place, or 3 prime end sequencing, where sampled RNA fragments are all essentially the same length).小心使用

5. 组装常用参数

-L | --label  default: CUFFCufflink以GTF格式来报告转录子片段(transfrags),该参数是GTF文件的前缀

-F/--min-isoform-fraction <0.0-1.0>在计算一个基因的isoform 丰度后,过滤了丰度极低的转录本,因为这些转录本不可以信任。也可以过滤一些read匹配极低的外显子。默认为0.1或者10% of the most abundant isoform (the major isoform) of the gene.(一个基因的主要isoform的丰度的10%)

-j/--pre-mrna-fraction <0.0-1.0>内含子被aligment覆盖的最低深度。若小于这个值则那些内含子的alignments被忽略掉。默认为15%。 The minimum depth of coverage in the intronic region covered      by the alignment is divided by the number of spliced reads, and if the          result is lower than this parameter value, the intronic alignments are          ignored. The default is 15%.

-I/--max-intron-length内含子的最大长度。若大于该值的内含子,cufflinks不会报告。默认为300000.Cufflinks will not report transcripts with    introns longer than this, and will ignore SAM alignments with REF_SKIP          CIGAR operations longer than this.  The default is 300,000.

-a/--junc-alpha <0.0-1.0>剪接比对过滤中假阳性的二项检验中的 alpha value。默认为 0.001

-A/--small-anchor-fraction <0.0-1.0>在junction中一个reads小于自身长度的这个百分比,会被怀疑,可能会在拼接前被过滤掉。默认为0.09

--min-frags-per-transfrag  default: 10组装出的transfrags被支持的RNA-seq的fragments数少于该值则不被报道。

--overhang-tolerance当决定一个reads或转录本与某个转录本兼容或匹配的时候,允许的能加入该转录本的外显子的延伸长度。默认是8bp和bowtie/tophat默认的一致。

--max-bundle-lengthMaximum genomic length allowed for a given bundle.  The default is 3,500,000bp.

--min-intron-length  default: 50最小的intron大小。

--trim-3-avgcov-thresh最小的3‘端的平均覆盖程度。小于该值,则删除其3’端序列。默认10Minimum average coverage required to attempt 3' trimming.  The default is 10.

--trim-3-dropoff-frac最低百分比的拼接的转录本的3‘端的平均覆盖程度。默认0.1The fraction of average coverage below which to trim the 3' end of an assembled          transcript.  The default is 0.1.

--max-multiread-fraction <0.0-1.0>若一个转录本Transfrags的reads能匹配到基因组的多个位置,其中该转录本的reads有超过该百分比是multireads,则不会报告这个转录本。默认为75%The fraction a transfrag's supporting reads that may be multiply mapped to the genome. A transcript composed of more than this fraction will not be reported by the assembler.  Default: 0.75 (75% multireads or more is suppressed).

--overlap-radius  default: 50Transfrags之间的距离少于该值,则将其连到一起。

Advanced Reference Annotation Based Transcript (RABT) Assembly Options:当你使用-g/--GTF-guide这个参数时,需要考虑的选项。

--3-overhang-tolerance当决定一个拼接的转录本(这个转录本可能不是新的转录本)和一个参考转录本是否合并时,参考转录本的3‘端允许延伸的长度。默认600bpThe number of bp allowed to overhang the 3' end of a reference transcript when determining      if an assembled transcript should be merged with it (ie, the assembled transcript is not novel).        The default is 600 bp.

--intron-overhang-tolerance当决定一个拼接的转录本(这个转录本可能不是新的转录本)和一个参考转录本是否合并时,参考转录本的外显子允许延伸的长度。默认50bpThe number of bp allowed to enter the intron of a reference transcript when determining if an    assembled transcript should be merged with it (ie, the assembled transcript is not novel).      The default is 50 bp.

--no-faux-readsThis option disables tiling of the reference transcripts with faux reads.  Use this if you only        want to use sequencing reads in assembly but do not want to output assembled transcripts that lay      within reference transcripts.  All reference transcripts in the input annotation will also      be included in the output.这一项将不能掩盖参考转录组中的假reads。当你只想在拼接中使用测序的reads而不想输出lay within reference transcripts的拼接的转录组。输入时注释的所有的参考转录组也将会输入到输出中。

其他参数(无关紧要)

-v/--verbose显示版本信息等等

-q/--quiet除了警告和错误外,其他信息将不会print

--no-update-check关系cufflinks自动更新的能力

6. Cufflinks输出结果

cufflinks的输入文件是sam或bam格式。并且sam或bam格式的文件必须排好序。(The SAM file supplied to Cufflinksmustbe sorted by          reference position.)Tophat的输出结果sam或bam已经排好了序。针对其他的未排序的sam或bam文件采用如下排序方式:

sort -k 3,3 -k 4,4n hits.sam > hits.sam.sorted

1. transcripts.gtf

该文件包含Cufflinks的组装结果isoforms。前7列为标准的GTF格式,最后一列为attributes。其每一列的意义:

列数  列的名称  例子        描述

1    序列名    chrX        染色体或contig名; 2    来源      Cufflinks  产生该文件的程序名; 3    类型      exon        记录的类型,一般是transcript或exon; 4    起始      1          1-base的值; 5    结束      1000        结束位置; 6    得分      1000        ; 7    链        +          Cufflinks猜测isoform来自参考序列的那一条链,一般是'+','-'或'.';8    frame    .          Cufflinks不去预测起始或终止密码子框的位置; 9    attributes  ...      详见下

每一个GTF记录包含如下attributes:

Attribute      例子      描述

gene_idCUFF.1Cufflinks的gene id;transcript_idCUFF.1.1  Cufflinks的转录子 id; FPKM          101.267  isoform水平上的丰度,FragmentsPerKilobase of exon model perMillion mapped fragments; frac          0.7647    保留着的一项,忽略即可,以后可能会取消这个;conf_lo        0.07      isoform丰度的95%置信区间的下边界,即 下边界值 = FPKM * ( 1.0 - conf_lo );conf_hi        0.1102    isoform丰度的95%置信区间的上边界,即 上边界值 = FPKM * ( 1.0 + conf_hi ); cov            100.765  计算整个transcript上read的覆盖度;full_read_support  yes  当使用 RABT assembly 时,该选项报告所有的introns和exons是否完全被reads所覆盖

2. ispforms.fpkm_tracking

isoforms(可以理解为gene的各个外显子)的fpkm计算结果

3. genes.fpkm_tracking

gene的fpkm计算结果

四.Cuffmerge的使用

1. Cuffmerge简介

Cuffmerge将各个Cufflinks生成的transcripts.gtf文件融合称为一个更加全面的transcripts注释结果文件merged.gtf。以利于用Cuffdiff来分析基因差异表达。

2. 使用方法

$ cuffmerge [options]*

输入文件为一个文本文件,是包含着GTF文件路径的list。常用例子:

$ cuffmerge -o ./merged_asm -p 8 assembly_list.txt

3. 使用参数

-h | --help

-o  default: ./merged_asm

将结果输出至该文件夹。

-g | --ref-gtf将该reference GTF一起融合到最终结果中。

-p | --num-threads  defautl: 1

使用的CPU线程数

-s | --ref-sequence /该参数指向基因组DNA序列。如果是一个文件夹,则每个contig则是一个fasta文件;如果是一个fasta文件,则所有的contigs都需要在里面。Cuffmerge将使用该ref-sequence来帮助对transfrags分类,并排除repeats。比如transcripts包含一些小写碱基的将归类到repeats.

4. Cuffmerge输出结果

输出的结果文件默认为 /merged.gtf

五.Cuffcompare的使用

1. Cuffcompare简介

Cuffcompare使用Cufflinks的GTF结果,对GTF结果进行比较。和reference gtf比较寻找novel转录本等。

2. Cuffcompare的使用方法

$ cuffcompare [options]*  [cuff2.gtf] ... [cuffN.gtf]

使用例子:

$ cuffcompare -o cuffcmp cuff1.gtf cuff2.gtf

3. 使用参数

-h-V显示进程

-C默认,表示"contained" transcripts 也会写入.combined.gtf中。

-o  default: cuffcmp输出文件的前缀

-r 参考的GFF文件。用来评估输入的gtf文件中gene models的精确性。每一个输入的gtf的isoforms将和该参考文件进行比较,并被标注为 overlapping, matching 或 novel。

-R当有了 -r 参数时,指定该参数时,将忽略参考GFF文件中的一些transcripts。这些transcripts不和任何输入的GTF文件overlapped。

-s

该参数指向基因组DNA序列。如果是一个文件夹,则每个contig则是一个fasta文件;如果是一个fasta文件,则所有的contigs都需要在里面。小写字母的碱基用来将相应的transcripts作为repeats处理。

4.输出结果

在当前目录下输出3个文件:

.stats, 报告与参考注释比较时,各种与准确性相关的数据。其中,Sn和Sp展示的是specificity and sensitivity values。fSnandfSp列展示的 "fuzzy" variants of these same accuracy calculations。允许存在变动。(-o 没有设置,默认为cuffcmp为文件前缀)

.combined.gtf报告每个样本的所有的 transfrags 的信息。若一个transfrag在多个样本中,它只报道一次。

.tracking匹配到样本间的转录本。this file matches transcripts up between samples.  Each row contains                a transcript structure that is present in one or more input GTF files.                Because the transcripts will generally have different IDs (unless you                assembled your RNA-Seq reads against a reference transcriptome),cuffcompareexamines the structure of each the transcripts,                matching transcripts that agree on the coordinates and order of all of                their introns, as well as strand.  Matching transcripts are allowed to                differ on the length of the first and last exons, since these lengths                will naturally vary from sample to sample due to the random nature of                sequencing.

例子;

TCONS_00000045 XLOC_000023 Tcea|uc007afj.1    j      \    q1:exp.115|exp.115.0|100|3.061355|0.350242|0.350207 \    q2:60hr.292|60hr.292.0|100|4.094084|0.000000|0.000000

In this example, a transcript present in the two input files,calledexp.115.0in the first and60hr.292.0inthe second, doesn't match any reference transcript exactly, butshares exons withuc007afj.1, an isoform of the gene Tcea,as indicated by theclass

codej. The first three columns are as follows:

其中,1 Cufflinks transfrag idTCONS_00000045内部的transfrag id;2Cufflinks locus idXLOC_000023内部的locus id; 3Reference gene idTcea参考的注释的gene的id或者“-”表示没有匹配到参考的转录本; 4Reference transcript iduc007afj.1参考的注释的转录本的id或者“-”表示没有匹配到参考的转录本; 5 Class codec转录本和参考转录本之间的匹配类型。第五列之后如下:

qJ: | | | | | | |

在输入的GTF的同目录下输出.refmap 和 .tmap 文件。

.refmap具体内容如下:

1Reference gene name参考注释的gtf中的基因名字 2 Reference transcript id 参考的转录本id3Class code 表示cufflinks拼接的转录本和参考转录本间的匹配情况:c 表示部分匹配;= 表示全部匹配

4Cufflinks matches匹配到参考转录本的cufflinks拼接的转录本的id

.tmap具体内容如下:

1Reference gene name参考注释的gtf中的基因名字 2 Reference transcript id 参考的转录本id3Class code 表示cufflinks拼接的转录本和参考转录本间的匹配情况:c 表示部分匹配;= 表示全部匹配

4 Cufflinks gene id; 5 Cufflinks transcript id;6 Fraction of major isofor m (FMI) ; 7FPKM ; 8 FPKM_conf_lo; 9FPKM_conf_hi; 10 Coverage ; 11 Length; 12Major isoform ID

class cord :

PriorityCodeDescription

1=Complete match of intron chain

2cContained

3jPotentially novel isoform (fragment): at least one

splice junction is shared with a reference transcript

4eSingle exon transfrag overlapping a reference exon

and at least 10 bp of a reference intron, indicating a possible

pre-mRNA fragment.

5iA transfrag falling entirely within a reference

intron

6oGeneric exonic overlap with a reference

transcript

7pPossible polymerase run-on fragment (within

2Kbases of a reference transcript)

8rRepeat. Currently determined by looking at the

soft-masked reference sequence and applied to transcripts where at

least 50% of the bases are lower case

9uUnknown, intergenic transcript

10xExonic overlap with reference on the opposite

strand

11sAn intron of the transfrag overlaps a reference

intron on the opposite strand (likely due to read mapping

errors)

12.(.tracking file only, indicates multiple

classifications)

六.Cuffdiff的使用

1. Cuffdiff简介

用于寻找转录子表达的显著性差异。

2. Cuffdiff使用方法

cuffdiff主要是发现转录本表达,剪接,启动子使用的明显变化。

cuffdiff [options]* ...

[sampleN.sam_replicate1.sam[,...,sample2_replicateM.sam]]

$ cuffdiff [options]*  ...[sampleN_1.sam[,...,sampleN_M.sam]]其中transcripts.gtf是由cufflinks,cuffcompare,cuffmerge所生成的文件,或是由其它程序生成的。一个样本有多个replicate,用逗号隔开。sample多于一个时,cuffdiff将比较samples间的基因表达的差异性。一个常用例子:$ cuffdiff --lables lable1,lable2 -p 8 --time-series --multi-read-correct --library-type fr-unstranded --poisson-dispersion transcripts.gtf sample1.sam sample2.sam

cuffdiff接受bam/sam或cuffquant的CXB文件,同时也可以接受bam与sam的混合文件,不能接受bam/sam和CXB的混合文件。

3. 使用参数

-h | --help

-o | --output-dir  default: ./

输出的文件夹目录。

-L | --lables  default: q1,q2,...qN

给每个sample一个样品名或者一个环境条件一个lable

-p | --num-threads  default: 1

使用的CPU线程数

-T | --time-series

让Cuffdiff来按样品顺序来比对样品,而不是对所有的samples都进行两两比对。即第二个SAM和第一个SAM比;第三个SAM和第二个SAM比;第四个SAM和第三个SAM比...

-N | --upper-quartile-form

使用75%分为数的值来代替总的值(比对到单一位点的fragments的数值),作normalize。这样有利于在低丰度基因和转录子中寻找差异基因。

--total-hits-norm

Cufflinks在计算FPKM时,算入所有的fragments和比对上的reads。和下一个参数对立。默认不激活该参数。

--compatible-hits-normCufflinks在计算FPKM时,只针对和reference transcripts兼容的fragments以及比对上的reads。该参数默认激活,使用该参数可以降低核糖体rna的reads对基因表达的干扰。

-b | --frag-bias-correct(一般是genome.fa)提供一个fasta文件来指导Cufflinks运行新的bias detection and correction algorithm。这样能明显提高转录子丰度计算的精确性。

-u | --multi-read-correct让Cufflinks来做initial estimation步骤,从而更精确衡量比对到genome多个位点的reads。

-c | --min-alignment-count  default: 10

如果比对到某一个位点的fragments数目少于该值,则不做该位点的显著性分析。认为该位点的表达量没有显著性差异。

-M | --mask-file

提供GFF文件。Cufflinks将忽略比对到该GTF文件的transcripts中的reads。该文件中常常是rRNA的注释,也可以包含线立体和其它希望忽略的transcripts的注释。将这些不需要的RNA去除后,对计算mRNA的表达量是有利的。

-FDR  default: 0.05允许的false discovery rate.

--library-type default:fr-unstranded处理的reads具有链特异性。比对结果中将会有个XS标签。一般Illumina数据的library-type为 fr-unstranded。

--dispersion-method

其他高级参数:

-m | --frag-len-mean default: 200插入片段的平均长度。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主设置此值。

-s | --frag-len-std-dev default: 80

插入片段长度的标准差。不过现在Cufflinks能learns插入片段的平均长度,因此不推荐自主设置此值。-v/--verbose显示版本信息等等

-q/--quiet除了警告和错误外,其他信息将不会print

--no-update-check关系cufflinks自动更新的能力

-F/--min-isoform-fraction <0.0-1.0>建议不要更改,主要的isorform丰度若低于这个分数,可变的isoform将四舍五入为0.默认为1e-5

--max-bundle-frags一个skipped locus/loci在skipped前可以拥有的最大的fragment片段。默认为1000000

--max-frag-count-draws (默认为100)和--max-frag-assign-draws (默认为50)--min-reps-for-js-test一个针对不同调控的基因做test的最小的复制次数。Cuffdiff won't test genes for differential regulation unless the

conditions in question have at least this many replicates.  Default: 3.

--no-effective-length-correctionCuffdiff will not employ its "effective" length normalization to transcript FPKM. Cufflinks将不会使用它的“effective” 长度标准化去计算转录的FPKM

--no-length-correctioncufflinks将根本不会使用转录本的长度去标准化fragment的数目。当fragment的数目和the

features being quantified的size是独立的,可以使用(例如for small RNA libraries,

where no fragmentation takes place, or 3 prime end sequencing, where

sampled RNA fragments are all essentially the same length).小心使用

--max-mle-iterations极大似然法的迭代次数,默认5000

--poisson-dispersionUse the Poisson fragment dispersion model instead of learning one in each condition.

4.Cuffdiff输出

1. FPKM tracking filescuffdiff计算每个样本中的转录本,初始转录本和基因的FPKM。其中,基因和初始转录本的FPKM的计算是在每个转录本group和基因group中的转录本的FPKM的求和。

isoforms.fpkm_trackingTranscript FPKMs

genes.fpkm_trackingGene FPKMs. Tracks the summed FPKM of transcriptssharing eachgene_id

cds.fpkm_trackingCoding sequence FPKMs. Tracks the summed FPKM oftranscripts sharing eachp_id, independent oftss_id

tss_groups.fpkm_trackingPrimary transcript FPKMs. Tracks the summed FPKMof transcripts sharing eachtss_id

2. Count tracking files评估每个样本中来自每个 transcript, primary transcript,                and gene的fragment数目。其中primary transcript,                and gene的fragment数目是每个primary transcript group或gene group中trancript的数目之和。

isoforms.count_trackingTranscript counts

genes.count_trackingGene counts. Tracks the summed counts oftranscripts sharing eachgene_id

cds.count_trackingCoding sequence counts. Tracks the summed countsof transcripts sharing eachp_id, independent oftss_id

tss_groups.count_trackingPrimary transcript counts. Tracks the summedcounts of transcripts sharing eachtss_id

3. Read group tracking files计算在每个repulate中每个transcript, primary transcript和gene的表达量和frage数目

isoforms.read_group_trackingTranscript read group tracking

genes.read_group_trackingGene read group tracking. Tracks the summedexpression and counts of transcripts sharing eachgene_idin each replicate

cds.read_group_trackingCoding sequence FPKMs. Tracks the summedexpression and counts of transcripts sharing eachp_id,independent oftss_idin each replicate

tss_groups.read_group_trackingPrimary transcript FPKMs. Tracks the summedexpression and counts of transcripts sharing eachtss_idin each replicate

4. Differential expression test对于splicing transcript,                primary transcripts, genes, and coding sequences.样本之间的表达差异检验。对于每一对样本x和y,都会有以下四个文件:

isoform_exp.diffTranscript differential FPKM.

gene_exp.diffGene differential FPKM. Tests difference sin thesummed FPKM of transcripts sharing eachgene_id

tss_group_exp.diffPrimary transcript differential FPKM. Testsdifferences in the summed FPKM of transcripts sharing eachtss_id

cds_exp.diffCoding sequence differential FPKM. Testsdifferences in the summed FPKM of transcripts sharing eachp_idindependent oftss_id

每个文件的样式如下:

Column numberColumn nameExampleDescription

1Tested idXLOC_000001A unique identifier describing the transcipt,

gene, primary transcript, or CDS being tested

2geneLypla1Thegene_name(s) orgene_id(s)being tested

3locuschr1:4797771-4835363Genomic coordinates for easy browsing to the genes

or transcripts being tested.

4sample 1LiverLabel (or number if no labels provided) of the

first sample being tested

5sample 2BrainLabel (or number if no labels provided) of the

second sample being tested

6Test statusNOTESTCan be one of OK (test successful), NOTEST (not

enough alignments for testing), LOWDATA (too complex or shallowly

sequenced), HIDATA (too many fragments in locus), or FAIL, when an

ill-conditioned covariance matrix or other numerical exception

prevents testing.

7FPKMx8.01089FPKM of the gene in samplex

8FPKMy8.551545FPKM of the gene in sampley

9log2(FPKMy/FPKMx)0.06531The (base 2) log of the fold changey/x

10test stat0.860902The value of the test statistic used to compute

significance of the observed change in FPKM

11p value0.389292Theuncorrectedp-value of thetest statistic

12q value0.985216TheFDR-adjustedp-value of thetest statistic

13significantnoCan be either "yes" or "no", depending on whetherpis greater then the FDRafterBenjamini-Hochbergcorrection for multiple-testing

5. Differential splicing tests – splicing.diff对于每个primary transcript,鉴定的不同的isoform的差异性。只有2个或2个以上的isoforms的primary transcript存在

Column numberColumn nameExampleDescription

1Tested idTSS10015A unique identifier describing the primary

transcript being tested.

2gene nameRtknThegene_nameorgene_idthatthe primary transcript being tested belongs to

3locuschr6:83087311-83102572Genomic coordinates for easy browsing to the genes

or transcripts being tested.

4sample 1LiverLabel (or number if no labels provided) of the

first sample being tested

5sample 2BrainLabel (or number if no labels provided) of the

second sample being tested

6Test statusOKCan be one of OK (test successful), NOTEST (not

enough alignments for testing), LOWDATA (too complex or shallowly

sequenced), HIDATA (too many fragments in locus), or FAIL, when an

ill-conditioned covariance matrix or other numerical exception

prevents testing.

7Reserved0

8Reserved0

9√JS(x,y)0.22115The splice overloading of the primary transcript,

as measured by the square root of the Jensen-Shannon divergence

computed on the relative abundances of the splice variants

10test stat0.22115The value of the test statistic used to compute

significance of the observed overloading, equal to √JS(x,y)

11p value0.000174982Theuncorrectedp-value of thetest statistic.

12q value0.985216TheFDR-adjustedp-value of thetest statistic

13significantyesCan be either "yes" or "no", depending on whetherpis greater then the FDRafterBenjamini-Hochbergcorrection for multiple-testing

6. Differential coding output – cds.diff对于每个基因,它的cds的鉴定。样本间的输出cds的差异性。只有2个或2个以上的cds(multi-protein genes)列举在文件中。

Column numberColumn nameExampleDescription

1Tested idXLOC_000002-[chr1:5073200-5152501]A unique identifier describing the gene being

tested.

2gene nameAtp6v1hThegene_nameorgene_id

3locuschr1:5073200-5152501Genomic coordinates for easy browsing to the genes

or transcripts being tested.

4sample 1LiverLabel (or number if no labels provided) of the

first sample being tested

5sample 2BrainLabel (or number if no labels provided) of the

second sample being tested

6Test statusOKCan be one of OK (test successful), NOTEST (not

enough alignments for testing), LOWDATA (too complex or shallowly

sequenced), HIDATA (too many fragments in locus), or FAIL, when an

ill-conditioned covariance matrix or other numerical exception

prevents testing.

7Reserved0

8Reserved0

9√JS(x,y)0.0686517The CDS overloading of the gene, as measured by

the square root of the Jensen-Shannon divergence computed on the

relative abundances of the coding sequences

10test stat0.0686517The value of the test statistic used to compute

significance of the observed overloading, equal to √JS(x,y)

11p value0.00546783Theuncorrectedp-value of thetest statistic

12q value0.985216TheFDR-adjustedp-value of thetest statistic

13significantyesCan be either "yes" or "no", depending on whetherpis greater then the FDRafterBenjamini-Hochbergcorrection for multiple-testing

7. Differential promoter use – promoters.diff样本间启动子使用的差异性。只有表达2个或2个以上isoform的基因列举在这里。

8. Read group info – read_groups.info每个repulate,在进行定量分析时,cuffdiff的关键属性会列出。

Column numberColumn nameExampleDescription

1filemCherry_rep_A/accepted_hits.bamBAM or SAM file containing the data for the read

group

2conditionmCherryCondition to which the read group belongs

3replicate_num0Replicate number of the read group

4total_mass4.72517e+06Total number of fragments for the read group

5norm_mass4.72517e+06Fragment normalization constant used during

calculation of FPKMs.

6internal_scale1.23916Internal scaling factor, used to transform

replicates of a single condition onto the "internal" common count

scale.

7external_scale0.96External scaling factor, used to transform counts

from different conditions onto an internal common count scale.

9. Run info – run.info运行的信息。

其中:输出文件FPKM Tracking file的格式如下:

1tracking_idTCONS_00000001内部唯一object的id(识别基因,转录本,CDS,初始转录本)A

unique identifier describing the object (gene, transcript, CDS,

primary transcript)

2class_code=内部定义的类别的id,“-”表明不是转录本。Theclass_codeattribute for the object, or "-" if not a transcript, or ifclass_codeisn't present

3nearest_ref_idNM_008866.1最接近的参考转录本The

reference transcript to which the class code refers, if

any

4gene_idNM_008866基因idThegene_id(s)

associated with the object

5gene_short_nameLypla1基因名字Thegene_short_name(s)

associated with the object

6tss_idTSS1初始转录本id,或者“-”表示没有初始转录本。Thetss_idassociated with the object, or "-" if not a transcript/primary

transcript, or iftss_idisn't

present

7locuschr1:4797771-4835363基因组上的位置Genomic

coordinates for easy browsing to the object

8length2447转录本的长度The

number of base pairs in the transcript, or '-' if not a

transcript/primary transcript

9coverage43.4279read覆盖深度的估测值Estimate for the absolute depth of read coverage across the

object

10q0_FPKM8.01089样本0中object的FPKMFPKMof

the object in sample 0

11q0_FPKM_lo7.03583object在样本0中FPKM的95%置信区间的下界the

lower bound of the 95% confidence interval on the FPKM of the

object in sample 0

12q0_FPKM_hi8.98595object在样本0中FPKM的95%置信区间的上界the

upper bound of the 95% confidence interval on the FPKM of the

object in sample 0

13q0_statusOKobject在样本0中的量化状态,0K表示成功,LOWDATA:太复杂或测序深度不够;HIDATA:在一个基因座上太多fragments,FAIL:失败的协方差矩阵或其他数值阻止了去卷积Quantification

status for the object in sample 0. Can be one of OK

(deconvolutionsuccessful), LOWDATA (too complex or shallowly sequenced), HIDATA

(too many fragments in locus), or FAIL, when an ill-conditioned

covariance matrix or other numerical exception preventsdeconvolution.

Count tracking files 格式如下:

1tracking_idTCONS_00000001A unique identifierdescribing the object (gene, transcript, CDS, primarytranscript)

2q0_count201.334Estimated (externally scaled) number of fragments generated by theobject in sample 0

3q0_count_variance5988.24Estimated variance inthe number of fragments generated by the object in sample0

4q0_count_uncertainty_var170.21Estimated variance inthe number of fragments generated by the object in sample 0 due tofragment assignment uncertainty.

5q0_count_dispersion_var4905.63Estimated variance inthe number of fragments generated by the object in sample 0 due tocross-replicate variability.

6q0_statusOKQuantification status for the object in sample 0. Can be one of OK(deconvolutionsuccessful), LOWDATA (too complex or shallowly sequenced), HIDATA

(too many fragments in locus), or FAIL, when an ill-conditioned

covariance matrix or other numerical exception preventsdeconvolution.

七. cufflinks使用中遇到的问题

1.

使用cuffdiff时候,在最新版本下,无重复的RNA-seq样作比较,结果中没有差异表达基因?

在v2.0.1及之后的版本中cuffdiff貌似不支持无重复的RNA-seq数据了。使用之前的版本即可。

八 Cuffquant

cuffquant是cuffquant能够对单个 BAM文件的基因转录本表达水平进行定量分析。生成的是CXB文件abundances.cxb,,可以作为cuffdiff的输入,这会加快cuffdiff的运行速度。也可以作为Cuffnorm的输入。

具体使用:Usage: cuffquant [options]*

它的参数:(和前面参数的含义是一样的)

-h/--help;-o/--output-dir

;-p/--num-threads ;-M/--mask-file

;-b/--frag-bias-correct

;-u/--multi-read-correct;--library-type;-m/--frag-len-mean

;-s/--frag-len-std-dev ;--max-mle-iterations

;--max-bundle-frags

;--no-effective-length-correction;--no-length-correction;-v/--verbose;-q/--quiet;--no-update-check;

九Cuffnorm

cuffnorm能够用 cuffquant

的输出文件作为输入文件,对基因和转录组,简单计算标准化过的表达水平。当你想要的是一系列可比较的基因、转录组、CDS 组和 TSS

组的表达值时,可是使用 cuffnorm。例如,当你仅仅想对单个基因的表达值做个热图或者点图时。

cuffnorm [options]* ...[sampleN.sam_replicate1.sam[,...,sample2_replicateM.sam]]

具体参数:它的参数和前面的类似,可以看前面的相关参数。

-h/--help;-o/--output-dir

;-L/--labels ;-p/--num-threads ;

--total-hits-norm(默认不激活);--compatible-hits-norm(默认激活);--library-type;--library-norm-method;--output-format;-v/--verbose;-q/--quiet;--no-update-check;

cuffnorm的输出文件是实验中的each gene, transcript, TSS group, andCDSgroup的标准化的表达水平。不做表达差异的分析。cuffnorm的输出文件默认是“simple-table”的文件。这些文件和cuffdiff输出的文件格式不同。若你想要cuffdiff格式的文件,你需要输入命令:--output-format cuffdiff

cuffnorm报道FPKM values andnormalized,estimates for the number of fragments that originate from eachgene, transcript, TSS group, and CDSgroup.这些结果已经做了标准化处理。对于某些下游软件需要原始文件,是不作为其输入的。

可以创建一个文件,例如sample_sheet.txt作为cuffdiff或cuffnorm的输入(存入sam文件的path)。文件格式如下:

sample_idgroup_label

C1_R1.samC1

C1_R2.samC1

C2_R1.samC2

C2_R2.samC2

输出结果文件如下:

FPKM tracking files:估测的基因的表达水平

Count tracking files:估测的基因的fragment count values

Read group tracking files:报道per-replicate expression and count

data.

对于每个genes, transcripts, TSS groups, and CDS

groups,cuffnorm会报道两种文件形式: *.fpkm_table files and *.count_table

files。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 193,812评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,626评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,144评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,052评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 60,925评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,035评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,461评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,150评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,413评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,501评论 2 307
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,277评论 1 325
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,159评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,528评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,868评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,143评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,407评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,615评论 2 335

推荐阅读更多精彩内容