Using RNA-Seq to quantify gene levels and assay for differential expression for transposable elements

Background

• For each sample, map reads to genome using splice-aware mapper.
• Count reads mapped to each gene (or other set of features).
• Use gene counts to identify differentially expressed genes.

General suggestions

• Preliminary issues
  • Statistics for all methods require a matrix of counts (positive integer values) for each gene for each sample.
  • Create a tab-delimited matrix of integer counts, with column labels for each sample.
  • Genes with no counts in any sample should generally be removed to permit higher statistical power to identify differential expression.
  • According to ​Bullard et al., 2010, differential expression analysis is influenced more by the normalization method than by the choice of differential expression statistic.
  • Note that without replication, one cannot make very strong conclusions. High-throughput sequencing, just like every other technology, needs biological replication.
    • One can conclude that certain genes in sample A have a different RNA abundance than in sample B, but the results cannot be generalized.
    • Example, using an extremely precise balance: If Dick weighs more than Sally, we cannot conclude that males weigh more than females because we know nothing about the variability of weights among males and among females. Even if we weighed several individuals together, we'd still be missing information about within-group variability.
  • Sample commands to get raw counts from an alignment file:
    • coverageBed -split -abam accepted_hits.bam -b transcripts.gtf > transcript.coverage.bed (See the ​bedTools coverage page for details)
    • htseq-count -m intersection-strict --stranded=no accepted_hits.sam transcripts.gff > transcript.coverage.txt (See the ​htseq-count page for details)
    • In our view, htseq-count is better at handling reads that map to a genome region with overlapping genes.

Step by step analysis

• Mapping
  • Use ​STAR or another spliced mapper to map short reads to the genome of choice.
  • See our ​mapping SOP for more details.

• Quantification of raw counts

• Is your sequencing library stranded or unstranded? This information is needed to help these tools accurately count features. Strandedness of some library prep methods:
  • TruSeq Stranded mRNA Kits ("TruSeqStrandedPolyA") reads are reverse stranded (stranded in the reverse direction relative to the transcript orientation).
  • SMART-Seq v4 Ultra Low Input RNA Kit ("SMARTerUltra-lowPOLYA-V4") reads are unstranded.
  • KAPA RNA HyperPrep Kits ("KAPAHyperPrepmRNA") reads are reverse stranded.
  • The Whitehead Genome Core has some more ​Library Prep Descriptions.
• See SAMBAMqc (and/or look at mapped reads in a genome browser) to determine or verify strandedness

# single-end reads (unstranded)
featureCounts -a gene_anotations.gtf -o MySample.featureCounts.txt MySample.bam
# single-end reads (forward stranded)
featureCounts -s 1 -a gene_anotations.gtf -o MySample.featureCounts.txt MySample.bam
# single-end reads (reverse stranded)
featureCounts -s 2 -a gene_anotations.gtf -o MySample.featureCounts.txt MySample.bam


# paired-end reads (unstranded)
featureCounts -p -a gene_anotations.gtf -o MySample.featureCounts.txt MySample.bam
# paired-end reads (forward stranded)
featureCounts -p -s 1 -a gene_annotations.gtf -o MySamples.featureCounts.txt *sortedByName.bam
# paired-end reads (reverse stranded)
featureCounts -p -s 2 -a gene_annotations.gtf -o MySamples.featureCounts.txt *sortedByName.bam

• Other
  • Review articles:
    • ​Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. - Rapaport F, Khanin R, Liang Y, Pirun M, Krek A, Zumbo P, Mason CE, Socci ND, Betel D. Genome Biol. 2013 Sep 10;14(9):R95.
    • ​A survey of statistical software for analyzing RNA-seq data - Gao D, Kim J, Kim H, Phang TL, Selby H, Tan AC, Tong T. Hum Genomics. 2010 Oct;5(1):56-60.
    • ​From RNA-seq reads to differential expression results - Oshlack A, Robinson MD, Young MD. Genome Biol. 2010;11(12):220. Epub 2010 Dec 22.
  • For more practical information, see the third session of ​An introduction to R and Bioconductor: A BaRC Short Course and the ​BaRC Hot Topic (under "Short Read Sequencing", see "Practical RNA-Seq analysis")