Skip to content

Commit

Permalink
Added libgsl to build process and added Bonf corr to approx
Browse files Browse the repository at this point in the history
  • Loading branch information
Andreas Wilm committed Dec 23, 2020
1 parent 3ce3d3a commit 7a6577c
Show file tree
Hide file tree
Showing 3 changed files with 22 additions and 15 deletions.
7 changes: 4 additions & 3 deletions configure.ac
Original file line number Diff line number Diff line change
Expand Up @@ -146,11 +146,12 @@ AX_PTHREAD([
AC_MSG_ERROR([No pthread support on this machine]))
#AX_PTHREAD()


# explicit libm check
# if any of these sit in unusual places use
# export LDFLAGS="-L$path" before calling configure
AC_CHECK_LIB(m, log,, AC_MSG_ERROR([Could not find libm]))
AC_CHECK_LIB(z, gzread,, AC_MSG_ERROR([Could not find libz]))

AC_CHECK_LIB([gslcblas],[cblas_dgemm])
AC_CHECK_LIB(gsl, gsl_cdf_poisson_P,, AC_MSG_WARN([libgsl not found. Not using fast approximation]))

# http://www.gnu.org/software/automake/manual/html_node/Python.html
AM_PATH_PYTHON([2.7])
Expand Down
3 changes: 2 additions & 1 deletion src/lofreq/Makefile.am
Original file line number Diff line number Diff line change
Expand Up @@ -41,4 +41,5 @@ LIBS_for_htslib = -lhts
endif

# note: order matters
lofreq_LDADD = $(LIBS_for_htslib) ../cdflib90/libcdf.a -l:libgsl.a -lcblas -lm
lofreq_LDADD = $(LIBS_for_htslib) ../cdflib90/libcdf.a
# -l:libgsl.a -lm
27 changes: 16 additions & 11 deletions src/lofreq/snpcaller.c
Original file line number Diff line number Diff line change
Expand Up @@ -1102,17 +1102,22 @@ snpcaller(long double *snp_pvalues,
goto free_and_exit;
}

// Only approximate if sufficient data available
if (num_err_probs > 100) {
long double mu = 0;
for (int i = 0; i < num_err_probs; ++i) {
mu += err_probs[i];
}
const long double poibin_approximation = 1 - gsl_cdf_poisson_P(max_noncons_count - 1, mu);
if (poibin_approximation > sig_level + 0.05) {
goto free_and_exit;
}
}
#ifdef HAVE_LIBGSL
#ifdef HAVE_LIBGSLCBLAS
// Only approximate if sufficient data available
if (num_err_probs > 100) {
long double mu = 0;
for (int i = 0; i < num_err_probs; ++i) {
mu += err_probs[i];
}
const long double poibin_approximation = 1 - gsl_cdf_poisson_P(max_noncons_count - 1, mu);
if (poibin_approximation * (double)bonf_factor > sig_level) {
goto free_and_exit;
}
}
#endif
#endif

probvec = poissbin(&pvalue, err_probs, num_err_probs,
max_noncons_count, bonf_factor, sig_level);

Expand Down

0 comments on commit 7a6577c

Please sign in to comment.