/usr/share/perl5/Genome/Model/Tools/Music/MutationRelation.pm.R is in libgenome-model-tools-music-perl 0.04-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 | #get command line arguments
mutation_file = as.character(commandArgs()[4]);
permutations = as.numeric(commandArgs()[5]);
output_file = as.character(commandArgs()[6]);
# FUNCTION prepare mutation matrix data for permutation
seq.vjoint=function(x)
{
tt=NULL;
for (i in 1:length(x))
{
tt=c(tt,rep(names(x[i]),x[i]));
}
tt;
}
# END seq.vjoint
# FUNCTION randomly permute the mutation matrix, keeping mutation-number/sample unchanged
matrix.sample.v1=function(x,keep.freq=T)
{
colSums(x)->pool;
rowSums(x)->freq;
freq=freq[order(freq,decreasing=T)];
tt=NULL;
for (i in c(1:length(freq)))
{
pooli=seq.vjoint(pool);
xi=NULL;
for (j in 1:freq[i])
{
if (length(pooli)>1) sample(pooli,1)->temp;
if (length(pooli)==1) pooli->temp;
xi=c(xi,temp);
pooli=pooli[pooli!=temp];
}
tt=rbind(tt,cbind(i,xi));
pool[xi]=pool[xi]-1;
pool=pool[pool>0];
}
tt=table(tt[,1],tt[,2]);
tt[tt>1]=1;
tt;
}
# END matrix.sample.v1
# FUNCTION calculate the probability of having no correlation between any two genes
prob.concur=function(x)
{
x>0->x;
cols=colnames(x);
cols=cols[order(cols)];
n=length(cols);
tt=NULL;
for (i in 1:(n-1))
{
for (j in (i+1):n)
{
vi=cols[i];
vj=cols[j];
sum(x[,vi] & x[,vj])->nand;
sum(x[,vi] != x[,vj])->nexc;
ni=sum(x[,vi]);
nj=sum(x[,vj]);
temp=cbind(vi,vj,ni,nj,nand,nexc);
tt=rbind(tt,temp);
}
}
tt=as.data.frame(tt,stringsAsFactors=F);
tt$nand=as.numeric(tt[,"nand"]);
tt$nexc=as.numeric(tt[,"nexc"]);
tt;
}
# END prob.concur
# FUNCTION performing mutation correlation
mut_cor_permu_test=function(x=NULL,n.permu=100,seed=NULL,mut.file=NULL,out.file=NULL,out.rdata=NULL)
{
if (!is.null(mut.file)) read.table(mut.file,header=T)->x
x=x[,-1];x[x>0]=1
x=x[rowSums(x)>0,colSums(x)>0]
prob.concur(x)->pc0
if (!is.null(seed)) set.seed(seed)
pp=0;en=0
for (i in 1:n.permu)
{
doit=1;
while(doit==1)
{
xi=matrix.sample.v1(x);
if (max(xi)==1) doit=0;
}
prob.concur(xi)->pci;
en=en+pci[,c("nand","nexc")];
as.numeric(pci$nand>=pc0$nand)->pand;
as.numeric(pci$nexc>=pc0$nexc)->pexc;
pp=pp+cbind(pand,pexc);
}
en=en/n.permu;pp=pp/n.permu;
pci[,c("nand","nexc")]=en;
pci=cbind(pci,pp);
tt=merge(pc0,pci,by.x=c("vi","vj"),by.y=c("vi","vj"));
if (!is.null(out.rdata)) save(tt,file=out.rdata,compress=T);
names(tt)= c("Gene1","Gene2","CntGene1","CntGene2","AndCnt","XorCnt","Perm_CntGene1","Perm_CntGene2","Perm_AndCnt","Perm_XorCnt","Pvalue_And","Pvalue_Xor")
if (!is.null(out.file)) write.table(tt,file=out.file,quote=F,sep="\t",row.names=F);
invisible (tt);
}
#END mut_cor_permu_test
#run test using mut_cor_permu_test
mut_cor_permu_test(mut.file=mutation_file,n.permu=permutations,out.file=output_file);
|