RNA-seq Data Analysis

RNAseq: Next Generation Sequencing

• Platforms:
  • Illumina's Genome Analyzer
  • Applied Biosystems' SOLiD
  • Roche's 454 Life Sciences
• Terminology
  • Sequencing Depth or Coverage: Total number of reads mapped to the genome/transcriptome, also known as library size.
  • Transcript/gene length: Number of bases in a gene.
  • Read counts: Number of reads mapping to that gene/transcript (expression measurement). The reads are typically 30 ~ 400 bp, depending on the DNA-sequencing technology used.
• Illumina's sequencing technology
  • One flow cell: 8 lanes, one of which is often used for the control sample.

  • Multiplexing (see Conceptual overview)

    • Pooling: Sequencing multiple samples on a single unit (an Illumina's lane or flow cell)
    • Barcoding: Individual barcode sequences are added to each sample so they can be distinguished and sorted during data analysis.
      • Barcodes (barcode sequences) are used to separate and de-multiplex reads from each sample.
      • Barcodes in a single unit: 12 different samples can be indexed with unique subsequences and loaded onto each lane. In total, 96 samples can be sequenced per run.
    • Highlights:
      1. Fast, High-Throughput Strategy: Large sample numbers can be simultaneously sequenced during a single experiment
      2. Cost-Effective Method: Multisample pooling improves productivity by reducing time and reagent use
      3. High-Quality Data: Accurate maintenance of read length of unknown sequences (why?)
      4. Balanced Blocked Designs: A feasibility to construct balanced blocked designs for the purpose of testing differential expression.
  • Sequencing modes
    • Single-end Read: One read sequenced from one end of each cDNA insert
    • Paired-end Read: two reads sequenced from each cDNA sample insert (one from each end)
      • The costs of paired end sequencing are higher than single end sequencing
  • Quantitative standards (spike-ins, see ENCODE Consortium 2011)
    • It is highly desirable to include a ladder of RNA spike-ins to calibrate quantification, sensitivity, coverage and linearity.
  • Questions about the sequencing
    • Is the cost calculated by the number of lanes or flow cells used?
    • About $2,000 per lane?
    • How many samples can be in one lane?
    • How many reads can be gotten in one lane?
    • How many reads per sample can be gotten in one lane or one flow cell?
• RNA Sequencing Pipeline

RNAseq: Experimental Design

• Sequencing variations
  • Different genes have different variances and are potentially subject to different errors and biases.
  • Sources of variation affecting only a minority of genes should be integrated into the design as well (PCR-based GC bias).
  • Technical variability (experimental errors and biases): Repeated measurements of the same biological sample in multiple lanes or flow cells. For example, the same biological sample is in different lanes, which provides information about the variability of lanes. Two main sources of variation that may contribute to confounding effects:
    1. Batch effects: errors that occur after random fragmentation of the RNA until it is input to the flow cell (PCR, reverse transcription).
    2. Lane effects: errors that occur from the flow cell until obtaining the data from the sequencing machine (bad sequencing cycles, base-calling)
  • Biological variability (see ENCODE Consortium 2011)
    1. A biological replicate is defined as an independent growth of cells/tissue and subsequent analysis.
    2. Experiments should be performed with two or more biological replicates, unless there is a compelling reason why this is impractical or wasteful (e.g. overlapping time points with high temporal resolution).
    3. A typical R2 (Pearson) correlation of gene expression (RPKM) between two biological replicates, for RNAs that are detected in both samples using RPKM or read counts, should be between 0.92 to 0.98. Experiments with biological correlations that fall below 0.9 should be either be repeated or explained.
• Sequencing depth (see ENCODE Consortium 2011)
  • The amount of sequencing needed for a given sample is determined by the goals of the experiment and the nature of the RNA sample.

  • Experiments whose purpose is to evaluate the similarity between the transcriptional profiles of two polyA+ samples may require only modest depths of sequencing (e.g. 30M pair-end reads of length > 30NT, of which 20-25M are mappable to the genome or known transcriptome)

  • Experiments whose purpose is discovery of novel transcribed elements and strong quantification of known transcript isoforms requires more extensive sequencing. The ability to detect reliably low copy number transcripts/isoforms depends upon the depth of sequencing and on a sufficiently complex library.
  • For experiments from a typical mammalian tissue or in which sensitivity of detection is important, a minimum depth of 100-200 M 2 x 76 bp or longer reads is currently recommended.
• Purposes of the experimental design
  • avoid or eliminate the technical variation (possible confounding factors)
  • estimate the biological variation
• Sequencing design
  • Sampling: subject sampling, RNA sampling, and fragment sampling
  • Randomization: assigning individuals at random to groups (reduce the sample variability or variation)
  • Replication: biological replicates allow for estimation of within-treatment group (biological) variability, which is needed for making inferences between treatment groups.
  • Blocking: experimental units are grouped into homogeneous clusters (blocks)
• Balanced block designs
  • Barcoding: DNA fragments can be labeled or barcoded with sample specific sequences that allow multiple samples to be included in the same sequencing.
  • Pooling: All the samples of RNA are pooled into the same batch and then sequenced in one lane of a flow cell.
  • Any batch and lane effects are the same for all the samples.
• Balanced incomplete block designs (BIBD)
  • Technical constraints and the scientific hypotheses:
    1. the number of treatments (I)
    2. the number of biological replicates per treatment (J)
    3. the number of unique barcodes (s) that can be included in a single lane (block)
    4. the number of lanes available for sequencing (L)

  • When s < I, i.e., the number of unique bar codes in one lane is less than the number of treatments or samples,

    • a complete block design is not be possible.
    • In these cases, a BIBD is suggested.
  • Balanced incomplete block:
    • Incomplete: cannot fit all treatments (samples) in each block
    • Balanced: each pair of treatments occur together in the same number of times, and then the variance of the difference between two treatments is constant.
  • BIBD:
    • The total number of possible technical replicates per biological replicate is T = sL/JI.
    • The number of times each pair of treatments occurs together is k = J(s-1)/(I-1), an integer.
    • Extensive list of BIBD can be found in Fisher and Yates (1963) and Cochran and Cox (1957)

RNAseq: Statistical Analysis

RNAseq data analysis using edgeR package

• Filtering very low expressed genes using a CPM cutoff
  • Very low expressed genes: Those genes should not be differentially expressed but could increase the false positive error in differential testings. Those genes could cause that the equivalently expressed (EE) genes are not uniformly distributed.
  • edgeR User Guide (page 11):
    1. Genes with very low counts across all libraries provide little evidence for differential expression. In the biological point of view, a gene must be expressed at some minimal level before it is likely to be translated into a protein or to be biologically important.
    2. The pronounced discreteness of these counts interferes with some of the statistical approximations that are used later in the pipeline. These genes should be filtered out prior to further analysis.
    3. As a rule of thumb, genes are kept if they are expressed in at least one condition. Usually a gene is required to have a count of 5-10 in a library to be considered expressed in that library. Users should also filter with count-per-million (CPM) rather than filtering on the counts directly, as the latter does not account for differences in library sizes between samples.
    4. It is also recommended to recalculate the library sizes of the DGEList object after the filtering though the difference is usually negligible.
  • CPM cutoff = min.counts*1e6/min(colSums(counts))
    • min.counts = 5 ~ 10
    • quite rough as the exact cutoff isn't important.
    • choosing a appropriate minimum counts by comparing the similarity of sample distributions from boxplots, density plots, and/or wilcox.test
• Normalization
• Differential expression testing
  • Design matrix, comparison or contract
  • Dispersion estimation
  • Statistical testings (model fitting)
    1. ET (exact tests): exactTest
    2. LRT: glmFit, glmLRT
    3. QLF: glmQLFit, glmQLFTest
• Interpreting the testing results
  • QQ-plot of pvalues: qqpvalue,
    • significant if there are some QQ-plot points that are above the confidence interval
    • conservative if the QQ-plot curve is under a confidence interval
    • liberal if the QQ-plot curve is far above a confidence interval
    • confounding or unknown effects if the QQ-plot curve straightly deviates the diagonal line
    • correlation: lots of correlated genes should be observed.
  • QQ-plots of z-values or chi-squared values: qqstats
  • adjusted pvalue (FDR): topTags
  • MAplot, volcano plot
  • pvalue or FDR scores

References

Auer and Doerge (2010) Statistical Design and Analysis of RNA Sequencing Data, Genetics
Illumina's technical note: Estimating Sequencing Coverage
Some Issues of Statistical Design & Analysis in RNA-seq Experiment, Shaheena Bashir (2012)
Standards, Guidelines and Best Practices for RNA-Seq V1.0. The ENCODE Consortium (June 2011)