|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2 args <- commandArgs(trailingOnly = TRUE)
|
|
|
3 file1=args[1]
|
|
|
4 file2=args[2]
|
|
|
5 config=args[3]
|
|
|
6 minsV=as.vector(as.integer(unlist(strsplit(args[4], split=" "))))
|
|
|
7 maxV=as.vector(as.integer(unlist(strsplit(args[5], split=" "))))
|
|
|
8 ofile=args[6]
|
|
|
9
|
|
|
10 header=c("disease_clinvar","score_AF","score_functional","score_NS","score_nIND","scoreeQTL","scoreG","scoreT","scoreM","scoreR","scoreSP","scoreGW","survP","totP","survA","totA","p-value","ratio","score");
|
|
|
11
|
|
|
12 cat(paste(header,collapse="\t"),"\n",file=ofile,append=T);
|
|
|
13 evalVINYL=function(x)
|
|
|
14 {
|
|
|
15
|
|
|
16 scoreN=c("-disease_clinvar","-score_AF","-score_functional",
|
|
|
17 "-score_NS","-score_nIND","-scoreeQTL",
|
|
|
18 "-scoreG","-scoreT","-scoreM",
|
|
|
19 "-scoreR","-scoreSP","-scoreGW")
|
|
|
20
|
|
|
21 params=unlist(paste(scoreN,x,sep=" ",collapse=" "))
|
|
|
22 onameprefix=as.integer(runif(1)*10000000)#paste(x,collapse="_");
|
|
|
23 oname_files1=unlist(paste(c("-ofile","-ovcfile","-osummary"),
|
|
|
24 paste(rep(file1,3), rep(onameprefix,3), c("ofile","ovcfile","osummary"),sep="."),sep=" ",collapse=" "))
|
|
|
25 oname_files2=unlist(paste(c("-ofile","-ovcfile","-osummary"),
|
|
|
26 paste(rep(file2,3), rep(onameprefix,3), c("ofile","ovcfile","osummary"),sep="."),sep=" ",collapse=" "))
|
|
|
27
|
|
|
28 command1=paste("perl ./score_complete_alt_M.pl -vcf",file1,config,params,oname_files1,sep=" ",collapse=" ")#,oname_files1)#,config)#,params)
|
|
|
29 command2=paste("perl ./score_complete_alt_M.pl -vcf",file2,config,params,oname_files2,sep=" ",collapse=" ")#,oname_files2)#,config)#,params)
|
|
|
30 Res1=system(command1,intern=FALSE)
|
|
|
31 Res2=system(command2,intern=FALSE)
|
|
|
32 fileR=paste(file1,onameprefix,"ofile",sep=".",collapse="")
|
|
|
33 fileT=paste(file2,onameprefix,"ofile",sep=".",collapse="")
|
|
|
34
|
|
|
35 data_R=read.table(fileR,header=T)
|
|
|
36 data_T=read.table(fileT,header=T)
|
|
|
37 range= rev(seq(min(data_R$VINYL_score),max(data_R$VINYL_score),0.5))
|
|
|
38 m=matrix(ncol=2,nrow=2)
|
|
|
39 totR=nrow(data_R)
|
|
|
40 totT=nrow(data_T)
|
|
|
41 score=0;
|
|
|
42 surv1=0;
|
|
|
43 surv2=0;
|
|
|
44 rat=0;
|
|
|
45 pval=1;
|
|
|
46 for (r in range)
|
|
|
47 {
|
|
|
48 posR=sum(data_R$VINYL_score>=r);
|
|
|
49 posT=sum(data_T$VINYL_score>=r)+1;
|
|
|
50 m[,1]=c(posR,totR);
|
|
|
51 m[,2]=c(posT,totT);
|
|
|
52
|
|
|
53 F=fisher.test(m,alternative="greater")
|
|
|
54 Fpv=F$p.value
|
|
3
|
55 Fodds=F$estimate
|
|
0
|
56 localScore=0.5*-log10(Fpv)+0.3*Fodds-0.2*posT #+0.175*posR-0.125*posT
|
|
|
57 if (localScore>score)
|
|
|
58 {
|
|
|
59
|
|
|
60 score=localScore
|
|
|
61 pval=Fpv
|
|
|
62 rat=Fodds
|
|
|
63 surv1=posR
|
|
|
64 surv2=posT
|
|
|
65
|
|
|
66 }
|
|
|
67
|
|
|
68 }
|
|
|
69 Command=system("rm *.ofile *.ovcfile *.osummary",intern=FALSE)
|
|
|
70 outV=paste(round(x,digits=2),collapse="\t");
|
|
3
|
71 if (surv1>60)
|
|
|
72 {
|
|
|
73 surv1=62;
|
|
|
74 }
|
|
0
|
75 cat(paste(outV,surv1,totR,surv2,totT,pval,rat,score,"\n",sep="\t"),file=ofile,append=T);
|
|
|
76 return(score*-1);
|
|
|
77 }
|
|
|
78 library(genalg)
|
|
|
79
|
|
3
|
80 G=rbga(stringMin=minsV,stringMax=maxV,popSize=50,iters=50,evalFunc=evalVINYL)
|
|
0
|
81 #cat(summary(G),file=ofile,append=T)
|
