RNA-Seq

Setup the Data and Workflow

1. Import raw fastq reads into Galaxy as a Collection
   • Get Data -> Upload File
     • If uploading data from another online location (e.g. Princeton's HTSEQ database), it is useful to use the "Rule Based Uploader" tab, and paste a table of sample names and URLs
     • If upload from your computer, use the "Collection" tab. Once you've uploaded the files, click the "Build" button to assign the appropriate sample name to each file and create a collection.
2. Import a copy of the RNA-Seq workflow
3. Run FastQC on fastq file(s) for each sample/replicate. Look for low quality, bias in sequence, etc.

• NGS: QC and manipulation -> Fastqc: Fastq QC
  • Use Per Base Sequence Content and Kmer Content to determine if adapter trimming is necessary

3 - Optional - Trim poor quality bases and/or adapter sequences from reads

• NGS: QC and manipulation -> Cutadapt
• Qualities are in Phred Scale, usually 20 or 30 are good thresholds

11 - Load the data into R

dds <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,

    directory = directory, design= ~ condition)

colData(dds)$condition <- factor(colData(dds)$condition,

    levels=c("untreated","treated"))

• Note: It is important to supply levels (otherwise the levels are chosen in alphabetical order) and to put the “control” or “untreated” level as the first element, so that the log2 fold changes and results will be most easily interpretable.

12 - Normalize samples

dds <- estimateSizeFactors(dds)

13 - Sample QC: Transform data to log space and visualize samples

rld <- rlogTransformation(dds, blind = TRUE)

print(plotPCA(rld, intgroup = c("condition")))

14 - Estimate biological variance and visualize

dds <- estimateDispersions(dds)

plotDispEsts(dds)

15 - Determine differential expression

dds <- nbinomWaldTest(dds)

16 - Export results to file

# Get results and write to output files

res <- results(dds)

res <- res[order(res$padj), ]

write.table(cbind('#ID'=rownames(as.data.frame(res)), as.data.frame(res)),

    file='deseq2_results.txt', quote=FALSE, sep='\t', row.names=FALSE) 

# Display the top few genes

head(res)

17 - Visualize with MA plot

plotMA(dds,ylim=c(-2,2),main="DESeq2")

Things to Consider

Mapping

Splice Junctions

Mapping reads across splice junctions is one of the main challenges in RNA-Seq. There are a few approaches:

1. Exon junction libraries - a reference sequence is constructed using annotated exon boundaries
   • Requires less coverage, but doesn't detect novel junctions. This approach can also lead to better reproducibility between samples. Suitable when the reference transcriptome is well known.
   • Using Tophat/Cufflinks with a reference GTF uses this approach.
2. De novo splice junction detection - the reads are used to determine splice junctions by mapping parts of reads independently to the reference
   • Requires greater coverage to ensure correct detection. Also slower due to increased computational demands.
   • This is the default method for TopHat and Cufflinks.
3. De novo transcriptome assembly - No reference is used, but the reads are completely assembled
   • Requires high coverage and high quality data. Most approaches require a large amount of memory.
   • Trinity is an example of software that takes this approach.

When using either of the de novo approaches, there are few additional considerations:

1. Results may be inconsistent between samples.
2. Since the search is done separately for each sample, a particular gene/junction may be found in one sample, but not in another simply due to differences in expression/coverage, etc. In order to compensate for this, it may be best to create a "master" gene model by combining the results of the novel gene search from multiple samples and then rerunning the alignments providing this master gene model as input and turning off the novel search. This obviously adds a bit of complexity and time to the analysis, but may be worthwhile. There are a couple of ways you could do this.
   1. Combining all of the data together for the initial "search phase"
   2. Combining the results of multiple searches (e.g. using cuffmerge)
3. Many "genes" will be unannotated.
4. The results you get for expression/etc. will often be for gene_#### which has no annotations other than location in the genome. While this obviously allows for a lot of discovery, interpretation can be difficult.

Unique vs. Multi-Mapped Reads

When mapping reads to the genome or transcriptome, many reads map best to one particular location. However, some map equally well to one or more places, usually due to repetitive regions in the genome (multiple copies of genes, etc.). Most mapping algorithm allow to control what happens with these "multi-mappers". By default, TopHat will map a read in up to 20 places in the genome. Reads that map to more locations are discarded. When calculating coverage, the count is divided by the number of places a read mapped to. For example, a read that maps to two places will contribute 0.5 read to the count for each position.

Cufflinks can also attempt to improve the estimation using the 'rescue' method described in Ali Mortazavi, Brian A Williams, Kenneth McCue, Lorian Schaeffer and Barbara Wold, Mapping and quantifying mammalian transcriptomes by RNA-Seq, Nature Methods, Volume 5, 621 - 628 (2008).

In a BAM file created by TopHat there is an auxiliary tag (NH) that specifies the number of hits each read has. For example, NH:i:2 says that there are two hits for that read. Other mapping programs may encode this information differently.

Samtools flagstat and FastQC read count - Both samtools flagstat and FastQC will report one "read" for each mapping and NOT each unique read. This means that you cannot simply use those tools to determine how many reads were mapped. To do that something like the following Linux command will work (though it's likely rather slow):

samtools view accepted_hits.bam | cut -f1 | sort | uniq | wc -l

Normalization

There are various methods of normalizing counts from RNA-Seq data. Sequencing depth, gene length, and count distribution are the main biases that must be accounted for in the normalization and/or differential expression calculations.

• RPKM (Mortazavi et al., 2008) - Counts divided by transcript length (kb) times the total number of millions of mapped reads.

          number or reads in the region

          -----------------------------

RPKM =    total reads       region length

          -----------   X   -------------

            1000000             1000

• FPKM (Trapnell et al., 2010) - Fragments per Kilobase of exon per Million mapped reads, analogous to RPKM.
• Upper-quartile (Bullard et al., 2010) - Counts are divided by upper-quartile of counts for transcripts with at least one read.
• TMM (Robinson and Oshlack, 2010) - Trimmed mean of M values.
• Quantiles, as in microarray normalization (Irizarry et al., 2003).

BEDTools and htseq-count can be used to get raw read counts.

Differential Expression

Parametric - Model counts using known probability distributions such as Poisson, Binomial, Negative Binomial, etc.

edgeR (Robinson et al. 2010) - Exact test based on negative binomial

DESeq (Anders and Huber, 2010) - Exact test based on negative binomial

DEGseq (Wang et al., 2010) - MA-plot based methods assuming normal distribution for M/A

baySeq (Hardcastle et al. 2010) - Estimation of posterior likelihood of differential expression via empirical Baysian methods using Poisson or NB distribution.

Non-Parametric - No assumptions about data distribution are made

Fisher's exact test (use normalized counts)

Cuffdiff (Trapnell et al., 2010) - Based on entropy divergence for relative transcript abundances. Divergence is a measurement of the "distance" between the relative abundances of transcripts between two conditions.

NOISeq - Compare signal and noise distributions to determine differentially expressed genes

Handling Replicates - Biological replicates are generally considered to be required in RNA-Seq

You can also determine normalized quantities for each sample and generate scatterplots to show correlation between each sample and produce a correlation matrix to determine the extent of correlation between replicates.

Most differential expression packages handle replicates explicitly.

Visualization

• CummeRbund - Advanced visualization of Cufflinks and Cuffdiff output using R
• IGV - Excellent Java based genome browser for viewing BAM files from the Broad Institute
• IGB - Another Java based genome browser that excels at viewing annotation data has can also view BAM files
• UCSC Genome Browser - Online genome browser that provides access to many types of annotation and existing data and allows you upload your own (even view directly from galaxy).
• Organism specific genome browsers - Many organisms provide users with custom genome browsers with the latest annotation and experimental data available online. Most also allow users to upload their own data.
  • Wormbase - Nematode Information Resource
  • PlasmoDB - Plasmodium Genomics Resource
  • SGD - Saccharomyces Genome Database

Annotation Resources

Finding gene annotations in well-formed GTF files can be a challenge, espcially if you wish to perform differential splicing analysis. You can find some explanation of the issue by Cole Trapnell in this SeqAnswers post. Below are some helpful resources.

• Illumina's iGenomes is an excellent source of such annotation data for many organisms. Illumina has provided the RNA-Seq user community with a set of genome sequence indexes (including Bowtie indexes) as well as GTF transcript annotation files. These files can be used with TopHat and Cufflinks to quickly perform expression analysis and gene discovery. The annotation files are augmented with the tss_id and p_id GTF attributes that Cufflinks needs to perform differential splicing, CDS output, and promoter user analysis. More information about Illumina's iGenomes project.
  • Many of the iGenomes files are already available in the Princeton Galaxy server's Data Libraries. Look in Shared Data -> Data Libraries -> Reference Genomes and Annotations for your genome.
• ENSEMBL provides both FASTA and GTF files for many organisms on their FTP site.
• UCSC is also provides GTF files, however, annotations in these files have a major limitation. The gene_id and transcript_id are identical, making it impossible to collect transcripts from a single gene. Using the UCSC Table Browser with the settings similar to below will provide a valid, but limited GTF file.
  • clade: Mammal
  • genome: Human
  • assembly: Feb. 2009 (GRCh37/hg19)
  • group: Genes and Gene Prediction Tracks
  • track: RefSeq Genes
  • table: refFlat
  • output format: GTF - gene transfer format
• UCSC GenePred to GTF conversion. A better way to get the UCSC gene annotations into GTF format is to use the genePredToGtf conversion program.
• Convert from GenBank. When all else fails, you can use BioPerl scripts (genbank2gff.pl or genbank2gff3.pl) to convert a GenBank record into a GFF file. From there, you might be able to convert it to a GTF file, though it will take at least some scripting and likely some hand tweaking to get it all right.

Other Resources

Overviews

• Review of RNA-Seq Analysis: Oshlack, A., Robinson, M. D. & Young, M. D. From RNA-seq reads to differential expression results. Genome Biology 11, 220 (2010).
• RNA-Seq Blog - Transcriptome Sequencing News
• The ENCODE Project has finalized Standards, Guidelines and Best Practices for RNA-Seq V1.0 (June 2011).
• The Lewis-Sigler Sequencing Core Facility has an excellent RNA-Seq Field Guide which provides an overview of the many variations of RNA-Seq library preparation, each suited to different experiments.

Tutorials

• UAB Galaxy RNA Seq Step by Step Tutorial - Nice step by step tutorial using Galaxy that explains the output files generated from Tophat/Cufflinks
• RNA-Seq Analysis Exercise in Galaxy by Jeremy Goecks - A walk through of RNA-Seq analysis using Galaxy tools. Includes the data used in the analysis, so you can try it yourself and explore further.
• How-to/RNASeq analysis - Guide to running a more customized analysis from the command line using various tools, including Bowtie and R.
• UC Davis RNA-Seq Course - Tutorial covering a few different workflows using including some QA and contamination removal as well as differential expression analysis.

Tuxedo Tools (Bowtie, Tophat, Cufflinks, etc.)

• Tophat short-read RNA mapper
• Cufflinks transcript assembler
• Tophat and Cufflinks Workflow: Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. (2012) Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7(3), 562-78.
• Tophat-cuff[links,compare,diff] tutorial (command-line) http://cufflinks.cbcb.umd.edu/tutorial.html

Other Tools

• MISO - Probabilistic analysis and design of RNA-Seq experiments for identifying isoform regulation
  • This software is installed on the cetus cluster (access restricted to Princeton researchers).
• edgeR (Robinson et al. 2010)
• DESeq (Anders and Huber, 2010)
• DEGseq (Wang et al., 2010)
• baySeq (Hardcastle et al. 2010)
• NOISeq - Estimation and removal of RNA-Seq noise, non-parametric differential expression
• Trinity - de novo transcriptome assembly
• HTSeq - Python package and scripts to work with high throughput sequencing data including fastq, gff/gtf, and vcf formatted files