Download strand-specific RNA-Seq data from Schizosaccharomyces pombe grown in four conditions45 (logarithmic growth, plateau phase, diauxic shift, and heat shock, each with 1M Illumina paired-end strand-specific RNA-Seq data, to a total of 4M paired-end reads) by visiting the URL reported below in a web browser, or directly from the command line using ‘wget’:

% wget http://sourceforge.net/projects/trinityrnaseq/files/misc/TrinityNatureProtocolTutorial.tgz/download

Name the file thus downloaded ‘TrinityNatureProtocolTutorial.tgz’, which should be 540 MB in size.

Unpack this file by:

tar –xvf TrinityNatureProtocolTutorial.tgz

Create a working folder and place the ‘TrinityNatureProtocolTutorial/’ directory contents there (as per the Materials section).

In order to facilitate downstream analyses, concatenate the RNA-Seq data across all samples into a single set of inputs to generate a single reference Trinity assembly. Combine all ‘left’ reads into a single file, and combine all ‘right’ reads into a single file by:

% cat 1M_READS_sample/*.left.fq > reads.ALL.left.fq

% cat 1M_READS_sample/*.right.fq > reads.ALL.right.fq

Assemble the reads into transcripts using Trinity by:

% $TRINITY_HOME/Trinity.pl --seqType fq --JM 10G \
    --left reads.ALL.left.fq --right reads.ALL.right.fq \
    --SS_lib_type RF --CPU 6

Use the script ‘$TRINITY_HOME/utilities/TrinityStats.pl’ to examine the statistic for the Trinity assemblies:

% $TRINITY_HOME/util/TrinityStats.pl trinity_out_dir/Trinity.fasta

Examine the breadth of genetic composition and transcript contiguity by leveraging a reference data set. The annotated reference transcriptome of Schizosaccharomyces pombe is included as file ‘S_pombe_refTrans.fasta’. Use megablast and our included analysis script to analyze its representation by the Trinity assembly as described below (Steps 9-13).

Prepare the reference transcriptome fasta file as a BLAST database:

% makeblastdb -in S_pombe_refTrans.fasta -dbtype nucl

Run megablast to align the known transcripts to the Trinity assembly:

% blastn -query trinity_out_dir/Trinity.fasta \
    -db S_pombe_refTrans.fasta \
    -out Trinity_vs_S_pombe_refTrans.blastn \
    -evalue 1e-20 -dust no -task megablast -num_threads 2 \
    -max_target_seqs 1 -outfmt 6

Once megablast has completed, run the script below to examine the length coverage of top database hits:

% $TRINITY_HOME/util/analyze_blastPlus_topHit_coverage.pl  \
          Trinity_vs_S_pombe_genes.blastn \
          trinity_out_dir/Trinity.fasta \
          S_pombe_refTrans.fasta

Examine the number of input RNA-Seq reads that are well represented by the transcriptome assembly. Trinity provides a script ‘alignReads.pl’ that executes Bowtie to align the left and right fragment reads separately to the Trinity contigs and then groups the reads together into pairs while retaining those single read alignments that are not found to be properly paired with their mate. Run ‘alignReads.pl’ by:

% $TRINITY_HOME/util/alignReads.pl --seqType fq \
                                   --left reads.ALL.left.fq --right reads.ALL.right.fq \
                                   --SS_lib_type RF --retain_intermediate_files \
                                   --aligner bowtie \
                                --target trinity_out_dir/Trinity.fasta -- -p 4

When ‘alignReads.pl’ is run using strand-specific data, as indicated above with the ‘--SS_lib_type RF’ parameter setting, it will separate the alignments into those that align to the sense strand (‘+’) from the antisense strand (‘-’). All output files including coordinate-sorted and read-name-sorted SAM files should exist in a ‘bowtie_out/’ directory. Count the number of reads aligning (at least once) to the sense strand of transcripts by running the utility below on the sense-strand read name-sorted alignment file as shown:

% $TRINITY_HOME/util/SAM_nameSorted_to_uniq_count_stats.pl   \
                           bowtie_out/bowtie_out.nameSorted.sam.+.sam

Obtain transcript abundance estimates by running RSEM separately for each sample, as shown below: The PERL script ‘run_RSEM_align_n_estimate.pl’ simply provides an interface to the RSEM software, translating the familiar Trinity command-line parameters to their RSEM equivalents and then executing the RSEM software. Each relevant step below (15–18) generates files ‘$prefix.isoforms.results’ and ‘$prefix.genes.results’, containing the abundance estimations for Trinity transcripts (Table 1) and components (Table 2), respectively. The $prefix in the filename is set based on the ‘--prefix’ setting in the commands below, which is unique to each sample. Please note that, with regard to the parameters reported in Table 1 and Table 2, Transcript_id is the Trinity transcript identifier; Gene_id is the Trinity component to which the reconstructed transcript was derived; Length is the length of the reconstructed transcript; Effective length is the mean number of 5' start positions from which an RNA-Seq fragment could have been derived from this transcript, given the distribution of fragment lengths inferred by RSEM — the value is equal to (transcript_length - mean_fragment_length + 1); Expected count is the number of expected RNA-Seq fragments assigned to the transcript given maximum likelihood transcript abundance estimates; TPM is the number of transcripts per million; FPKM is the number of RNA-Seq fragments per kilobase of transcript effective length per million fragments mapped to all transcripts; and IsoPct is the percent of expression for a given transcript compared to all expression from that Trinity component. Note also that the gene’s length and effective length are defined as the IsoPct weighted sum of transcript lengths and effective lengths. The gene’s expected counts, TPM, and FPKM are defined as the sum of its transcripts' expected counts, TPM and FPKM.

RSEM for log phase

% $TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl \
   --transcripts trinity_out_dir/Trinity.fasta \
   --left Sp.log.1M.left.fq \
   --right Sp.log.1M.right.fq \
   --seqType fq \
   --SS_lib_type RF \
   --prefix LOG

RSEM for diauxic shift

% $TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl \
  --transcripts trinity_out_dir/Trinity.fasta \
  --left Sp.ds.1M.left.fq \
  --right Sp.ds.1M.right.fq \
  --seqType fq \
  --SS_lib_type RF \
  --prefix DS

RSEM for heat shock

% $TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl \
  --transcripts trinity_out_dir/Trinity.fasta \
  --left Sp.hs.1M.left.fq \
  --right Sp.hs.1M.right.fq \
  --seqType fq \
  --SS_lib_type RF \
  --prefix HS

RSEM for plateau phase

% $TRINITY_HOME/util/RSEM_util/run_RSEM_align_n_estimate.pl \
   --transcripts trinity_out_dir/Trinity.fasta \
   --left Sp.plat.1M.left.fq \
   --right Sp.plat.1M.right.fq \
   --seqType fq \
   --SS_lib_type RF \
   --prefix PLAT

You will see that each of the RSEM ‘*.isoforms.results’ files has a number of columns, but we only need the one called ‘expected_count’. Create a matrix containing the counts of RNA-seq fragments per feature in a simple tab-delimited text file using the expected fragment count data produced by RSEM.

% $TRINITY_HOME/util/RSEM_util/merge_RSEM_frag_counts_single_table.pl\
    LOG.isoforms.results DS.isoforms.results HS.isoforms.results \
    PLAT.isoforms.results > Sp_isoforms.counts.matrix

Use edgeR to identify differentially expressed transcripts for each pair of samples. The following script automates many of the tasks of running edgeR or DESeq; in this procedure, we only leverage edgeR. Use the matrix created in Step 19 as input.

% $TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl \
  --matrix Sp_isoforms.counts.matrix \
  --method edgeR \
  --output edgeR_dir

To perform TMM normalization and to generate a matrix of expression values measured in FPKM, first extract the transcript length values from any one of RSEM’s *.isoform.results files:

% cut -f1,3,4 DS.isoforms.results > Trinity.trans_lengths.txt

Now, perform the TMM normalization:

% $TRINITY_HOME/Analysis/DifferentialExpression/run_TMM_normalization_write_FPKM_matrix.pl \
  --matrix Sp_isoforms.counts.matrix \
  --lengths Trinity.trans_lengths.txt

To study expression patterns of transcripts or genes across samples, it is often useful to restrict analysis to those transcripts that are significantly differentially expressed in at least one pairwise sample comparison. Given a set of differentially expressed transcripts, extract their normalized expression values and perform hierarchical clustering to group together transcripts with similar expression patterns across samples, and to group together those samples that have similar expression profiles according to transcripts. For example, enter the ‘edgeR_dir/’ output directory and extract those transcripts that are at least fourfold differentially expressed with false discovery–corrected statistical significance of at most 0.001 by using the following commands:

% cd edgeR_dir/
% $TRINITY_HOME/Analysis/DifferentialExpression/analyze_diff_expr.pl \
  --matrix ./Sp_isoforms.counts.matrix.TMM_normalized.FPKM \
  -C 2 -P 0.001

Determine the number of such differentially expressed transcripts by counting the number of lines in the file by using the command:

% wc -l diffExpr.P0.001_C2.matrix

Extract clusters of transcripts with common expression profiles from the earlier generated hierarchical clusters by running the script below, which uses R to cut the tree representing the hierarchically clustered transcripts based on specified criteria, such as to generate a specific number of clusters or by cutting the tree at a certain height. For example, run the following to partition transcripts by cutting the tree at 20% of its height:

% $TRINITY_HOME/Analysis/DifferentialExpression/define_clusters_by_cutting_tree.pl \
  --Ptree 20 -R diffExpr.P0.001_C2.matrix.R.all.Rdata

(Optional) Run the script in Step 25 several times with different values of ‘--Ptree’ in order to increase or decrease the number of clusters generated. Automating the required sections of the PROCEDURE ● TIMING 2–3 h

(Optional) If you are interested in executing the sections of the PROCEDURE without manually typing in each command, run the script below, which executes the required sections of the PROCEDURE. These sections include concatenating all samples’ reads into a single input data set, assembling the reads using Trinity, performing abundance estimation separately for each sample and running edgeR to identify differentially expressed transcripts. Run the automated procedure with the following command, including the -I (optional) parameter for an interactive experience, in which the system will pause and wait for a user response before proceeding to the next step.

% $TRINITY_HOME/util/run_Trinity_edgeR_pipeline.pl --samples_file samples_n_reads_described.txt -I