转录组分析一条龙

Trinity 进行转录组分析的一条龙服务

1. Trinity进行转录组组装

Trinity 进行转录组组装的典型命令如下:

$ /opt/biosoft/trinityrnaseq_r20131110/Trinity.pl --seqType fq --JM 50G\ --left sample1_1.clean.fastq sample2_1.clean.fastq\

--right sample1_2.clean.fastq sample2_2.clean.fastq\

--jaccard_clip --CPU 6 --SS_lib_type FR

–JM 后的参数设定与转录组的大小有关，在内存足够的情况下，设定大点能节约时间；

–left 和 –right 后可以接多个样平的数据，并用空格隔开，值得注意的是，left reads name 以/1结尾，rigth reads name以/2结尾；

–jaccard_clip 适合于基因稠密的真菌物种；

–SS_lib_type 适合于链特异性测序

大数据量(>300M pairs)的RNA-seq 数据，最好使用

TRINITY_RNASEQ_ROOT/util/normalize_by_kmer_coverage.pl对reads 进行处理后再使用trinity 进行组装，以降低内存消耗和大量时间。

也可以设置–min_kmer_cov 2，丢弃uniquely occurring kmer, 从而降低内存消耗。 参考文献：

1. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev

A. Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol. 2011 May 15;29(7):644-52. doi: 10.1038/nbt.1883. PubMed PMID: 21572440.

2. Borodina T, Adjaye J, Sultan M. A strand-specific library preparation protocol for RNA sequencing. Methods Enzymol. 2011;500:79-98. PubMed PMID: 21943893.

2. Trinity输出结果的统计

Trinity 默认的输出结果为：trinity_out_dir/Trinity.fasta。

该fasta 格式文件中序列名例如：

>comp6749_c0_seq1 len=328 path=[471:0-83 388:84-208 679:209-327] >comp6749_c0_seq2 len=328 path=[304:0-83 388:84-208 679:209-327] >comp6749_c0_seq3 len=245 path=[901:0-125 679:126-244]

可以看到，trinity 生成的结果为components, 而一个components 可能有多个seq 。这相当于一个gene 能有多个transcripts 。

可以使用trinity 自带的程序TrinityStats.pl 对components 和transcripts 的数目，大小和N50等进行统计。

$ $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta Total trinity transcripts: 40138

Total trinity components:

Percent GC: 61.31 31067

3. 将reads 比对到转录组, 并进行可视化

TRINITY_RNASEQ_ROOT/util/alignReads.pl能调用bowtie 将reads map到转录组，并可以设置链特异性参数。

$ TRINITY_RNASEQ_ROOT/util/alignReads.pl --left left.fq --right right.fq --seqType fq\

--target Trinity.fasta --aligner bowtie --retain_intermediate_files 结果中生成coordSorted 和nameSorted 的sam 和bam 文件。如果设置了链特异性参数，则额外生成+链和-链的比对结果文件。

TRINITY_RNASEQ_ROOT/util/SAM_nameSorted_to_uniq_count_stats.pl用于统计比对结果

$ $TRINITY_HOME/util/SAM_nameSorted_to_uniq_count_stats.pl

bowtie_out.nameSorted.sam.+.sam

#read_type count pct

proper_pairs 21194964 93.22 both read pairs align to a single contig and point toward each other.

left_only 836213 3.68 only the left (/1) read is reported in an alignment

right_only 687576 3.02 only the right (/2) read is reported in an alignment

improper_pairs 16640 0.07 both left and right reads align, but to separate contigs, or to a single contig in the wrong expected relative orientations.

可以将Trinity.fasta 导入到IGV 中作为genome ，上载bam 文件，从而可视化比对结果。

4. 使用RSEM 进行表达量计算

首先，需要下载最新版本的RSEM ，安装并将程序加入到$PATH中。

$ wget http://deweylab.biostat.wisc.edu/rsem/src/rsem-1.2.8.tar.gz $ tar zxf rsem-1.2.8.tar.gz

$ cd rsem-1.2.8

$ make

$ echo "PATH=$PWD:\$PATH" >> ~/.bashrc

使用$TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl可以调用RSEM ，从而计算表达量。如果是链特异性测序，则加入–SS_lib_type参数。 $TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl

--transcripts Trinity.fasta \

--seqType fq --left left.reads.fq --right right.reads.fq --SS_lib_type FR \

--prefix RSEM --thread_count 4 -- --bowtie-phred64-quals --no-bam-output

将rsem-calculate-expression 程序的参数–bowtie-phred64-quals 和

–no-bam-output 加入到run_RSEM_align_n_estimate.pl中，则如上所示。这两个参数分别代表fastq 的质量格式是phred64，不输出bam 文件(节约大量时间) 。 若运行出现问题，点击：RSEM 的README 文件。

结果生成两个abundance estimation information文件:

RSEM.isoforms.results : EM read counts per Trinity transcript

RSEM.genes.results : EM read counts on a per-Trinity-component (aka… gene) basis, ‘gene’ used loosely here.

可以根据得到的结果，去除掉IsoPct 低于1%的transcripts 。可以依据

RSEM.isoforms.results 使用

TRINITY_RNASEQ_ROOT/util/filter_fasta_by_rsem_values.pl过滤掉trinity 组装结果中的lowly supported transcripts。

但不推荐过滤掉这些序列。

5. 鉴定差异表达transcripts

Trinity 可以使用Bioconductor package中的edgeR 或DESeq 来鉴定差异表达trancripts 。因此，需要安装R 和相关的一些包。

source("http://bioconductor.org/biocLite.R")

biocLite('edgeR')

biocLite('DESeq')

biocLite('ctc')

biocLite('Biobase')

install.packages('gplots’)

install.packages(‘ape’)

5.1 使用上一节中的RSEM 来分别对每个样品的每个生物学重复进行表达量计算

5.2 将每个样的RSEM 的结果进行合并

$ $TRINITY_HOME/util/RSEM_util/merge_RSEM_frag_counts_single_table.pl \

sampleA.RSEM.isoform.results sampleB.RSEM.isoform.results ... \ > transcripts.counts.matrix

$ TRINITY_HOME/util/RSEM_util/merge_RSEM_frag_counts_single_table.pl \ sampleA.RSEM.gene.results sampleB.RSEM.gene.results ... \ > genes.counts.matrix

然后修改生成的两个matrix 文件的column headers（代表着样品和重复的名字），有利于下游的分析。如果要分析transcripts 水平的差异表达，则使用transcripts.counts.matrix 文件；若要分析gene 水平的差异表达，则使用genes.counts.matrix 。

5.3 无生物学重复进行差异表达分析

$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl用于调用edgeR 或DESeq 进行差异表达基因分析。直接输入该命令查看其用法。 Trinty 推荐使用edgeR 进行差异表达分析。

$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl \ --matrix counts.matrix --method edgeR

注意输入的matrix 是counts 的数据，而不要是FPKM 的数据。

5.4 有生物学重复进行差异表达分析

首先，要建立文件samples_described.txt，内容为：

conditionA condA-rep1

conditionA condA-rep2

conditionB condB-rep1

conditionB condB-rep2

conditionC condC-rep1

conditionC condC-rep2

condA-rep1, condA-rep2, condB-rep1… 等对应着counts.matrix 文件中的column names 。

命令如下：

$TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl \ --matrix SP2.rnaseq.counts.matrix --method edgeR \

--samples_file samples_described.txt

结果文件中 logFC 是 log2 Fold Change; logCPM 是 log2-counts-per-million 。 值得注意的是：程序默认去除counts 数都少于10的transcripts 或genes ，不对其进行差异分析。所以有差异分析的genes 或transcripts 数目低于原始的数目。

5.5 提取差异表达基因，对其进行聚类分析

5. 5. 1 表达量的 N O R M A L I Z E D

使用TMM 方法将counts 转换为FPKM 。

首先从1个样平的RSEM 结果中提取长度数据：

$ cut -f 1,3,4 sampleA.RSEM.isoforms.results > feature_lengths.txt 然后使用TMM 方法将counts 数据转换为FPKM 数据:

$ $TRINITY_HOME/Analysis/DifferentialExpression/run_TMM_normalization_write_FPKM_matrix.pl \

--matrix counts.matrix --lengths feature_lengths.txt

5. 5. 2 提取差异表达转录子

注意的是，这一步要在edgeR 的结果文件中运行程序：

$ $TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl \ --matrix matrix.TMM_normalized.FPKM -P 0.001 -C 2

默认下选择FDR 值低于0.001，log2fold-change 的绝对值>=2为差异表达基因。 程序输出差异表达基因FPKM 、log2FC 、FDR 等值 和 聚类图 Heat Map.

5. 5. 3 根据聚类图提取子类

根据聚类结果，可以自动或手动确定子类。

自动确定子类：

$ $TRINITY_HOME/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl \

--Ptree 20 -R file.all.RData

上例中从数的20%处来自动划分子类。

手动确定子类：

$ R

> load("all.RData") # check for your corresponding .RData file name to use here, replace all.RData accordingly

>

source("$TRINITY_HOME/Analysis/DifferentialExpression/R/manually_define_clusters.R")

> manually_define_clusters(hc_genes, centered_data)

然后左键点击选择子类，右键结束选择

6. 提取蛋白编码区

使用transdecoder 从trinity 的转录子中提取coding region。最新版的transdecoder 貌似有点问题。

$ $TRINITY_HOME/trinity-plugins/transdecoder/transcripts_to_best_scoring_ORFs.pl \

-t transcripts.fasta -m 100

默认下允许的最小的protein 长度为100.

提取出了coding region，得出对应的protein 序列，有利于于下一步的功能注释。