|
1
|
1 #!/usr/bin/env Rscript
|
|
|
2 .libPaths('~/R_libs_36')
|
|
|
3 args = commandArgs(trailingOnly=TRUE)
|
|
|
4
|
|
|
5 library(plyr)
|
|
|
6 library("dplyr", warn.conflicts = FALSE)
|
|
|
7 library(reshape2)
|
|
|
8 library(readr)
|
|
|
9 getwd()
|
|
|
10
|
|
|
11 exp<-read.delim(args[1],sep=" ",row.names = 1)
|
|
|
12 samples<-read.delim(args[2],sep="\t",row.names = 1)
|
|
|
13 exp<-exp %>% dplyr::filter(row.names(exp) %in% row.names(samples))
|
|
|
14 desired<-args[3]
|
|
|
15 undesired<-args[4]
|
|
|
16 filter<-args[5]
|
|
|
17
|
|
|
18
|
|
|
19 dfdes<- exp %>% filter(rownames(exp) %in% rownames(samples[samples$phenotype %in% desired,,drop=FALSE]))
|
|
|
20 meandes<- apply(dfdes,2,mean)
|
|
|
21 dfundes<- exp %>% filter(rownames(exp) %in% rownames(samples[samples$phenotype %in% undesired,,drop=FALSE]))
|
|
|
22 meandundes<-apply(dfundes,2,mean)
|
|
|
23 test<-meandes-meandundes
|
|
|
24 disc<-test
|
|
|
25 disc[disc<0]<-0
|
|
|
26 disc[disc>0]<-1
|
|
|
27 nodes<-colnames(dfdes)
|
|
|
28
|
|
|
29 is<-read.delim(args[6],sep=" ") #internal_markers.txt
|
|
|
30 isdisc<-is
|
|
|
31 if (filter=="strict"){
|
|
|
32 for (node in nodes){
|
|
|
33 if (disc[node]==0){
|
|
|
34 val=max(dfdes[,node])
|
|
|
35 isdisc[,node]<-replace(isdisc[,node],is[,node]<val, 0)
|
|
|
36 isdisc[,node]<-replace(isdisc[,node],is[,node]>val, 1)
|
|
|
37 }
|
|
|
38 if (disc[node]==1){
|
|
|
39 val=min(dfdes[,node])
|
|
|
40 isdisc[,node]<-replace(isdisc[,node],is[,node]>val, 1)
|
|
|
41 isdisc[,node]<-replace(isdisc[,node],is[,node]<val, 0)
|
|
|
42 }
|
|
|
43 }
|
|
|
44 }else { #filter=="relaxed"
|
|
|
45 for (node in nodes){
|
|
|
46 if (disc[node]==0){
|
|
|
47 val=min(dfundes[,node])
|
|
|
48 isdisc[,node]<-replace(isdisc[,node],is[,node]<val, 0)
|
|
|
49 isdisc[,node]<-replace(isdisc[,node],is[,node]>val, 1)
|
|
|
50 }
|
|
|
51 if (disc[node]==1){
|
|
|
52 val=max(dfdes[,node])
|
|
|
53 isdisc[,node]<-replace(isdisc[,node],is[,node]>val, 1)
|
|
|
54 isdisc[,node]<-replace(isdisc[,node],is[,node]<val, 0)
|
|
|
55 }
|
|
|
56 }
|
|
|
57
|
|
|
58 }
|
|
|
59
|
|
|
60
|
|
|
61
|
|
|
62 hammings<-data.frame(isdisc$replicate,isdisc$perturbation)
|
|
|
63 colnames(hammings)=c("replicate","perturbation")
|
|
|
64 hammings$hammingdist<-sapply(isdisc[,3],function(x) sum(disc != x))
|
|
|
65
|
|
|
66 if (length(nodes)<10){
|
|
|
67 hams1_1<-hammings %>%filter(hammingdist<1)
|
|
|
68 }else{
|
|
|
69 hams1_1<-hammings %>%filter(hammingdist<=length(nodes)-round(.9*length(nodes)))
|
|
|
70 }
|
|
|
71
|
|
|
72 write.table(hams1_1[,c(1,2)],'crit2_perturbations.txt',quote=FALSE,row.names = FALSE,sep=" ")
|
|
|
73
|
|
|
74
|
