|
34
|
1 sink(file="/tmp/none")
|
|
|
2 sink("/dev/null")
|
|
|
3 options(warn=-1)
|
|
|
4 options(echo=F)
|
|
|
5
|
|
|
6 invisible("EBSeq")
|
|
|
7 suppressMessages(library("EBSeq"))
|
|
|
8
|
|
|
9 args <- commandArgs(trailingOnly = T)
|
|
|
10 inputfile <- args[1]
|
|
|
11 WhetherSampleName <- args[2]
|
|
|
12 CondIn <- args[3]
|
|
|
13 FDR <- args[4]
|
|
|
14 outputfile <- args[5]
|
|
|
15 Sort.out<-args[6]
|
|
|
16 Sort.out.FDR <-args[7]
|
|
|
17 Sizesout <-args[8]
|
|
|
18
|
|
|
19 Conditions=strsplit(CondIn,split=",")[[1]]
|
|
|
20
|
|
|
21 if(WhetherSampleName=="y"){
|
|
|
22 ReadIn=read.table(inputfile,stringsAsFactors=F,header=T,sep="\t")
|
|
|
23 Names=names(ReadIn)[-1]
|
|
|
24 }
|
|
|
25 if(WhetherSampleName=="n"){
|
|
|
26 ReadIn=read.table(inputfile,stringsAsFactors=F,header=F,sep="\t")
|
|
|
27 Names=paste0("S",1:length(Conditions))
|
|
|
28 }
|
|
|
29
|
|
|
30 if(class(ReadIn[[1]])=="character"){
|
|
|
31 GeneMat=do.call(cbind,ReadIn[-1])
|
|
|
32 rownames(GeneMat)=ReadIn[[1]]
|
|
|
33 colnames(GeneMat)=Names
|
|
|
34 }
|
|
|
35 if(class(ReadIn[[1]])=="numeric"){
|
|
|
36 GeneMat=data.matrix(ReadIn)
|
|
|
37 colnames(GeneMat)=Names
|
|
|
38 }
|
|
|
39
|
|
|
40
|
|
|
41 Sizes=MedianNorm(GeneMat)
|
|
|
42 EBOut=EBTest(Data=GeneMat,Conditions=as.factor(Conditions),sizeFactors=Sizes, maxround=5)
|
|
|
43 PP=GetPP(EBOut)
|
|
|
44 PP.sort=sort(PP,decreasing=T)
|
|
|
45 PP.sort.FDR=PP.sort[which(PP.sort>=1-as.numeric(FDR))]
|
|
|
46
|
|
|
47 Data.norm=GetNormalizedMat(GeneMat, Sizes)
|
|
|
48 FC=PostFC(EBOut)
|
|
|
49 realFC=FC[[2]]
|
|
|
50 postFC=FC[[1]]
|
|
|
51
|
|
|
52 Mat=cbind(PP, realFC[names(PP)], postFC[names(PP)],Data.norm[names(PP),])
|
|
|
53 Mat.sort=cbind(PP.sort, realFC[names(PP.sort)], postFC[names(PP.sort)],Data.norm[names(PP.sort),])
|
|
|
54
|
|
|
55
|
|
|
56 if(length(PP.sort.FDR)>1)Mat.sort.FDR=cbind(PP.sort.FDR, realFC[names(PP.sort.FDR)], postFC[names(PP.sort.FDR)],Data.norm[names(PP.sort.FDR),])
|
|
|
57
|
|
|
58 if(length(PP.sort.FDR)==1)Mat.sort.FDR=matrix(
|
|
|
59 c(PP.sort.FDR, realFC[names(PP.sort.FDR)], postFC[names(PP.sort.FDR)],Data.norm[names(PP.sort.FDR),])
|
|
|
60 ,nrow=1)
|
|
|
61
|
|
|
62 colnames(Mat)=colnames(Mat.sort)=
|
|
|
63 c("PPDE","RealFC","PosteriorFC",colnames(Data.norm))
|
|
|
64 if(length(PP.sort.FDR)>0)colnames(Mat.sort.FDR)=
|
|
|
65 c("PPDE","RealFC","PosteriorFC",colnames(Data.norm))
|
|
|
66
|
|
|
67 ##cms - col.names changed below to have a header for the row names. This way the output of RSEM can be joined with EBSeq output
|
|
|
68 #write.table(round(Mat,2),file=outputfile,quote=F,col.names=T,row.names=T,sep = "\t")
|
|
|
69 #write.table(round(Mat,2),file=outputfile,quote=F,col.names=NA,row.names=T,sep = "\t") ##cms - made correction on 11/11/13
|
|
|
70 write.table(round(Mat,2),file=outputfile,quote=F,col.names=c("gene_id\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - made correction on 11/11/13
|
|
|
71 #write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=T,row.names=T,sep = "\t")
|
|
|
72 #write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=NA,row.names=T,sep = "\t") ##cms - made correction on 11/11/13
|
|
|
73 write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=c("gene_id\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - made correction on 11/11/13
|
|
|
74 #if(length(PP.sort.FDR)>0)write.table(round(Mat.sort.FDR,2),file=Sort.out.FDR,quote=F,col.names=T,row.names=T,sep = "\t")
|
|
|
75 #if(length(PP.sort.FDR)>0)write.table(round(Mat.sort.FDR,2),file=Sort.out.FDR,quote=F,col.names=NA,row.names=T,sep = "\t") ##cms - made correction on 11/11/13
|
|
|
76 if(length(PP.sort.FDR)>0)write.table(round(Mat.sort.FDR,2),file=Sort.out.FDR,quote=F,col.names=c("gene_id\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - made correction on 12/17/13; this preserves the header if merged with RSEM output.
|
|
|
77 write.table(Sizes,file=Sizesout,quote=F,col.names=F,row.names=F,sep = "\t")
|
|
|
78
|
