HOMER

Software for motif discovery and next-gen sequencing analysis

Analysis of RNA-Seq data with R/Bioconductor

Installing Bioconductor and packages in R

To install R, go to the R homepage and install the appropriate version for your computer (CRAN download page).

Go here to get a full description about how what bioconductor is and how to install it (below is the cheat sheet):
After you start R, type:

# load the script from the internet that is used in install bioconductor
source("http://bioconductor.org/biocLite.R")

# Each of these commands tells Bioconductor to download and install each package
biocLite("GenomicRanges")
biocLite("GenomicFeatures")
biocLite("Rsamtools")
biocLite("DESeq")
biocLite("edgeR")
biocLite("org.Mm.eg.db")

Couple of R basics

Main R tutorial

This one isn't too bad either.

Also, each bioconductor package has it's own tutorial/documentation (usually they offer a lot of explanation)

Also, to change your current working directory in R, got to the top menu "Misc -> Change Working Directory"

Counting Reads in annotated genes (using R)

Before checking differential expression, we need to quantify reads in transcripts.

# load library for genomic annotations
library(GenomicFeatures)
library(GenomicRanges)

# load the transcript annotation file from UCSC.  Make sure to enter the correct genome version
txdb=makeTranscriptDbFromUCSC(genome='mm9',tablename='refGene')

# Use the function transcriptsBy(txdb,'gene') for the whole genic region instead of just exons
ex_by_gene=exonsBy(txdb,'gene')


# load the samtools library for R
library(Rsamtools)

# read the sequencing read alignment into R (combine with next step to save memory)
reads1r1=readBamGappedAlignments("exp1.rep1.bam")
reads1r2=readBamGappedAlignments("exp1.rep2.bam")
reads2r1=readBamGappedAlignments("exp2.rep1.bam")
reads2r2=readBamGappedAlignments("exp2.rep2.bam")
#repeat as necessary for more samples)

# count reads overlapping the exons
counts1r1 = countOverlaps(ex_by_gene,reads1r1)
counts1r2 = countOverlaps(ex_by_gene,reads1r2)
counts2r1 = countOverlaps(ex_by_gene,reads2r1)
counts2r2 = countOverlaps(ex_by_gene,reads2r2)

# create count table
countTable = data.frame(condition1r1=counts1r1,condition1r2=counts1r2,condition2r1=counts2r1,
    condition2r2=counts2r2,stringsAsFactors=FALSE)
# set the gene IDs to the table row names
rownames(countTable)=names(ex_by_gene)

#output tag counts to a file
write.table(countTable,file="countTable.txt",sep="\t")

#removing rows that are zero for all genes (edgeR and DESeq have trouble with these)
x <- rowSums(countTable==0)!=ncol(countTable)
newCountTable <- countTable[x,]

Differential Expression Analysis with EdgeR

# install edgeR
biocLite("edgeR")

# load the edgeR library
library(edgeR)

# either use the input data from above, or load your own table of read counts.
data <- as.matrix(read.table("countTable.txt"))
# - or -
data = newCountTable

# assign groups to the samples in the counts table
g <- c(0,0,1,1)
# get the library sizes (total counts in annotated genes)
libSizes <- as.vector(colSums(data))


# edgeR command pipeline (basically the same for each sample)
d <- DGEList(counts=data,group=g,lib.size=libSizes)
d <- calcNormFactors(d)
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
de.com <- exactTest(d)
results <- topTags(de.com,n = length(data[,1]))

# write the output to a text file
write.table(as.matrix(results$table),file="outputFile.txt",sep="\t")



# No replicates
# If you don't have replicates, edgeR will throw a fit.  In place of:
d <- estimateCommonDisp(d)
d <- estimateTagwiseDisp(d)
you'll want to use:
d$common.dispersion <- 0.05

DESeq

# install DESeq
biocLite("DESeq")

# load the DESeq library
library(DESeq)

# create a DESeq analysis object
cds = newCountDataSet(data,c("notx","notx","lps","lps"))

# DESeq pipeline (similar to the edgeR one)
cds = estimateSizeFactors( cds )
cds = estimateDispersions( cds )
res = nbinomTest( cds, "notx", "lps" )

# print out the results
write.table(res,file="DESeqOutput.txt",sep="\t")

Adding Annotation

Lets say you have a table named "dataTable" (must be data table, i.e. dataTable <- as.data.table(x)).

# for human, replace org.Mm.eg.db with org.Hs.eg.db
biocLite("org.Mm.eg.db")
library(org.Mm.eg.db)

# Use this command to see which types of IDs you can convert:
keytypes(org.Mm.eg.db)

fbids <- rownames(dataTable)
cols <- "SYMBOL"
annots <- select(org.Mm.eg.db,keys=fbids,columns=cols,keytype="ENTREZID")
dataTable <- cbind(entrezid=rownames(as.data.frame(dataTable)), as.data.frame(dataTable))
newTable = merge(dataTable,annots,by.x="entrezid",by.y="ENTREZID")

write.table(newTable,file="annotatedOutput.txt",sep="\t")