eveilyeverafter / HMMancestry

R package using the Forward-Backward algorithm to infer genotypes, recombination hotspots, and gene conversion tracts from low-coverage next-generation sequence data.
1 stars 0 forks source link

Error in if (sums[a] != 2) { : missing value where TRUE/FALSE needed #7

Closed eveilyeverafter closed 9 years ago

eveilyeverafter commented 9 years ago

Example:

set.seed(1234567) # For reproducability l2 <- 10000 # number of loci to simulate rec2 <- 0.01 # recombination rate between each snp r2 <- recombine_index(rep(rec2, l2-1)) # recombination rate between each snp (vector form) p_a2 <- .999 # probability of correct sequencing assignment (1-sequence error rate) p2 <- make_parents(l2) # make the parent recomb_sim2 <- recombine(parents=p2, r.index=r2, mu.rate=0) # recombine parents sim_reads2 <- simulate_coverage(a=recomb_sim2, p_assign=p_a2, coverage=20) # simulate sequencing coverage fbres <- estimate_anc_fwd_back(snp_dat=sim_reads2, chr_name="I", p_assign=p_a2, p_trans=rec2) infer_tracks(data=fbres, tetrad=1) Error: could not find function "infer_tracks"

eveilyeverafter commented 9 years ago

The error traced back to the estimate_anc_fwd_back() function. In some cases there is a tie where the posterior probability of belonging to either ancestral allele is equal (0.5). I fixed this by randomly picking one of the states using the sample() function. There are other ways to go about this fix. For example, I could have called the state with the highest posterior probability on either side but that would create additional computation time and would start getting into problems if there were long regions along the chromosome that are equally likely to come from either parent (e.g., in extremely low coverage scenarios).