wiki:SOPs/rna-seq-diff-expressions_TE

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

Background

  • Repetitive elements comprise a substantial portion of many eukaryotic genomes. In humans, for example, estimates of the repetitive fraction of the genome range from ​1/2 to more than ​2/3. Moreover, repetitive elements are known to play an important role in cellular function and disease, yet are typically excluded from the analysis of high-throughput sequencing experiments. This exclusion is due to the ambiguity that accompanies assigning multi-mapping reads (i.e. short reads that cannot be mapped uniquely to genomic loci). Strategies to address this ambiguity include assigning fractional reads to multiple matching loci or assigning such reads to ​subfamilies of repetitive elements. ​TEtranscripts leverages subfamily annotations to include transposable elements (TEs), along with the customary genes annotations, in the analysis of short read sequencing data. Transposable elements (e.g. LTRs, LINEs and SINEs) make up most of the repetitive DNA in the human genome, with the remainder being tandem repeats (e.g. satellites and microsatellites) that characterize heterochromatin and centromeres. The workflow below illustrates how to use ​TEtranscripts on the resources at the Whitehead Institute.

Step by step analysis

  • Mapping
    • Use ​STAR or another spliced mapper to map short reads to the genome of choice.
    • See our ​mapping SOP to search for details on running STAR.
    • When assessing transcription of TEs, it is essential to include multi-mapping reads. When using STAR, in particular, the winAnchorMultimapNmax and outFilterMultimapNmax flags are used to control multimapping, with the former setting a lower bound on how many loci must have a matching seed and the latter defining the upper bound on how many loci a read maps to in order to report it. A command for STAR mapping paired end reads using gzipped fastq input can look like: 
bsub STAR --genomeDir /path/to/STAR/index/for/organism --readFilesIn /path/to/reads_1.fastq.gz /path/to/reads_2.fastq.gz --outFileNamePrefix somePrefix --sjdbScore 2 --runThreadN 8 --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --winAnchorMultimapNmax 100 --outFilterMultimapNmax 100
  • Quantification of raw counts
    • TEtranscripts uses the BAM file(s) of aligned reads (from STAR in this example) as input.
    • TEtranscripts relies on separate gene annotation files (GTFs) for genes and TEs. A curated collection of TE GTFs can be found ​here and in the genome resources on the cluster at the Whitehead Institute (for human, mouse and fly).
    • Before assigning reads, it is important to know whether they are stranded (see the Quantification of raw counts section of our ​best practices page for details on how to determine this).
    • The best way to use the resources on the cluster to assign reads to genes and TEs is by running TEcount separately on each experiment (reverse stranded reads are shown in the example below, for forward stranded reads use --stranded forward and for unstranded reads use --stranded no (the default)): 
# Reverse stranded reads
bsub TEcount --sortByPos --format BAM --stranded reverse -b /path/to/alignment.bam --GTF /path/to/gene.gtf --TE /path/to/TE.gtf --mode multi --project projectName -i 100
  • The --sortByPos flag is necessary here because this was the sorting used in the STAR mapping, above.
  • The -I 100 (default) flag sets the maximum number of expectation maximization steps to take in computing maximum likelihood estimates of counts for repetitive elements.
  • Assessing differential expression for genes and TEs
    • After running TEcount on each sample in your experiment, the reported counts (i.e. a list of raw counts per gene/TE for each sample) can be combined into a counts matrix and analyzed following the steps outlined in the Statistics for differential expression, Identifying differentially expressed genes and Accounting for a batch effect in a differential expression model sections of our ​best practices page.
    • If the number of samples is not too large, the counting and analysis of differential expression can be carried out using a single execution of TEtranscripts instead of using TEcount (reverse stranded reads are shown in the example below, for forward stranded reads use --stranded forward and for unstranded reads use --stranded no (the default)): 
# Reverse stranded reads
bsub TEtranscripts --format BAM --stranded reverse -t /path/to/treat1.bam /path/to/treat2.bam -c /path/to/control1.bam /path/to/control2.bam --GTF /path/to/gene.gtf --TE /path/to/TE.gtf --mode multi --project treat_vs_control --minread 1 -i 100 --padj 0.05 --norm DESeq_default --sortByPos
  • The design for tests of differential expression above is a comparison between two biological contexts (e.g. treatment versus control, samples listed after the -t flag versus samples listed after the -c flag). If your experimental design is more complex, you should use TEcounts with a subsequent custom analysis of differential expression.
  • The output from TEtranscripts includes tests of differential expression carried out using ​DESeq2, in addition to a (sample x transcript) matrix of counts.
Note: See TracWiki for help on using the wiki.