Mercurial > repos > elixir-it > vinyl_optimizer

comparison GENEO_VINYL.R @ 0:4c6529d120c3 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author elixir-it
date Tue, 09 Jun 2020 16:07:19 +0000
parents
children 35c308dd6420
comparison
equal deleted inserted replaced
-1:000000000000 0:4c6529d120c3
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
55 Fodds=F$estimate
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");
71 cat(paste(outV,surv1,totR,surv2,totT,pval,rat,score,"\n",sep="\t"),file=ofile,append=T);
72 return(score*-1);
73 }
74 library(genalg)
75
76 G=rbga(stringMin=minsV,stringMax=maxV,popSize=20,iters=25,evalFunc=evalVINYL)
77 #cat(summary(G),file=ofile,append=T)
78