Cauchy scorer

Author: Dan-RAI

python/PascalX/genescorer.py

@@ -1810,7 +1810,6 @@ def __init__(self,window=50000,varcutoff=0.99,MAF=0.05,genome=None,gpu=False):
             self._GENESYMB = genome._GENESYMB
             self._GENEIDtoSYMB = genome._GENEIDtoSYMB
             self._CHR = genome._CHR
-
             self._SKIPPED = genome._SKIPPED
 
 
@@ -1832,7 +1831,6 @@ def _getSNPs(self,cr,gene,REF,useAll=False):
         else:
             if gene in self._MAP:
                 P = set(REF[str(cr)][0].getSNPsPos(list(self._MAP[gene].keys())))
-                useAll = True
             else:
                 P = []
 
@@ -1872,7 +1870,6 @@ def _getSNPs_wAlleles(self,cr,gene,REF,useAll=False):
         else:
             if gene in self._MAP:
                 P = set(REF[str(cr)][0].getSNPsPos( list(self._MAP[gene].keys()) ))
-                useAll = True
             else:
                 P = []
 
@@ -1893,9 +1890,29 @@ def _getSNPs_wAlleles(self,cr,gene,REF,useAll=False):
         # Calc corr
         RID = list(filtered.keys())
 
-        return C,np.array(RID)
+        return np.array(RID)
 
 
+            
+    def _getChi2Sum_mapper(self,RIDs,gene):
+        #print([GWAS[x] for x in RIDs])
+        ps = np.zeros(len(RIDs))
+        for i in range(0,len(ps)):
+            #ps = chi2.ppf(1- np.array([GWAS[x] for x in RIDs]),1)
+            if RIDs[i] in self._MAP[gene] and self._MAP[gene][RIDs[i]][4] is not None:
+                ps[i] = self._MAP[gene][RIDs[i]][4]
+            else:
+                ps[i] = self._GWAS[RIDs[i]]
+                
+        return ps
+        
+    def _getChi2Sum(self,RIDs):
+        #print([GWAS[x] for x in RIDs])
+        ps = np.zeros(len(RIDs))
+        for i in range(0,len(ps)):
+            #ps = chi2.ppf(1- np.array([GWAS[x] for x in RIDs]),1)
+            ps[i] = self._GWAS[RIDs[i]]
+        return ps      
 
     def _scoremain(self,gene,unloadRef,label='',baroffset=0,nobar=False,lock=None):
 
@@ -1938,19 +1955,15 @@ def _scoremain(self,gene,unloadRef,label='',baroffset=0,nobar=False,lock=None):
                     R = self._getSNPs_wAlleles(cr,G[i],REF)
 
                 if len(R) > 0:
-                    
-                    #hpstats.cauchytest_100d()
-                      
-                    return
-                        
+                                    
                     # Score
                     if self._MAP is not None and self._joint == False:
 
                         S = self._getChi2Sum_mapper(R,G[i]) 
                     else:
                         S = self._getChi2Sum(R)
 
-                    #RES = self._scoreThread(self._calcAndFilterEV(C),S,G[i],method,mode,reqacc,intlimit)
+                    RES = [G[i],[hpstats.cauchytest_300d(S)]]
 
                     RESULT.append( [self._GENEIDtoSYMB[RES[0]],float(RES[1][0]),len(R)])
 

python/PascalX/hpstats.py

@@ -96,4 +96,38 @@ def onemin_norm_cdf_100d(x, m, s):
 
     return lib.onemin_normcdf_100d(x,m,s)
 
-    
+    
+def cauchytest_100d(x):
+    """
+    cauchy test statistics
+    
+    x : p-value input array
+    
+    """
+    
+    _L = np.ascontiguousarray(x, dtype='float64')
+
+    res = lib.cauchytest_100d(
+            ffi.cast("double *",_L.ctypes.data),
+            len(_L))
+  
+        
+    return max(res,1e-100)
+
+
+
+def cauchytest_300d(x):
+    """
+    cauchy test statistics
+    
+    x : p-value input array
+    
+    """
+    
+    _L = np.ascontiguousarray(x, dtype='float64')
+
+    res = lib.cauchytest_300d(
+            ffi.cast("double *",_L.ctypes.data),
+            len(_L))
+  
+    return max(res,1e-300)

python/hpstats.cpp

@@ -26,7 +26,7 @@
 #include <boost/math/distributions/normal.hpp>
 using namespace boost::multiprecision;
 using namespace boost::math;
-using namespace boost::constants;
+using namespace boost::math::constants;
 
 extern "C"
 double invchi2cdf_1mx(double x, int dof) {
@@ -118,5 +118,21 @@ double cauchytest_100d(double* x, int N) {
 
     cpp_bin_float_100 t = 0.5 - atan(sum)/MP;
 
+    return t.convert_to<double>();
+}
+
+extern "C"
+double cauchytest_300d(double* x, int N) {
+
+    number<cpp_bin_float<300>>  MP = pi<number<cpp_bin_float<300>>>();
+    
+    number<cpp_bin_float<300>>  sum = 0;
+    for(int i = 0; i < N; i++) {
+        number<cpp_bin_float<300>>  X = number<cpp_bin_float<300>> (x[i]);
+        sum += tan((0.5-X)*MP);
+    }
+    
+    number<cpp_bin_float<300>>  t = 0.5 - atan(sum)/MP;
+    
     return t.convert_to<double>();
 }

python/hpstats.hpp

@@ -30,4 +30,5 @@ extern double onemin_chi2cdf_100d(double x, int dof);
 extern double normcdf_100d(double x, double m, double s);
 extern double onemin_normcdf_100d(double x, double m, double s);
 
-extern double cauchytest_100d(double* x, int N);
+extern double cauchytest_100d(double* x, int N);
+extern double cauchytest_300d(double* x, int N);