-
Notifications
You must be signed in to change notification settings - Fork 3
expression_quantification
To estimate the expression levels of the Trinity-reconstructed transcripts, we use the strategy supported by the RSEM software involving read alignment followed by expectation maximization to assign reads according to maximum likelihood estimates. In essence, we first align the original rna-seq reads back against the Trinity transcripts, then run RSEM to estimate the number of rna-seq fragments that map to each contig. Because the abundance of individual transcripts may significantly differ between samples, the reads from each sample (and each biological replicate) must be examined separately, obtaining sample-specific abundance values.
We include a script to facilitate running of RSEM on Trinity transcript assemblies. The script we execute below will run the Bowtie aligner to align reads to the Trinity transcripts, and RSEM will then evaluate those alignments to estimate expression values. Again, we need to run this separately for each sample and biological replicate (ie. each pair of fastq files).
Let's start with one of the wt37 fastq pairs like so:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left wt37_rep1_1.fastq \
--right wt37_rep1_2.fastq \
--transcripts Trinity.fasta \
--output_prefix wt37_rep1 \
--est_method RSEM --aln_method bowtie \
--trinity_mode \
--output_dir wt37_rep1_rsem \
--prep_reference
The outputs generated from running the command above will exist in the wt37_rep1_rsem/ directory, as we indicate with the --output_dir parameter above. Since this is the first time we're running RSEM, we need to prepare the Trinity.fasta file for alignment with the bowtie software, and so we'll simply include the --prep_reference parameter to invoke that functionality.
The primary output generated by RSEM is the file containing the expression values for each of the transcripts. Examine this file like so:
% head wt37_rep1_rsem/wt37_rep1.isoforms.results | column -t
and you should see the top of a tab-delimited file:
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct
TRINITY_DN0_c0_g1_i1 TRINITY_DN0_c0_g1 2010 1881.38 225.00 72.01 65.65 100.00
TRINITY_DN0_c0_g2_i1 TRINITY_DN0_c0_g2 2011 1882.38 0.00 0.00 0.00 100.00
TRINITY_DN1000_c0_g1_i1 TRINITY_DN1000_c0_g1 1664 1535.38 208.00 81.57 74.37 100.00
TRINITY_DN1001_c0_g1_i1 TRINITY_DN1001_c0_g1 2948 2819.38 8.00 1.71 1.56 100.00
TRINITY_DN1002_c0_g1_i1 TRINITY_DN1002_c0_g1 351 222.54 1.00 2.71 2.47 100.00
TRINITY_DN1003_c0_g1_i1 TRINITY_DN1003_c0_g1 1913 1784.38 91.00 30.71 28.00 100.00
TRINITY_DN1004_c0_g1_i1 TRINITY_DN1004_c0_g1 2395 2266.38 167.00 44.37 40.45 100.00
TRINITY_DN1005_c0_g1_i1 TRINITY_DN1005_c0_g1 3076 2947.38 425.00 86.82 79.16 100.00
TRINITY_DN1006_c0_g1_i1 TRINITY_DN1006_c0_g1 1141 1012.38 253.00 150.47 137.19 100.00
The key columns in the above RSEM output are the transcript identifier, the 'expected_count' corresponding to the number of RNA-Seq fragments predicted to be derived from that transcript, and the 'TPM' or 'FPKM' columns, which provide normalized expression values for the expression of that transcript in the sample.
The FPKM expression measurement normalizes read counts according to the length of transcripts from which they are derived (as longer transcripts generate more reads at the same expression level), and normalized according to sequencing depth. The FPKM acronym stand for 'fragments per kilobase of cDNA per million fragments mapped'.
TPM 'transcripts per million' is generally the favored metric, as all TPM values should sum to 1 million, and TPM nicely reflects the relative molar concentration of that transcript in the sample. FPKM values, on the other hand, do not always sum to the same value, and do not have the similar property of inherrently representing a proportion within a sample, making comparisons between samples less straightforward. TPM values can be easily computed from FPKM values like so: TPMi = FPKMi / (sum all FPKM values) * 1 million.
Note, a similar RSEM outputted gene expression file exists for the 'gene' level expression data. For example:
% head wt37_rep1_rsem/wt37_rep1.genes.results | column -t
.
gene_id transcript_id(s) length effective_length expected_count TPM FPKM
TRINITY_DN0_c0_g1 TRINITY_DN0_c0_g1_i1 2010.00 1881.38 225.00 72.01 65.65
TRINITY_DN0_c0_g2 TRINITY_DN0_c0_g2_i1 2011.00 1882.38 0.00 0.00 0.00
TRINITY_DN1000_c0_g1 TRINITY_DN1000_c0_g1_i1 1664.00 1535.38 208.00 81.57 74.37
TRINITY_DN1001_c0_g1 TRINITY_DN1001_c0_g1_i1 2948.00 2819.38 8.00 1.71 1.56
TRINITY_DN1002_c0_g1 TRINITY_DN1002_c0_g1_i1 351.00 222.54 1.00 2.71 2.47
TRINITY_DN1003_c0_g1 TRINITY_DN1003_c0_g1_i1 1913.00 1784.38 91.00 30.71 28.00
TRINITY_DN1004_c0_g1 TRINITY_DN1004_c0_g1_i1 2395.00 2266.38 167.00 44.37 40.45
TRINITY_DN1005_c0_g1 TRINITY_DN1005_c0_g1_i1 3076.00 2947.38 425.00 86.82 79.16
Kallisto is a newly developed tool for doing rna-seq based quantitation of transcript expression. It is both very accurate and very fast. Try running kallisto on this same set of reads for the sake of comparison.
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left wt37_rep1_1.fastq \
--right wt37_rep1_2.fastq \
--transcripts Trinity.fasta \
--output_prefix wt37_rep1 \
--est_method kallisto \
--trinity_mode \
--output_dir wt37_rep1_kallisto \
--prep_reference
In the case of kallisto, you'll find the transcript expression results in the file 'wt37_rep1_kallisto/abundance.tsv'. View it like so:
% head wt37_rep1_kallisto/abundance.tsv | column -t
.
target_id length eff_length est_counts tpm
TRINITY_DN144_c0_g1_i1 4833 4703.42 138 16.266
TRINITY_DN144_c0_g2_i1 2228 2098.42 0.000103136 2.72479e-05
TRINITY_DN179_c0_g1_i1 1524 1394.42 227 90.2502
TRINITY_DN159_c0_g1_i1 659 529.534 7.75713 8.12123
TRINITY_DN159_c0_g2_i1 247 119.949 0.24287 1.12251
TRINITY_DN153_c0_g1_i1 2378 2248.42 16 3.9451
TRINITY_DN130_c0_g1_i1 215 89.2898 776 4818.09
TRINITY_DN130_c1_g1_i1 295 166.986 216 717.115
TRINITY_DN106_c0_g1_i1 4442 4312.42 390 50.137
and similarly, we have a file containing the 'gene' level expression estimates:
% head wt37_rep1_kallisto/abundance.tsv.genes | column -t
.
target_id length eff_length est_counts tpm
TRINITY_DN2774_c0_g1 2926.00 2796.42 31.00 6.15
TRINITY_DN5482_c0_g1 3064.00 2934.42 344.00 64.99
TRINITY_DN6803_c0_g1 1439.00 1309.42 1379.00 583.85
TRINITY_DN386_c0_g2 4279.00 4149.42 3.23 0.43
TRINITY_DN23_c0_g2 632.00 502.53 9.99 11.02
TRINITY_DN5348_c0_g1 2091.00 1961.42 264.00 74.62
TRINITY_DN5222_c0_g1 2416.00 2286.42 148.00 35.89
TRINITY_DN4680_c0_g1 1420.00 1290.42 167.00 71.75
TRINITY_DN2900_c0_g1 283.00 155.12 1.00 3.57
Running this on all the samples can be monotonous, and with many more samples, advanced users would generally write a short script to fully automate this process. We won't be scripting here, but instead just directly execute abundance estimation just as we did above but on the other five pairs of fastq files. We'll examine the outputs of each in turn, as a sanity check and also just for fun.
Abundance estimates for wt37_rep2:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left wt37_rep2_1.fastq \
--right wt37_rep2_2.fastq \
--transcripts Trinity.fasta \
--output_prefix wt37_rep2 \
--est_method kallisto \
--trinity_mode \
--output_dir wt37_rep2_kallisto
Abundance estimates for wt37_rep3:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left wt37_rep3_1.fastq \
--right wt37_rep3_2.fastq \
--transcripts Trinity.fasta \
--output_prefix wt37_rep3 \
--est_method kallisto \
--trinity_mode \
--output_dir wt37_rep3_kallisto
Now we're done with processing the wt37 biological replicates, and we'll proceed to now run abundance estimations for the GSNO samples.
Abundance estimates for GSNO_rep1:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left GSNO_rep1_1.fastq \
--right GSNO_rep1_2.fastq \
--transcripts Trinity.fasta \
--output_prefix GSNO_rep1 \
--est_method kallisto \
--trinity_mode \
--output_dir GSNO_rep1_kallisto
Abundance estimates for GSNO_rep2:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left GSNO_rep2_1.fastq \
--right GSNO_rep2_2.fastq \
--transcripts Trinity.fasta \
--output_prefix GSNO_rep2 \
--est_method kallisto \
--trinity_mode \
--output_dir GSNO_rep2_kallisto
Abundance estimates for GSNO_rep3:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left GSNO_rep3_1.fastq \
--right GSNO_rep3_2.fastq \
--transcripts Trinity.fasta \
--output_prefix GSNO_rep3 \
--est_method kallisto \
--trinity_mode \
--output_dir GSNO_rep3_kallisto
And now let's tackle abundance estimation for each of the ph8 sample replicates.
Abundance estimates for ph8_rep1:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left ph8_rep1_1.fastq \
--right ph8_rep1_2.fastq \
--transcripts Trinity.fasta \
--output_prefix ph8_rep1 \
--est_method kallisto \
--trinity_mode \
--output_dir ph8_rep1_kallisto
Abundance estimates for ph8_rep2:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left ph8_rep2_1.fastq \
--right ph8_rep2_2.fastq \
--transcripts Trinity.fasta \
--output_prefix ph8_rep2 \
--est_method kallisto \
--trinity_mode \
--output_dir ph8_rep2_kallisto
Abundance estimates for ph8_rep3:
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--left ph8_rep3_1.fastq \
--right ph8_rep3_2.fastq \
--transcripts Trinity.fasta \
--output_prefix ph8_rep3 \
--est_method kallisto \
--trinity_mode \
--output_dir ph8_rep3_kallisto
For each of these abundance estimation runs, you should find the corresponding output directories and abundance estimation tables within them.
Now, given the expression estimates for each of the transcripts in each of the samples, we're going to pull together expression values into matrices containing transcript IDs in the rows, and sample replicate names in the columns. We'll make two matrices, one containing the estimated counts, and another containing the TPM expression values that are cross-sample normalized using the TMM method. This is all done for you by the following script in Trinity, indicating the method we used for expression estimation and providing the list of individual sample abundance estimate files:
% $TRINITY_HOME/util/abundance_estimates_to_matrix.pl \
--est_method kallisto --out_prefix Trinity_trans \
--name_sample_by_basedir \
wt37_rep1_kallisto/abundance.tsv \
wt37_rep2_kallisto/abundance.tsv \
wt37_rep3_kallisto/abundance.tsv \
GSNO_rep1_kallisto/abundance.tsv \
GSNO_rep2_kallisto/abundance.tsv \
GSNO_rep3_kallisto/abundance.tsv \
ph8_rep1_kallisto/abundance.tsv \
ph8_rep2_kallisto/abundance.tsv \
ph8_rep3_kallisto/abundance.tsv
You should find a matrix file called 'Trinity_trans.counts.matrix', which contains the counts of RNA-Seq fragments mapped to each transcript. Examine the first few lines of the counts matrix:
% head -n20 Trinity_trans.counts.matrix | column -t
.
wt37_rep1_kallisto wt37_rep2_kallisto wt37_rep3_kallisto GSNO_rep1_kallisto GSNO_rep2_kallisto GSNO_rep3_kallisto ph8_rep1_kallisto ph8_rep2_kallisto ph8_rep3_kallisto
TRINITY_DN323_c0_g1_i1 846 782 792 1403 1397 1062 1296 1291 1382
TRINITY_DN2438_c0_g1_i1 418 364 353 13 10 11 72 74 58
TRINITY_DN4819_c0_g1_i1 136 128 165 58 64 52 84 94 66
TRINITY_DN1223_c0_g1_i1 7 4 6 6 9 2 1 1 3
TRINITY_DN1807_c0_g1_i1 632 677 579 164 156 135 117 69 82
TRINITY_DN3825_c0_g1_i1 270 242 291 105 123 81 178 178 163
TRINITY_DN669_c0_g4_i1 0.00139273 0 0 0.0122585 0 22 0.00611495 0.000680272 0.009266
TRINITY_DN3817_c0_g1_i1 104 130 157 344 361 250 1051 1065 1126
TRINITY_DN4696_c0_g1_i1 13 5 4 27 31 27 41 25 26
TRINITY_DN5720_c0_g1_i1 395 428 464 699 737 607 347 382 382
TRINITY_DN246_c0_g2_i1 22.0612 0 0 0 9.98855 0 0 0 0
TRINITY_DN2531_c0_g1_i1 0 1 0 13 17 11 2 4 3
TRINITY_DN7039_c0_g1_i1 6 4 7 6 4 6 6 6 3
TRINITY_DN1166_c0_g1_i1 30 37 23 8 9 17 13 16 12
TRINITY_DN5294_c0_g1_i1 96 70 125 243 244 186 627 537 560
TRINITY_DN2674_c0_g1_i1 3 1 6 4 3 3 3 0 0
TRINITY_DN707_c0_g1_i1 191 190 193 62 74 55 193 181 154
TRINITY_DN4048_c0_g1_i1 1418 1173 1292 8430 8173 7243 3845 4091 4051
TRINITY_DN513_c0_g1_i1 90 57 87 191 209 161 162 136 147
You'll see that the above matrix has integer values representing the number of RNA-Seq paired-end fragments that are estimated to have been derived from that corresponding transcript in each of the samples. Don't be surprised if you see some values that are not exact integers but rather fractional read counts. This happens if there are multiply-mapped reads (such as to common sequence regions of different isoforms), in which case the multiply-mapped reads are fractionally assigned to the corresponding transcripts according to their maximum likelihood.
The counts matrix will be used by edgeR (or other tools in Bioconductor, eg. DESeq2) for statistical analysis and identifying significantly differentially expressed transcripts.
Now take a look at the first few lines of the normalized expression matrix:
% head -n20 Trinity_trans.TMM.EXPR.matrix | column -t
.
wt37_rep1_kallisto wt37_rep2_kallisto wt37_rep3_kallisto GSNO_rep1_kallisto GSNO_rep2_kallisto GSNO_rep3_kallisto ph8_rep1_kallisto ph8_rep2_kallisto ph8_rep3_kallisto
TRINITY_DN323_c0_g1_i1 48.683 46.605 41.263 86.390 84.111 80.252 61.934 65.842 70.628
TRINITY_DN2438_c0_g1_i1 171.798 155.077 134.165 5.737 4.305 5.934 25.409 27.508 21.668
TRINITY_DN4819_c0_g1_i1 174.458 170.392 206.578 80.622 86.288 87.534 100.517 114.944 81.709
TRINITY_DN1223_c0_g1_i1 9.443 5.597 7.932 8.778 12.765 3.540 1.265 1.291 3.924
TRINITY_DN1807_c0_g1_i1 460.899 512.190 398.106 128.816 119.278 129.182 75.552 46.366 55.530
TRINITY_DN3825_c0_g1_i1 55.670 51.695 54.804 23.199 26.546 21.927 30.905 32.804 30.138
TRINITY_DN669_c0_g4_i1 0.003 0.000 0.000 0.030 0.000 63.483 0.013 0.002 0.021
TRINITY_DN3817_c0_g1_i1 64.563 83.701 91.291 229.767 234.888 203.669 571.770 605.334 644.361
TRINITY_DN4696_c0_g1_i1 8.421 3.359 2.431 18.823 21.049 22.951 23.334 14.851 15.554
TRINITY_DN5720_c0_g1_i1 379.156 426.502 426.032 724.456 742.406 764.476 301.634 342.615 345.989
TRINITY_DN246_c0_g2_i1 6.725 0.000 0.000 0.000 3.188 0.000 0.000 0.000 0.000
TRINITY_DN2531_c0_g1_i1 0.000 5.021 0.000 72.715 90.670 71.765 10.387 21.234 16.164
TRINITY_DN7039_c0_g1_i1 13.139 9.054 15.705 14.395 9.252 17.241 13.042 13.197 6.704
TRINITY_DN1166_c0_g1_i1 18.340 23.459 13.162 5.262 5.767 13.639 6.958 8.951 6.758
TRINITY_DN5294_c0_g1_i1 42.611 32.210 51.411 115.848 113.451 108.366 239.727 215.991 226.430
TRINITY_DN2674_c0_g1_i1 8.249 2.830 17.239 12.135 8.750 10.830 8.399 0.000 0.000
TRINITY_DN707_c0_g1_i1 55.410 57.119 51.400 19.290 22.477 20.947 47.524 47.161 40.285
TRINITY_DN4048_c0_g1_i1 1623.714 1394.333 1431.606 10443.192 9829.451 10882.371 4056.834 4427.397 4434.405
TRINITY_DN513_c0_g1_i1 17.091 11.214 15.076 38.860 41.542 40.142 25.866 23.063 25.006
These are the normalized expression values, which have been further cross-sample normalized using TMM normalization to adjust for any differences in sample composition. TMM normalization assumes that most transcripts are not differentially expressed, and linearly scales the expression values of samples to better enforce this property. TMM normalization is described in A scaling normalization method for differential expression analysis of RNA-Seq data, Robinson and Oshlack, Genome Biology 2010.
We use the TMM-normalized expression matrix when plotting expression values in heatmaps and other expression analyses.
In addition to the transcript count matrix we created earlier, let's now create a 'gene' level counts matrix like so:
% $TRINITY_HOME/util/abundance_estimates_to_matrix.pl \
--est_method kallisto --out_prefix Trinity_genes \
--name_sample_by_basedir \
wt37_rep1_kallisto/abundance.tsv.genes \
wt37_rep2_kallisto/abundance.tsv.genes \
wt37_rep3_kallisto/abundance.tsv.genes \
GSNO_rep1_kallisto/abundance.tsv.genes \
GSNO_rep2_kallisto/abundance.tsv.genes \
GSNO_rep3_kallisto/abundance.tsv.genes \
ph8_rep1_kallisto/abundance.tsv.genes \
ph8_rep2_kallisto/abundance.tsv.genes \
ph8_rep3_kallisto/abundance.tsv.genes
You should now find files 'Trinity_genes.counts.matrix' and 'Trinity_genes.TMM.EXPR.matrix' in your workspace area. We'll use these later when we perform differential expression analysis at both the transcript and at the gene levels separately.
% $TRINITY_HOME/util/align_and_estimate_abundance.pl --seqType fq \
--samples_file samples.txt \
--transcripts Trinity.fasta \
--est_method kallisto \
--trinity_mode \
--prep_reference
For those replicates that you want to view in IGV, you'll need to run bowtie on them separately. First build the bowtie index for your Trinity.fasta file:
% bowtie2-build Trinity.fasta Trinity.fasta
Next align the reads to the Trinity assembly and generate a readname-sorted bam file: (replicate.left.fastq and replicate.right.fastq should correspond to your replicate paired-end files of interest)
% bowtie2 --local --no-unal -x Trinity.fasta -q \
-1 replicate.left.fastq -2 replicate.right.fastq \
| samtools view -Sb - | samtools sort -o replicate.bowtie2.coordSorted.bam
Run the 'find' command to locate the output files containing the expression data:
% find rep_dir_prefix* -name "*abundance.tsv" | tee quant_files.txt
Now, build a matrix using that 'quant_files.txt' containing list of abundance files.
% $TRINITY_HOME/util/abundance_estimates_to_matrix.pl \
--est_method kallisto --name_sample_by_basedir \
--quant_files quant_files.txt \
--out_prefix isoforms
Note, you can do the same for the 'genes' files, use --out_prefix genes