Page Comparison

Objectives

Once we've obtained abundance counts for our genes/exons/transcripts, we are usually interested in identifying those genes/exons/transcripts that are differentially expressed. 

...

• Learn about DESeq, edgeR, DEXSeq and cuffdiff packages and the differences among these packages.
• Become familiar with basic R usage and installing Bioconductor modules.
• Learn how to use edgeR/DESeq to identify differentially expressed genes.
• Learn how to use cuffdiff pacakge to identify differentially expressed genes.

Get set up

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

...

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.

...

• Normalization of gene counts
  
  
• Represent the gene counts by a distribution that defines the relation between mean and variance (dispersion).
  
  
• Perform a statistical test on this distribution to identify genes that are significantly different between the conditions.
  
  
• Provide fold change,  P-value information, false discovery rate for each gene.

Why Normalize?

Normalization smooths out technical variations among the samples we are comparing so that we can more confidently attribute variations we see to biological reasons. 

...

1. Sequencing depth: Say we are comparing gene counts in sample A against sample B.  If you start out with 10 million reads in sample A  vs 1 million reads in sample B, a 10 fold increase in expression in sample A is going to be purely due to its sequencing depth.
2. Gene length: A gene that is twice as long is likely to have twice as many reads sampling it.

Most commonly done normalization

RPKM: Normalizes for sequencing depth and gene length.

...

DEG tools compared

From paper: Wesolowski et al, Biosensors, 2013

Image Modified

R and Bioconductor, very briefly...

R is a very common scripting language used in statistics. There are whole courses on using R going on in other SSI classrooms as we speak! Inside the R universe, you have access to an incredibly large number of useful statistical functions (Fisher's exact test, nonlinear least-squares fitting, ANOVA ...). R also has advanced functionality for producing plots and graphs as output. 

...

Hints for working with R

• Don't forget: it's q() to quit.
• For help with a function, type ?command. Try ?read.table. The q key gets you out of help, just like for a man page.
• The left arrow <- (less-than-dash) is the same as an equals sign =. You can use them interchangeably.
• The prompt we will sometimes be showing for R is >. Don't type this for a command. It is like the login1$ at the beginning of the bash prompt when you log in to Lonestar. It just means that you are in the R shell.

...

> x <- 10 + 5 + 6
> x
[1] 21

Bioconductor packages for R

Like other languages, R can be expanded by loading packages. The R equivalent of Bioperl or Biopython is Bioconductor. Bioconductor can theoretically do things for you like convert sequences (none of us use it for that), but where it really shines is in doing statistical tests (where is it second-to-none in this list of languages). Many functions for analyzing microarray data are implemented in R, and this strength has now carried over to the analysis of RNAseq data.

...

login1$ R
> library("DESeq")
> library("edgeR")

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

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

head DESeq-wt-vs-mut.csv
head edgeR-wt-vs-mut.csv

#DESeq results
sed 's/,/\t/g' DESeq-C1-vs-C2.csv|sort -n -r -k6,6|cut -f 2,6|head

#edgeR results
sed 's/,/\t/g' edgeR-wt-vs-mut.csv |sort -n -r -k2,2|cut -f 1,2|head

#We will use log 2 fold change instead of fold change because edgeR doesnt report fold change
#DESeq results
sed 's/,/\t/g' DESeq-C1-vs-C2.csv|awk '{if (($7>=1)&&($8<=0.05)) print $1,$7,$8}'|head
sed 's/,/\t/g' DESeq-C1-vs-C2.csv|awk '{if (($7>=1)&&($8<=0.05)) print $1,$7,$8}'|wc -l

#edgeR results
sed 's/,/\t/g' edgeR-wt-vs-mut.csv |awk '{if (($2>=1)&&($3<=0.05)) print $1,$2,$3}'|head
sed 's/,/\t/g' edgeR-wt-vs-mut.csv |awk '{if (($2>=1)&&($3<=0.05)) print $1,$2,$3}'|wc -l