Mercurial > repos > ning > ebseq

comparison EBIsoformTwoCondTest.R @ 38:52d502275e1f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author ning
date Mon, 16 Jun 2014 16:23:41 -0400
parents
children
comparison
equal deleted inserted replaced
37:e11c511a2bb5 38:52d502275e1f
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 Ig.file <-args[4]
14 FDR <- args[5]
15 outputfile <- args[6]
16 Sort.out<-args[7]
17 Sort.out.FDR <-args[8]
18 Sizesout <-args[9]
19
20
21 Conditions=strsplit(CondIn,split=",")[[1]]
22 if(WhetherSampleName=="y"){
23 ReadIn=read.table(inputfile,stringsAsFactors=F,header=T,sep="\t")
24 Names=names(ReadIn)[-1]
25 }
26 if(WhetherSampleName=="n"){
27 ReadIn=read.table(inputfile,stringsAsFactors=F,header=F, sep="\t")
28 Names=paste0("S",1:length(Conditions))
29 }
30
31 if(class(ReadIn[[1]])=="character"){
32 GeneMat=do.call(cbind,ReadIn[-1])
33 rownames(GeneMat)=ReadIn[[1]]
34 colnames(GeneMat)=Names
35 }
36 if(class(ReadIn[[1]])=="numeric"){
37 GeneMat=data.matrix(ReadIn)
38 colnames(GeneMat)=Names
39 }
40
41 IgVIn=read.csv(Ig.file,stringsAsFactors=F,header=F)
42 IgV=IgVIn[[1]]
43 Conditions=strsplit(CondIn,split=",")[[1]]
44
45 Sizes=MedianNorm(GeneMat)
46 EBOut=EBTest(Data=GeneMat,NgVector=IgV,Conditions=as.factor(Conditions),sizeFactors=Sizes, maxround=5)
47 PP=GetPP(EBOut)
48 PP.sort=sort(PP,decreasing=T)
49 PP.sort.FDR=PP.sort[which(PP.sort>=1-as.numeric(FDR))]
50
51 Data.norm=GetNormalizedMat(GeneMat, Sizes)
52 FC=PostFC(EBOut)
53 realFC=FC[[2]]
54 postFC=FC[[1]]
55
56 Mat=cbind(PP, realFC[names(PP)], postFC[names(PP)],Data.norm[names(PP),])
57 Mat.sort=cbind(PP.sort, realFC[names(PP.sort)], postFC[names(PP.sort)],Data.norm[names(PP.sort),])
58
59
60 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),])
61
62 if(length(PP.sort.FDR)==1)Mat.sort.FDR=matrix(
63 c(PP.sort.FDR, realFC[names(PP.sort.FDR)], postFC[names(PP.sort.FDR)],Data.norm[names(PP.sort.FDR),])
64 ,nrow=1)
65
66 colnames(Mat)=colnames(Mat.sort)=
67 c("PPDE","RealFC","PosteriorFC",colnames(Data.norm))
68 if(length(PP.sort.FDR)>0)colnames(Mat.sort.FDR)=
69 c("PPDE","RealFC","PosteriorFC",colnames(Data.norm))
70
71 #write.table(round(Mat,2),file=outputfile,quote=F,col.names=T,row.names=T,sep = "\t")
72 #write.table(round(Mat,2),file=outputfile,quote=F,col.names=NA,row.names=T,sep = "\t") ##cms - changed on 11/15/13
73 write.table(round(Mat,2),file=outputfile,quote=F,col.names=c("Gene/IsoformID\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - changed on 11/15/13
74 #write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=T,row.names=T,sep = "\t")
75 #write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=NA,row.names=T,sep = "\t") ##cms - changed on 11/15/13
76 write.table(round(Mat.sort,2),file=Sort.out ,quote=F,col.names=c("Gene/IsoformID\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - changed on 11/15/13
77 #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")
78 #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 - changed on 11/15/13
79 if(length(PP.sort.FDR)>0)write.table(round(Mat.sort.FDR,2),file=Sort.out.FDR,quote=F,col.names=c("Gene/IsoformID\tPPDE","RealFC","PosteriorFC",colnames(Data.norm)),row.names=T,sep = "\t") ##cms - changed on 11/15/13
80 write.table(Sizes,file=Sizesout,quote=F,col.names=F,row.names=F,sep = "\t")
81