############################################################################################################################################ # 2014.03.11 # 2M, 4M, and 6M are replicates. Run the HS-WT comparison through limma library(Ringo) library(limma) todaysDate<-substr(format(Sys.time(), "%Y.%m.%d."),1,11) userName="coloson" RG<-readNimblegen("targets.txt", "spottypes.txt", headerPattern="# software=DEVA", path="pair") RG<-RG[RG$gene$Status != "RANDOM",]; rg<-RG; nrow(RG) #697752 (23835 RANDOM controls were discarded) MA<-normalizeWithinArrays(RG, method="none") #Normalization: use NimbleGen processing method (Tukey biweight mean) NG<-preprocess(RG,method="nimblegen", returnMAList=TRUE) Source<-c(2,2,4,4,6,6) Treatment<-c("wt","hs","wt","hs","wt","hs") design<-model.matrix(~Source+Treatment) fit<-lmFit(NG,design); fit<-eBayes(fit) tt<-topTable(fit, number=nrow(x), adjust.method="BH", coef="Treatmentwt") m<-merge(tt,anno,by.x="PROBE_ID",by.y="ProbeName"); m<-m[,!names(m) %in% c("GENE_EXPR_OPTION","POSITION","X","Y","Status","ID","t","B")]; m<-m[order(m$adj.P.Val),] names(m)[1]<-"ProbeName" #swap sign on ratio so it's HS-WT m$logFC<--m$logFC #add individual ratios to results normData.M<-data.frame(NG$M, check.names=F, row.names=NG$genes$PROBE_ID); names(normData.M)<-gsub("$",".M",names(normData.M)) normData.A<-data.frame(NG$A, check.names=F, row.names=NG$genes$PROBE_ID); names(normData.A)<-gsub("$",".A",names(normData.A)) normData<-data.frame(normData.M, normData.A, check.names=F) names(normData)<-gsub("[0-9]+A0[1-3]_","",names(normData)); names(normData)<-gsub("_meth_635","",names(normData)) normData<-data.frame(row.names(normData),normData,check.names=F); names(normData)[1]<-"ProbeName" m2<-merge(m,normData,by="ProbeName"); m2<-m2[order(m2$adj.P.Val),] write.table(m2,file="2014.03.11.colson_HS-WT_limma.txt",sep="\t",quote=F,row.names=F) #generate bedGraphs of limma results, and just sig limma results todaysDate<-"2014.03.11" x<-data.frame(m[6:8],m[2],check.names=F) trackName<-"limma logFC" i<-"limmaLogFC" trackLine<-paste("track type=bedGraph name=\"", trackName , "\" description=\"", trackName, "\" visibility=full color=100,100,0 altColor=0,100,200 priority=20", sep="") write.table(trackLine,file=paste(todaysDate, userName, "_", i,".bedGraph", sep=""), sep="\t", quote=F, row.names=F, col.names=F) write.table(x,file=paste(todaysDate, userName, "_", i,".bedGraph", sep=""), sep="\t", quote=F, row.names=F, col.names=F, append=TRUE) #sig sig<-m[abs(m$logFC) >=1 & m$adj.P.Val<.05,]; nrow(sig) #211174 x<-data.frame(sig[6:8], sig[2],check.names=F) trackName<-"limma sig logFC" i<-"limmaSigLogFC" trackLine<-paste("track type=bedGraph name=\"", trackName , "\" description=\"", trackName, "\" visibility=full color=100,100,0 altColor=0,100,200 priority=20", sep="") write.table(trackLine,file=paste(todaysDate, userName, "_", i,".bedGraph", sep=""), sep="\t", quote=F, row.names=F, col.names=F) write.table(x,file=paste(todaysDate, userName, "_", i,".bedGraph", sep=""), sep="\t", quote=F, row.names=F, col.names=F, append=TRUE)