...
| Code Block | ||
|---|---|---|
| ||
cds cd my_rnaseq_course cd gene_expression_exercise # you should have already copied this over #modify our gene counts file a bit grep '^F' gene_counts.htseq.gff > gene_counts.htseq.final.gff |
| Warning |
|---|
When following along here, please switch to your idev session for running these example commands. |
Introduction
Most RNA-Seq experiments are conducted with the aim of identifying genes/exons that are differentially expressed between two or more conditions. Many computational tools are available for performing the statistical tests required to identify these genes/exons.
...
| Code Block | ||
|---|---|---|
| ||
login1$ R
> library("DESeq")
> library("edgeR")
|
These commands will work for any Bioconductor package!
DESeq
Input:
DESeq takes as input count data in table form, with each column representing a biological replicate/biological condition. The count data must be raw counts of sequencing reads, not already normalized data.
Example:
untreated1 untreated2 untreated3 untreated4 treated1 treated2 treated3
FBgn0000003 0 0 0 0 0 0 1
FBgn0000008 92 161 76 70 140 88 70
FBgn0000014 5 1 0 0 4 0 0
Our gene_counts.gff file is almost there- we can remove some extra gene related information in the counts file to make it less complex.
| Code Block | ||
|---|---|---|
| ||
login1$ R
#LOAD LIBRARY, READ IN GENE COUNTS, ADD METADATA
> library("DESeq")
> counts = read.delim("gene_counts.htseq.final.gff", header=F, row.names=1)
> head(counts)
> colnames(counts) = c("C11","C12","C13","C21","C22","C23")
> head(counts)
> my.design <- data.frame(
row.names = colnames( counts ),
condition = c( "C1", "C1", "C1", "C2", "C2", "C2"),
libType = c( "paired-end", "paired-end", "paired-end", "paired-end", "paired-end", "paired-end" )
)
> my.design
> conds <- factor(my.design$condition)
#CREATE A DESEQ DATA OBJECT FROM COUNTS
> cds <- newCountDataSet( counts, conds )
> cds
#NORMALIZATION: ESTIMATE SIZE FACTORS
> cds <- estimateSizeFactors( cds )
> sizeFactors( cds )
> head( counts( cds, normalized=TRUE ) )
#ESTIMATE DISPERSION/VARIANCE
> cds <- estimateDispersions( cds )
#DO TEST FOR DIFFERENTIAL EXPRESSION AND WRITE RESULTS INTO FILE
> result <- nbinomTest( cds, "C1", "C2" )
> head(result)
> result = result[order(result$pval), ]
> head(result)
> write.csv(result, "DESeq-C1-vs-C2.csv")
#GENERATE MA PLOT
> pdf("DESeq-MA-plot.pdf")
> plot(
result$baseMean,
result$log2FoldChange,
log="x", pch=20, cex=.3,
col = ifelse( result$padj < .1, "red", "black" ) )
> dev.off()
> q()
Save workspace image? [y/n/c]: n
login1$ head DESeq-wt-vs-mut.csv
|
edgeR
These commands use the negative binomial model, calculate the false discovery rate (FDR ~ adjusted p-value), and make a plot similar to the one from DESeq.
| Code Block | ||
|---|---|---|
| ||
login1$ R
#LOAD LIBRARY, READ IN COUNTS, ADD SOME METADATA
> library("edgeR")
> counts = read.delim("gene_counts.htseq.gff", header=F, row.names=1)
> colnames(counts) = c("C11", "C12", "C13", "C21", "C22", "C23")
> head(counts)
> group <- factor(c("C1", "C1", "C1", "C2", "C2", "C2"))
#CREATE EDGER SPECIFIC DATA OBJECT FROM COUNTS DATA
> dge = DGEList(counts=counts,group=group)
#ESTIMATE DISPERSIONS(VARIANCE)
> dge <- estimateCommonDisp(dge)
> dge <- estimateTagwiseDisp(dge)
#DO TEST FOR DIFFERENTIAL EXPRESSION
> et <- exactTest(dge)
> etp <- topTags(et, n=100000)
#GENERATE MA PLOT
> pdf("edgeR-MA-plot.pdf")
> plot(
etp$table$logCPM,
etp$table$logFC,
xlim=c(-3, 20), ylim=c(-12, 12), pch=20, cex=.3,
col = ifelse( etp$table$FDR < .1, "red", "black" ) )
> dev.off()
#WRITE OUT GENE RESULTS
> write.csv(etp$table, "edgeR-wt-vs-mut.csv")
> q()
Save workspace image? [y/n/c]: y
login1$ head edgeR-wt-vs-mut.csv |
DEXSeq
This package is meant for finding differential exon usage between samples from different conditions.
Relative usage of an exon = transcripts from the gene that contain this exon / all transcripts from the gene
For each exon (or part of an exon) and each sample :
- count how many reads map to this exon
- count how many reads map to other exons of the same gene.
- calculate ratio of 1 to 2.
- Look for changes in this ratio across conditions
- Look for statistically significant changes in this ratio across conditions, by using replicates.
This lets you identify changes in alternative splicing, changes in usage of alternative transcript start sites.
Cuffdiff
Cuffdiff (a part of the tuxedo suite) is a popular tool for testing for differential expression. We will cover this along with the rest of the tuxedo suite tomorrow.
Lets look at our results
| Code Block | ||
|---|---|---|
| ||
head DESeq-wt-vs-mut.csv
head edgeR-wt-vs-mut.csv
|
| Code Block | ||
|---|---|---|
| ||
login1$ R > head DESeq-wt-vs-mut.csv > head edgeR-wt-vs-mut.csv |
The graphs we generated are located at:
http://web.corral.tacc.utexas.edu/BioITeam/rnaseq_course/DESeq-MA-plot.pdf
http://web.corral.tacc.utexas.edu/BioITeam/rnaseq_course/edgeR-MA-plot.pdf