Fixed a bug with assigning NAs and 1's to missing Vg's and balanced ehs.

R/ANEVA_DOT.R

@@ -70,7 +70,6 @@ ANEVA_DOT<-function(ASEdat, output_columns = c("refCount","altCount"), eh1 = "re
   return(output)
 }
 
-
 #' @param Eg_std Standard deviation in the log2 transformed total gene expression in a healthy
 #' population.
 #' @param eh1 Integer value for of expression of haplotype 1.
@@ -154,6 +153,27 @@ Binom_test_ctm_dbl<-function(X,N,p1,p2,log_BinCoeff,r0){
   return(p.val)
 }
 
+#' This function calculates Binomial coefficients for choosing 0:n out of n,
+#' in log scale. It is useful for precalculating the coefficients and
+#' providing them to "Binom_test_fast.R" or "pdf_binom_fast.R" when they
+#' are redone on the same N, over and over again. This function was
+#' designed by Pejman Mohammadi, 2018, Scripps Research, San Diego, CA.
+
+log_BinCoeffs<-function(n){
+  log_z<-numeric(n+1) #note first and last element will remain 0
+
+  lgamm<-lgamma(0:n+1) #precalculate gamma function
+
+  m1<-floor(n/2) #this is the middle. binom coeffs symmetric so we calculate only half
+  m2<-ceiling(n/2) #this is the "other" middle
+  x<- 1:m1
+  xc<- n-x
+  log_z[2:(m1+1)]<- (lgamm[n+1]-lgamm[x+1]-lgamm[xc+1])
+
+  log_z[n:(m2+1)]<-log_z[2:(m1+1)] #fill in the other side
+  return(log_z)
+}
+
 #' #################################################################################################
 #' Generate plot of ASE data showing significant outliers. (This Function not currently implemented)
 #'

R/ANEVAdot.R

@@ -58,10 +58,10 @@ ANEVAdot<-function(filepath, output_columns = c("refCount","altCount"), eh1 = "r
       output$p.val[i]<-NA
       next
     }
-    if (dat[i,eh1]==dat[i,eh2]){
+    if (dat[i,eh1]==dat[i,eh2] & (dat[i,eh1]+dat[i,eh2])>10){
       output$p.val[i]<-1
     }
-    else if ((dat[i,eh1]+dat[i,eh2])>8){
+    else if ((dat[i,eh1]+dat[i,eh2])>10){
       Hh<-max(dat[i,eh1],dat[i,eh2]) #higher expressed haplotype
       Lh<-min(dat[i,eh1],dat[i,eh2]) #lower expressed haplotype
       rad<-Eg_std[i]*4 #integration radius

R/Evaluate_ANEVAdot.R

@@ -23,14 +23,9 @@ test.filepath<-"testdata.txt"
 
 ####add column titles manually####
 
-<<<<<<< HEAD
 start_time <- Sys.time()
 test.result<-ANEVAdot(filepath = test.filepath, output_columns = c("eh1","eh2"),eh1 = "eh1", eh2 = "eh2",Eg_std=Sg)
 end_time <- Sys.time()
 runtime<-end_time - start_time
 qqplot(runif(10000),test.result$p.val)
-=======
-test.result<-ANEVAdot(filepath = test.filepath, output_columns = c("eh1","eh2"), eh1 = "eh1", eh2 = "eh2",Eg_std=Sg)
-qqplot(runif(1000),test.result$p.val)
->>>>>>> e615cc169111775374b560ddc35e4c3f964d07b2