Merge pull request #1 from TCLamnidis/skipJK

Added option to skip jackknife error estimation

Author: TCLamnidis

RASCalculator.py

@@ -27,6 +27,7 @@
 parser.add_argument("-C", "--NrChroms", type=int, metavar="<INT>", default=22, required=False, help="The number of chromosomes in the dataset. [22]")
 parser.add_argument("-d", "--details", action='store_true', help="Print RAS calculations for each allele frequency, in addition to the total.")
 parser.add_argument("--f3FreqCutoff", type=float, metavar="<FREQ>", default=0.05, help="minimum minor allele frequency for f3 values. Defaults to 0.05. This is recommended to skip rare mutations and reduce the sensitivity of the F3 statistics on sequencing errors and DNA damage.")
+parser.add_argument("--skipJackknife", action='store_true', default=False, help="When provided, the jackknife error is not estimated, but instead set to 0. [False]")
 args = parser.parse_args()
 
 if args.minAF<1:
@@ -42,6 +43,8 @@
 NumBins=args.NrChroms
 minAF=args.minAF
 maxAF=args.maxAF
+skipJK=args.skipJackknife
+
 # RightIndex={} #Holds the INDICES of the Right pops
 # LeftsIndex={} #Holds the INDICES of the Left pops
 
@@ -160,7 +163,7 @@ def getTotalDerivedAC(afDict):
 for x in range(len(LeftPops)):
     for j in range(len(TestPops)):
         for i in range(minAF-1,maxAF+2):
-            thetaJ,sigma2 = ras.getJackknife(RAS[x][j][i], mj[x][j], blockSizes)
+            thetaJ,sigma2 = ras.getJackknife(RAS[x][j][i], mj[x][j], blockSizes, skipJK)
             ThetaJ[x][j][i] = thetaJ
             Sigma2[x][j][i] = sigma2
 

RASUtils.py

@@ -40,31 +40,35 @@ def __readHeader(self):
             self.sizes[popName] = popSize
         print("#Available populations in Input File and their respective sizes: ", self.sizes, file=self.output)
 
-def getJackknife(blockValues, totalObservations, blockSizes):
+def getJackknife(blockValues, totalObservations, blockSizes, skipJackknife):
     thetaminus=[0 for x in range(len(blockSizes))]
     sum1=0
     sum2=0
     jackknifeStdErr=0
     if sum(totalObservations)==0:
+        ## If no observations were made, the rare allele sharing rate is 0.
         thetahat=0
     else:
+        ## thetahat is normalised by number of observations
         thetahat=sum(blockValues)/sum(totalObservations)
-    
-    ## Normalise blockValues.
-    normalisedValues = [0.0 for c in range(len(blockSizes))]
+
     for c in range(len(blockSizes)):
-        if totalObservations[c] == 0:
-            continue
+        if totalObservations[c] == sum(totalObservations):
+            ## If all rare alleles are on a single chunk, the ras/site without that chunk is 0.
+            thetaminus[c]=0
         else:
-            normalisedValues[c] = blockValues[c]/totalObservations[c]
-    
-    for c in range(len(blockSizes)):
-        thetaminus[c]=( (sum(blockValues)-blockValues[c]) / (sum(totalObservations)-totalObservations[c]) )
+            ## thetaminus is normalised by number of observations
+            thetaminus[c]=( (sum(blockValues)-blockValues[c]) / (sum(totalObservations)-totalObservations[c]) )
+        
+        ## Jackknife estimator calculation
         sum1+=thetahat-thetaminus[c]
         sum2+=(blockSizes[c]*thetaminus[c])/sum(blockSizes)
     jackknifeEstimator=sum1+sum2
-    for c in range(len(blockSizes)):
-        hj=sum(blockSizes)/blockSizes[c]
-        pseudoval = (hj*thetahat)-((hj-1)*thetaminus[c])
-        jackknifeStdErr+=(1/len(blockSizes))*(((pseudoval-jackknifeEstimator)**2)/(hj-1))
+    
+    ## Standard error calculation
+    if not skipJackknife:
+        for c in range(len(blockSizes)):
+            hj=sum(blockSizes)/blockSizes[c]
+            pseudoval = (hj*thetahat)-((hj-1)*thetaminus[c])
+            jackknifeStdErr+=(1/len(blockSizes))*(((pseudoval-jackknifeEstimator)**2)/(hj-1))
     return (jackknifeEstimator,jackknifeStdErr)