mmollina / MAPpoly

Compare 95236d0 ... +0 ... d9f81d1

No flags found

Use flags to group coverage reports by test type, project and/or folders.
Then setup custom commit statuses and notifications for each flag.

e.g., #unittest #integration

#production #enterprise

#frontend #backend

Learn more about Codecov Flags here.

Showing 2 of 7 files from the diff.

@@ -22,7 +22,10 @@
Loading
22 22
#'    time. However, it will require less memory. See examples.
23 23
#'    
24 24
#' @param est.type  Indicates whether to use the discrete ("disc") or the probabilistic ("prob") dosage scoring 
25 -
#'                  when estimating the two-point recombination fractions. 
25 +
#'                  when estimating the two-point recombination fractions. The option 
26 +
#'                  "sparse_grid" is for development purposes use only. 
27 +
#'                  It compares the likelihoods of a set of recombination fractions and returns the most likely.
28 +
#'                  Currently, it tests rf = c(0.0, 0.2, 0.5).
26 29
#'
27 30
#' @param verbose If \code{TRUE} (default), current progress is shown; if
28 31
#'     \code{FALSE}, no output is produced
@@ -84,7 +87,7 @@
Loading
84 87
#' @importFrom dplyr filter arrange
85 88
est_pairwise_rf <- function(input.seq, count.cache = NULL, ncpus = 1L,
86 89
                            mrk.pairs = NULL, n.batches = 1L,
87 -
                            est.type = c("disc","prob"),
90 +
                            est.type = c("disc","prob","sparse_grid"),
88 91
                            verbose = TRUE, memory.warning = TRUE, 
89 92
                            parallelization.type = c("PSOCK", "FORK"), 
90 93
                            tol = .Machine$double.eps^0.25)
@@ -130,12 +133,12 @@
Loading
130 133
  }
131 134
  est.type = match.arg(est.type)
132 135
  ## Checking for genotype probability 
133 -
  if(!exists('geno', where = get(input.seq$data.name, pos = 1)) & est.type != "disc"){
136 +
  if(!exists('geno', where = get(input.seq$data.name, pos = 1)) & est.type == "prob"){
134 137
    warning("There is no probabilistic dosage scoring in the dataset. Using est.type = 'disc'")
135 138
    est.type <- "disc"
136 139
  }
137 140
  ## get genotypes
138 -
  if(est.type  ==  "disc"){
141 +
  if(est.type  ==  "disc" || est.type == "sparse_grid"){
139 142
    geno <- as.matrix(get(input.seq$data.name, pos = 1)$geno.dose)
140 143
  } else {
141 144
    d1 <- get(input.seq$data.name, pos = 1)$geno 
@@ -169,6 +172,8 @@
Loading
169 172
        parallel::clusterExport(cl, "paralell_pairwise_discrete")
170 173
      if(est.type  ==  "prob")
171 174
        parallel::clusterExport(cl, "paralell_pairwise_probability")
175 +
      if(est.type  ==  "sparse_grid")
176 +
        parallel::clusterExport(cl, "paralell_pairwise_given_sparse_grid")
172 177
      on.exit(parallel::stopCluster(cl))
173 178
      if(est.type  ==  "disc"){
174 179
        res <- parallel::parLapply(cl,
@@ -192,6 +197,17 @@
Loading
192 197
                                   count.cache = count.cache,
193 198
                                   tol = tol)
194 199
      } 
200 +
      else if(est.type  ==  "sparse_grid"){
201 +
        res <- parallel::parLapply(cl,
202 +
                                   input.list,
203 +
                                   paralell_pairwise_given_sparse_grid,
204 +
                                   input.seq = input.seq,
205 +
                                   geno = geno,
206 +
                                   dP = get(input.seq$data.name)$dosage.p1,
207 +
                                   dQ = get(input.seq$data.name)$dosage.p2,
208 +
                                   count.cache = count.cache)
209 +
      } 
210 +
      
195 211
      end <- proc.time()
196 212
      if (verbose) {
197 213
        cat("INFO: Done with",
@@ -228,6 +244,14 @@
Loading
228 244
                      dQ = get(input.seq$data.name)$dosage.p2,
229 245
                      count.cache = count.cache,
230 246
                      tol = tol)
247 +
      } else if(est.type  ==  "sparse_grid"){
248 +
        res <- lapply(input.list,
249 +
                      paralell_pairwise_given_sparse_grid,
250 +
                      input.seq = input.seq,
251 +
                      geno = geno,
252 +
                      dP = get(input.seq$data.name)$dosage.p1,
253 +
                      dQ = get(input.seq$data.name)$dosage.p2,
254 +
                      count.cache = count.cache)
231 255
      } 
232 256
    }
233 257
    res <- unlist(res,
@@ -362,13 +386,28 @@
Loading
362 386
  return(lapply(res, format_rf))
363 387
}
364 388
365 -
366 -
367 -
368 -
369 -
370 -
371 -
389 +
#' Wrapper function to discrete-based pairwise likelihood computation given rf = 1e-5
390 +
#'
391 +
#' @param void internal function to be documented
392 +
#' @keywords internal
393 +
#' @export
394 +
paralell_pairwise_given_sparse_grid <- function(mrk.pairs,
395 +
                                       input.seq,
396 +
                                       geno,
397 +
                                       dP,
398 +
                                       dQ,
399 +
                                       count.cache)
400 +
{
401 +
  res <- .Call("pairwise_given_0rf",
402 +
               input.seq$ploidy,
403 +
               as.matrix(mrk.pairs),
404 +
               as.matrix(geno),
405 +
               as.vector(dP),
406 +
               as.vector(dQ),
407 +
               count.cache$cond,
408 +
               PACKAGE = "mappoly")
409 +
  return(lapply(res, format_rf))
410 +
}
372 411
373 412
#' Format results from pairwise two-point estimation in C++
374 413
#'

@@ -499,4 +499,78 @@
Loading
499 499
  delete [] G;
500 500
  return(out);
501 501
}
502 +
RcppExport SEXP pairwise_given_0rf(SEXP m_R,
503 +
                                   SEXP mrk_pairs_R,
504 +
                                   SEXP geno_R,
505 +
                                   SEXP dP_R,
506 +
                                   SEXP dQ_R,
507 +
                                   SEXP count_cache_R)
508 +
{
509 +
  Rcpp::NumericMatrix mrk_pairs = Rcpp::as<Rcpp::NumericMatrix>(mrk_pairs_R);
510 +
  Rcpp::NumericMatrix geno = Rcpp::as<Rcpp::NumericMatrix>(geno_R);
511 +
  Rcpp::NumericVector dP = Rcpp::as<Rcpp::NumericVector>(dP_R);
512 +
  Rcpp::NumericVector dQ = Rcpp::as<Rcpp::NumericVector>(dQ_R);
513 +
  Rcpp::List count_cache = Rcpp::as<Rcpp::List>(count_cache_R);
514 +
  Rcpp::NumericVector d_pair(4);
515 +
  Rcpp::List out(mrk_pairs.ncol());
516 +
  int m = Rcpp::as<int>(m_R);
517 +
  int n_ind = geno.ncol();
518 +
  for(int k=0; k < mrk_pairs.ncol(); k++)
519 +
  {
520 +
    //Rcpp::Rcout << mrk_pairs(0,k) << " - " << mrk_pairs(1,k) <<  std::endl;
521 +
    int id = (m+1)*(m+1)*(m+1)*dQ[mrk_pairs(1,k)]+(m+1)*(m+1)*dQ[mrk_pairs(0,k)]+(m+1)*dP[mrk_pairs(1,k)]+dP[mrk_pairs(0,k)]+1;
522 +
    //Rcpp::Rcout << "id: " << id <<  std::endl;
523 +
    Rcpp::List temp_list = count_cache[(id-1)];
524 +
    //Rcpp::List temp_list = count_cache[k];
525 +
    if(temp_list.size() > 1)
526 +
    {
527 +
      NumericVector gen_1 = geno( mrk_pairs(0,k), _);
528 +
      NumericVector gen_2 = geno( mrk_pairs(1,k), _);
529 +
      Rcpp::NumericMatrix res(3, temp_list.size());
530 +
      Rcpp::CharacterVector zn = temp_list.attr( "names" ) ;
531 +
      colnames(res)=zn;
532 +
      for(int i=0; i < temp_list.size(); i++)
533 +
      {
534 +
        Rcpp::NumericMatrix count_mat = temp_list[i] ;
535 +
        Rcpp::List dimnames = count_mat.attr( "dimnames" ) ;
536 +
        Rcpp::CharacterVector z = dimnames[0];
537 +
        Rcpp::NumericVector dk(z.size()), dk1(z.size());
538 +
        std::string delimiter = " ";
539 +
        for(int j=0; j < z.size(); j++)
540 +
        {
541 +
          std::string lnames = Rcpp::as<std::string>(z(j));
542 +
          dk(j) = std::stoi(lnames.substr(0,lnames.find(delimiter)));
543 +
          dk1(j) = std::stoi(lnames.substr(lnames.find(delimiter)+1, lnames.length()));
544 +
        }
545 +
        vector<double> x{0.0, 0.2, 0.5};
546 +
        double lltemp, ll, curx;
547 +
        ll = twopt_likelihood_dosage(x[0], m, n_ind, dP[mrk_pairs(0,k)], dQ[mrk_pairs(0,k)], dk, dk1, gen_1, gen_2, count_mat);
548 +
        curx = x[0];
549 +
        for(int j=1; j < x.size(); j++)
550 +
        {
551 +
          lltemp = twopt_likelihood_dosage(x[j], m, n_ind, dP[mrk_pairs(0,k)], dQ[mrk_pairs(0,k)], dk, dk1, gen_1, gen_2, count_mat);
552 +
          if(ll - lltemp < 0.0) break;
553 +
          else{
554 +
            curx = x[j];
555 +
            ll = lltemp;
556 +
          }
557 +
        }
558 +
          res(0,i) = curx;
559 +
          res(1,i) = ll;
560 +
          res(2,i) = ll;
561 +
      }
562 +
      out(k)=res;
563 +
    }
564 +
    else
565 +
    {
566 +
      Rcpp::NumericVector d_out(4);
567 +
      d_out(0)=dP[mrk_pairs(0,k)];
568 +
      d_out(1)=dP[mrk_pairs(1,k)];
569 +
      d_out(2)=dQ[mrk_pairs(0,k)];
570 +
      d_out(3)=dQ[mrk_pairs(1,k)];
571 +
      out(k)=d_out;
572 +
    }
573 +
  }
574 +
  return(out);
575 +
}
502 576
//end of file two_pts_est.cpp

Learn more Showing 1 files with coverage changes found.

Changes in R/filters.R
-1
+1
Loading file...
Files Coverage
R -0.31% 69.07%
src -2.18% 74.53%
Project Totals (59 files) 70.32%
Loading