library(edgeR) #Read in the file cds <- read.delim("C:/Users/rhorn83/Documents/Transcriptome Seq/T-longicaudatus-20130826/T-longicaudatus-20130826/T-longicaudatus-20130826/readcounts/cds.dat", row.names=1, header=TRUE) #group the correct columns together from the same replicate group <- factor(rep(c("NonNative", "Native"), times=c(3,3))) #makes a DGEList from the count data y <- DGEList(counts=cds,group=group) summarized=DGEList(cds,lib.size=colSums(cds),group=group) #Normalization z <- calcNormFactors(y) #estimate dispersion disp=estimateCommonDisp(summarized) disp$common.dispersion D <- estimateTagwiseDisp(summarized) D$tagwise.dispersion #Test for DE genes tested=exactTest(disp) topTags(tested) #will return a summary of how many genes are up-, down- or normally regulated summary(de <- decideTestsDGE(tested, p=0.05, adjust="BH")) #create an MA plot detags <- rownames(disp)[as.logical(de)] plotSmear(tested, de.tags=detags) abline(h = c(-2, 2), col = "blue")