转录组分析一条龙

1、转录组分析一条龙 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后的参数设定与转录组的大小有关，在内存

2、足够的情况下，设定大点能节约时间； 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，丢弃uni

3、quely 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 assemb

4、ly 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. tr

5、inity输出结果的统计 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可能有多

6、个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

7、/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_nam

8、esorted_to_uniq_count_stats.pl用于统计比对结果 $ $trinity_home/util/sam_namesorted_to_uniq_count_stats.pl bowtie_sorted.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 re

9、ported 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文件，从而可视化比对结果。

10、 4. 使用rsem进行表达量计算 首先，需要下载最新版本的rsem，安装并将程序加入到$path中。 $ wget /rsem/src/rsem-1.2.8.tar.gz $ tar zxf rsem-1.2.8.tar.gz $ cd rsem-1.2.8 $ make $ echo 使用$trinity_home/util/rsem_util/run_rsem_align_n_estimate.pl可以调用rsem，从而计算表达量。如果是链特异性测序，则加入ss_lib_type参数。 $trinity_home/util/rsem_ut

11、il/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中，则如上所

12、示。这两个参数分别代表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. 可以根据得到的结果，去除掉isop

13、ct低于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(bioclite(edger) bioclite(deseq) b

14、ioclite(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 $ trin

15、ity_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/differentialexpressio