update bpmsim -- only check similarity based on pathways, not interac…

…tions
wangx677 · Mar 4, 2020 · c4efcaa · c4efcaa
1 parent 7641583
commit c4efcaa
Show file tree

Hide file tree

Showing 3 changed files with 15 additions and 36 deletions.
diff --git a/corefuns/bpmsim.m b/corefuns/bpmsim.m
@@ -1,4 +1,4 @@
-function sim = bpmsim(BPMind1x,BPMind2x,BPMind1y,BPMind2y,option,ssM)
+function sim = bpmsim(BPMind1x,BPMind2x,BPMind1y,BPMind2y)
 
 % FUNCTION sim = bpmsim(BPMind1x,BPMind2x,BPMind1y,BPMind2y,option,ssM)
 %
@@ -10,9 +10,6 @@
 %    BPMind2x - BPMind2 for BPM set X
 %    BPMind1y - BPMind1 for BPM set Y
 %    BPMind2y - BPMind2 for BPM set Y
-%    option - option=1 calculates overlap based on interactions;
-%         option=2 calculates overlap based on SNPs
-%    ssM - binarized SNP-SNP interaction matrix. It is required when option=1
 % OUTPUTS:
 %    sim - pairwise similarity matrix (overlap coefficient)
 
@@ -35,38 +32,20 @@
      BPMind2y = cellsquareform(BPMind2y);
 end
 
-if option==1
 
-	mm =cellfun(@(x,y,z,r) min(nnz(ssM(x,y)),nnz(ssM(z,r))),BPMind1x,BPMind2x,BPMind1y,BPMind2y);
+mm = cellfun(@(x,y,z,r) min(length(x)*length(y),length(z)*length(r)),BPMind1x,BPMind2x,BPMind1y,BPMind2y);
 
-	ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind1y,'UniformOutput',false);
-	ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind2y,'UniformOutput',false);
-
-	simtmp1 = cellfun(@(x,y,z) nnz(ssM(x,y)),ind1,ind2);
-	simtmp1 = simtmp1./mm;
+ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind1y,'UniformOutput',false);
+ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind2y,'UniformOutput',false);
 
-     ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind2y,'UniformOutput',false);
-     ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind1y,'UniformOutput',false);
+simtmp1 = cellfun(@(x,y) (length(x)*length(y)),ind1,ind2);
+simtmp1 = simtmp1./mm;
 
-     simtmp2 = cellfun(@(x,y,z) nnz(ssM(x,y)),ind1,ind2);
-	simtmp2 = simtmp2./mm;	
+ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind2y,'UniformOutput',false);
+ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind1y,'UniformOutput',false);
 
-elseif option==2
-
-	mm = cellfun(@(x,y,z,r) min(length(x)*length(y),length(z)*length(r)),BPMind1x,BPMind2x,BPMind1y,BPMind2y);
-
-	ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind1y,'UniformOutput',false);
-     ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind2y,'UniformOutput',false);
-
-	simtmp1 = cellfun(@(x,y) (length(x)*length(y)),ind1,ind2);
-	simtmp1 = simtmp1./mm;
-
-	ind1 = cellfun(@(x,y) intersect(x,y),BPMind1x,BPMind2y,'UniformOutput',false);
-     ind2 = cellfun(@(x,y) intersect(x,y),BPMind2x,BPMind1y,'UniformOutput',false);
-
-	simtmp2 = cellfun(@(x,y) (length(x)*length(y)),ind1,ind2);
-	simtmp2 = simtmp2./mm;
-end
+simtmp2 = cellfun(@(x,y) (length(x)*length(y)),ind1,ind2);
+simtmp2 = simtmp2./mm;
 
 sim = max(simtmp1,simtmp2);
 

diff --git a/corefuns/check_BPM_WPM_redundancy.m b/corefuns/check_BPM_WPM_redundancy.m
@@ -1,7 +1,6 @@
 function [BPM_nosig_noRD,WPM_nosig_noRD,PATH_nosig_noRD,BPM_group,WPM_group,PATH_group] = check_BPM_WPM_redundancy(fdrBPM,fdrWPM,fdrPATH,bpmindfile,FDRcut)
 
 load(bpmindfile)
-ssM{1} = 0;ssM{2}=0; % arbitrary -- it's not actully used because the overlap is based on pathway genes not actual interactions.  
 j = 1;
 for fdrcut=0.05:0.05:FDRcut
         ind = find(fdrBPM<=fdrcut);
@@ -29,7 +28,7 @@
         end
 
         if (length(ind1)>1)
-                BPM_sim{1} = bpmsim(BPM.ind1(ind1),BPM.ind2(ind1),BPM.ind1(ind1),BPM.ind2(ind1),2,ssM{1});
+                BPM_sim{1} = bpmsim(BPM.ind1(ind1),BPM.ind2(ind1),BPM.ind1(ind1),BPM.ind2(ind1));
                 TTT=BPM_sim{1}>=0.25;
                 [noRD{1},group{1}] = graphconncomp(sparse(TTT));
                 %[num2cell(group') path1(ind1)' path2(ind1)']
@@ -43,7 +42,7 @@
         end
 
         if (length(ind2)>1)
-                BPM_sim{2} = bpmsim(BPM.ind1(ind2),BPM.ind2(ind2),BPM.ind1(ind2),BPM.ind2(ind2),2,ssM{1});
+                BPM_sim{2} = bpmsim(BPM.ind1(ind2),BPM.ind2(ind2),BPM.ind1(ind2),BPM.ind2(ind2));
                 TTT=BPM_sim{2}>=0.25;
                 [noRD{2},group{2}] = graphconncomp(sparse(TTT));
                 %[num2cell(group') path1(ind2)' path2(ind2)']
@@ -57,7 +56,7 @@
         end
 
         if (length(ind3)>1)
-                BPM_sim{3} = bpmsim(WPM.ind(ind3),WPM.ind(ind3),WPM.ind(ind3),WPM.ind(ind3),2,ssM{1});
+                BPM_sim{3} = bpmsim(WPM.ind(ind3),WPM.ind(ind3),WPM.ind(ind3),WPM.ind(ind3));
                 TTT=BPM_sim{3}>=0.25;
                 [noRD{3},group{3}] = graphconncomp(sparse(TTT));
                 %[num2cell(group') path1(ind1)' path2(ind1)']
@@ -71,7 +70,7 @@
         end
 
          if (length(ind4)>1)
-                BPM_sim{4} = bpmsim(WPM.ind(ind4),WPM.ind(ind4),WPM.ind(ind4),WPM.ind(ind4),2,ssM{1});
+                BPM_sim{4} = bpmsim(WPM.ind(ind4),WPM.ind(ind4),WPM.ind(ind4),WPM.ind(ind4));
                 TTT=BPM_sim{4}>=0.25;
                 [noRD{4},group{4}] = graphconncomp(sparse(TTT));
                 %[num2cell(group') path1(ind1)' path2(ind1)']

diff --git a/scripts/run_computeFDR.sh b/scripts/run_computeFDR.sh
@@ -60,6 +60,7 @@ if [ -z "${projectDir}" ]; then projectDir=`pwd`; fi
 if [ -z "${minPath}" ]; then minPath=10; fi
 if [ -z "${pvalueCutoff}" ]; then pvalueCutoff=0.005; fi
 if [ -z "${bpmindFile}" ]; then bpmindFile=BPMind.mat; fi
+if [ -z "${samplePerms}" ]; then samplePerms=10; fi
 
 # check required inuts
 if [ -z "${ssmFile}" ]; then echo "Please specify ssmFile to proceed."; exit; fi