mmollina / MAPpoly

Compare d8ea82f ... +10 ... b1cb497

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 47 of 173 files from the diff.
Other files ignored by Codecov
man/reest_rf.Rd has changed.
man/est_rf_hmm.Rd has changed.
data/hexafake.rda was deleted.
NAMESPACE has changed.
man/merge_maps.Rd has changed.
man/get_LOD.Rd has changed.
man/plot_GIC.Rd has changed.
R/sysdata.rda has changed.
man/draw_cross.Rd has changed.
man/read_geno.Rd has changed.
man/add_marker.Rd has changed.
man/imf_m.Rd has changed.
man/mf_m.Rd has changed.
man/imf_h.Rd has changed.
man/perm_tot.Rd has changed.
man/print_mrk.Rd has changed.
man/get_counts.Rd has changed.
man/rev_map.Rd has changed.
man/imf_k.Rd has changed.
man/mf_k.Rd has changed.
R/data.R has changed.
man/format_rf.Rd has changed.
man/read_vcf.Rd has changed.
man/get_submap.Rd has changed.
README.md has changed.
man/create_map.Rd has changed.
man/mf_h.Rd has changed.
man/hexafake.Rd has changed.
man/select_rf.Rd has changed.
man/get_ij.Rd has changed.
man/elim_equiv.Rd has changed.
man/get_w_m.Rd has changed.
man/perm_pars.Rd has changed.
DESCRIPTION has changed.

@@ -30,7 +30,7 @@
Loading
30 30
#' 
31 31
#' @examples
32 32
#'  \donttest{
33 -
#'      probs.error<-calc_genoprob_error(input.map = solcap.err.map[[1]],
33 +
#'      probs.error <- calc_genoprob_error(input.map = solcap.err.map[[1]],
34 34
#'                                 error = 0.05,
35 35
#'                                 verbose = TRUE)
36 36
#'  }
@@ -45,7 +45,7 @@
Loading
45 45
#'     \doi{10.1534/g3.119.400378} 
46 46
#'
47 47
#' @export calc_genoprob_error
48 -
calc_genoprob_error<-function(input.map,  step = 0, phase.config = "best", error = 0.01, 
48 +
calc_genoprob_error <- function(input.map,  step = 0, phase.config = "best", error = 0.01, 
49 49
                              th.prob = 0.95, restricted = TRUE, verbose = TRUE)
50 50
{
51 51
  if (!inherits(input.map, "mappoly.map")) {
@@ -56,75 +56,75 @@
Loading
56 56
  }
57 57
  ## choosing the linkage phase configuration
58 58
  LOD.conf <- get_LOD(input.map, sorted = FALSE)
59 -
  if(phase.config == "best") {
59 +
  if(phase.config  ==  "best") {
60 60
    i.lpc <- which.min(LOD.conf)
61 61
  } else if (phase.config > length(LOD.conf)) {
62 62
    stop("invalid linkage phase configuration")
63 63
  } else i.lpc <- phase.config
64 -
  m<-input.map$info$m
65 -
  n.ind<-get(input.map$info$data.name, pos=1)$n.ind
66 -
  Dtemp<-get(input.map$info$data.name, pos=1)$geno.dose[input.map$maps[[1]]$seq.num,]
64 +
  ploidy <- input.map$info$ploidy
65 +
  n.ind <- get(input.map$info$data.name, pos = 1)$n.ind
66 +
  Dtemp <- get(input.map$info$data.name, pos = 1)$geno.dose[input.map$maps[[1]]$seq.num,]
67 67
  map.pseudo <- create_map(input.map, step, phase.config = i.lpc)
68 -
  mrknames<-names(map.pseudo)
69 -
  n.mrk<-length(map.pseudo)
70 -
  indnames<-colnames(Dtemp)
71 -
  D<-matrix(m+1, nrow = length(map.pseudo), ncol = ncol(Dtemp), 
72 -
            dimnames = list(mrknames, indnames))
68 +
  mrk.names <- names(map.pseudo)
69 +
  n.mrk <- length(map.pseudo)
70 +
  ind.names <- colnames(Dtemp)
71 +
  D <- matrix(ploidy+1, nrow = length(map.pseudo), ncol = ncol(Dtemp), 
72 +
            dimnames = list(mrk.names, ind.names))
73 73
  D[rownames(Dtemp), ] <- as.matrix(Dtemp)
74 -
  dptemp <- get(input.map$info$data.name)$dosage.p[input.map$maps[[i.lpc]]$seq.num]
75 -
  dqtemp <- get(input.map$info$data.name)$dosage.q[input.map$maps[[i.lpc]]$seq.num]
76 -
  dq<-dp<-rep(m/2, length(mrknames))
77 -
  names(dp)<-names(dq)<-mrknames
78 -
  dp[names(dptemp)]<-dptemp
79 -
  dq[names(dqtemp)]<-dqtemp
74 +
  dptemp <- get(input.map$info$data.name)$dosage.p1[input.map$maps[[i.lpc]]$seq.num]
75 +
  dqtemp <- get(input.map$info$data.name)$dosage.p2[input.map$maps[[i.lpc]]$seq.num]
76 +
  dq <- dp <- rep(ploidy/2, length(mrk.names))
77 +
  names(dp) <- names(dq) <- mrk.names
78 +
  dp[names(dptemp)] <- dptemp
79 +
  dq[names(dqtemp)] <- dqtemp
80 80
 
81 81
  phP <- phQ <- vector("list", n.mrk)
82 82
  for(i in 1:length(phP)){
83 -
    phP[[i]] <- phQ[[i]] <- c(0:(m/2 - 1))
83 +
    phP[[i]] <- phQ[[i]] <- c(0:(ploidy/2 - 1))
84 84
  }
85 -
  names(phP) <- names(phQ) <- mrknames
85 +
  names(phP) <- names(phQ) <- mrk.names
86 86
  phP[rownames(Dtemp)] <- input.map$maps[[i.lpc]]$seq.ph$P
87 87
  phQ[rownames(Dtemp)] <- input.map$maps[[i.lpc]]$seq.ph$Q
88 88
89 89
  ## Including error  
90 -
  gen<-vector("list", length(indnames))
91 -
  names(gen)<-indnames
92 -
  d.pq<-data.frame(dp = dp, 
90 +
  gen <- vector("list", length(ind.names))
91 +
  names(gen) <- ind.names
92 +
  d.pq <- data.frame(dp = dp, 
93 93
                   dq = dq)
94 -
  d.pq$mrk<-rownames(d.pq)
94 +
  d.pq$mrk <- rownames(d.pq)
95 95
  for(i in names(gen))
96 96
  {
97 -
    a<-matrix(0, nrow(D), input.map$info$m+1, dimnames = list(mrknames, 0:input.map$info$m))
97 +
    a <- matrix(0, nrow(D), input.map$info$ploidy+1, dimnames = list(mrk.names, 0:input.map$info$ploidy))
98 98
    for(j in rownames(a)){
99 -
      if(D[j,i] == input.map$info$m+1){
100 -
        a[j,]<-segreg_poly(m = input.map$info$m, dP = dp[j], dQ = dq[j])
99 +
      if(D[j,i]  ==  input.map$info$ploidy+1){
100 +
        a[j,] <- segreg_poly(ploidy = input.map$info$ploidy, dP = dp[j], dQ = dq[j])
101 101
      } else {
102 -
        a[j,D[j,i]+1]<-1          
102 +
        a[j,D[j,i]+1] <- 1          
103 103
      }
104 104
    }
105 -
    a<-as.data.frame(a)
106 -
    a$mrk<-rownames(a)
107 -
    a.temp<-t(merge(a, d.pq, sort = FALSE)[,-c(1)])
105 +
    a <- as.data.frame(a)
106 +
    a$mrk <- rownames(a)
107 +
    a.temp <- t(merge(a, d.pq, sort = FALSE)[,-c(1)])
108 108
    if(!is.null(error))
109 -
      a.temp<-apply(a.temp, 2, genotyping_global_error, 
110 -
                    m = input.map$info$m, error=error, 
109 +
      a.temp <- apply(a.temp, 2, genotyping_global_error, 
110 +
                    ploidy = input.map$info$ploidy, error = error, 
111 111
                    th.prob = th.prob, 
112 112
                    restricted = restricted)
113 113
    else
114 -
      a.temp <- a.temp[1:(input.map$info$m+1), ]
115 -
    colnames(a.temp)<-a[,1]
116 -
    gen[[i]]<-a.temp
114 +
      a.temp <- a.temp[1:(input.map$info$ploidy+1), ]
115 +
    colnames(a.temp) <- a[,1]
116 +
    gen[[i]] <- a.temp
117 117
  }
118 118
  g = as.double(unlist(gen))
119 119
  p = as.numeric(unlist(phP))
120 120
  dp = as.numeric(cumsum(c(0, sapply(phP, function(x) sum(length(x))))))
121 121
  q = as.numeric(unlist(phQ))
122 122
  dq = as.numeric(cumsum(c(0, sapply(phQ, function(x) sum(length(x))))))
123 123
  rf = as.double(mf_h(diff(map.pseudo)))
124 -
  res.temp <-
124 +
  res.temp  <- 
125 125
    .Call(
126 126
      "calc_genoprob_prior",
127 -
      as.numeric(m),
127 +
      as.numeric(ploidy),
128 128
      as.numeric(n.mrk),
129 129
      as.numeric(n.ind),
130 130
      as.numeric(p),
@@ -133,15 +133,15 @@
Loading
133 133
      as.numeric(dq),
134 134
      as.double(g),
135 135
      as.double(rf),
136 -
      as.numeric(rep(0, choose(m, m/2)^2 * n.mrk * n.ind)),
136 +
      as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)),
137 137
      as.double(0),
138 138
      as.numeric(verbose),
139 139
      PACKAGE = "mappoly"
140 140
    )
141 141
  if (verbose) cat("\n")
142 -
  dim(res.temp[[1]])<-c(choose(m,m/2)^2,n.mrk,n.ind)
143 -
  dimnames(res.temp[[1]])<-list(kronecker(apply(combn(letters[1:m],m/2),2, paste, collapse=""),
144 -
                                          apply(combn(letters[(m+1):(2*m)],m/2),2, paste, collapse=""), paste, sep=":"),
145 -
                                mrknames, indnames)
146 -
  structure(list(probs = res.temp[[1]], map = map.pseudo), class="mappoly.genoprob")
142 +
  dim(res.temp[[1]]) <- c(choose(ploidy,ploidy/2)^2,n.mrk,n.ind)
143 +
  dimnames(res.temp[[1]]) <- list(kronecker(apply(combn(letters[1:ploidy],ploidy/2),2, paste, collapse = ""),
144 +
                                          apply(combn(letters[(ploidy+1):(2*ploidy)],ploidy/2),2, paste, collapse = ""), paste, sep = ":"),
145 +
                                mrk.names, ind.names)
146 +
  structure(list(probs = res.temp[[1]], map = map.pseudo), class = "mappoly.genoprob")
147 147
}

@@ -4,7 +4,7 @@
Loading
4 4
#' and the dosage of the locus in both parents. It does not consider
5 5
#' double reduction.
6 6
#'
7 -
#' @param m the ploidy level
7 +
#' @param ploidy the ploidy level
8 8
#'
9 9
#' @param dP the dosage in parent P
10 10
#'
@@ -16,8 +16,8 @@
Loading
16 16
#' @examples
17 17
#' # autohexaploid with two and three doses in parents P and Q,
18 18
#' # respectively
19 -
#' seg<-segreg_poly(m=6, dP=2, dQ=3)
20 -
#' barplot(seg, las=2)
19 +
#' seg <- segreg_poly(ploidy = 6, dP = 2, dQ = 3)
20 +
#' barplot(seg, las = 2)
21 21
#'
22 22
#' @author Marcelo Mollinari, \email{mmollin@ncsu.edu}
23 23
#'
@@ -37,21 +37,21 @@
Loading
37 37
#' @importFrom stats dhyper
38 38
#' @export segreg_poly
39 39
#'
40 -
segreg_poly <- function(m, dP, dQ) {
41 -
    if (m%%2 != 0)
40 +
segreg_poly <- function(ploidy, dP, dQ) {
41 +
    if (ploidy%%2 != 0)
42 42
        stop("m must be an even number")
43 -
    p.dose <- numeric((m + 1))
44 -
    p.names <- character((m + 1))
45 -
    seg.p1 <- dhyper(x = c(0:(m + 1)), m = dP, n = (m - dP), k = m/2)
46 -
    seg.p2 <- dhyper(x = c(0:(m + 1)), m = dQ, n = (m - dQ), k = m/2)
43 +
    p.dose <- numeric((ploidy + 1))
44 +
    p.names <- character((ploidy + 1))
45 +
    seg.p1 <- dhyper(x = c(0:(ploidy + 1)), m = dP, n = (ploidy - dP), k = ploidy/2)
46 +
    seg.p2 <- dhyper(x = c(0:(ploidy + 1)), m = dQ, n = (ploidy - dQ), k = ploidy/2)
47 47
    M <- tcrossprod(seg.p1, seg.p2)
48 48
    for (i in 1:nrow(M)) {
49 49
        for (j in 1:ncol(M)) {
50 50
            p.dose[i + j - 1] <- p.dose[i + j - 1] + M[i, j]
51 51
        }
52 52
    }
53 53
    p.dose <- p.dose[!is.na(p.dose)]
54 -
    for (i in 0:m) p.names[i + 1] <- paste(paste(rep("A", i), collapse = ""), paste(rep("a", (m - i)), collapse = ""), sep = "")
54 +
    for (i in 0:ploidy) p.names[i + 1] <- paste(paste(rep("A", i), collapse = ""), paste(rep("a", (ploidy - i)), collapse = ""), sep = "")
55 55
    names(p.dose) <- p.names
56 56
    return(p.dose)
57 57
}

@@ -50,17 +50,17 @@
Loading
50 50
                                      filter.non.conforming = TRUE,
51 51
                                      verbose = TRUE){
52 52
  input.type <- match.arg(input.type)
53 -
  if(input.type == "discrete"){
53 +
  if(input.type  ==  "discrete"){
54 54
    geno.dose <- input.data[,-match(c(parent1, parent2), colnames(input.data)), drop = FALSE]
55 -
    mappoly.data <- structure(list(m = ploidy,
55 +
    mappoly.data <- structure(list(ploidy = ploidy,
56 56
                                   n.ind = ncol(geno.dose),
57 57
                                   n.mrk = nrow(geno.dose),
58 58
                                   ind.names = colnames(geno.dose),
59 59
                                   mrk.names = rownames(geno.dose),
60 -
                                   dosage.p = input.data[,parent1],
61 -
                                   dosage.q = input.data[,parent2],
62 -
                                   sequence = NA,
63 -
                                   sequence.pos = NA,
60 +
                                   dosage.p1 = input.data[,parent1],
61 +
                                   dosage.p2 = input.data[,parent2],
62 +
                                   chrom = NA,
63 +
                                   genome.pos = NA,
64 64
                                   seq.ref = NULL,
65 65
                                   seq.alt = NULL,
66 66
                                   all.mrk.depth = NULL,
@@ -76,7 +76,7 @@
Loading
76 76
    if(is.null(pardose)) 
77 77
      stop("provide parental dosage.")
78 78
    rownames(pardose) <- pardose$MarkerName
79 -
    dat<-input.data[,c("MarkerName", "SampleName",paste0("P", 0:ploidy))]
79 +
    dat <- input.data[,c("MarkerName", "SampleName",paste0("P", 0:ploidy))]
80 80
    p1 <- unique(sapply(parent1, function(x) unique(grep(pattern = x, dat[,"SampleName"], value = TRUE))))
81 81
    p2 <- unique(sapply(parent2, function(x) unique(grep(pattern = x, dat[,"SampleName"], value = TRUE))))
82 82
    if(is.null(offspring)){
@@ -95,17 +95,17 @@
Loading
95 95
    ## get individual names ------------------
96 96
    ind.names <- offspring
97 97
    ## get dosage in parent P ----------------
98 -
    dosage.p <- as.integer(pardose[mrk.names,"parent1"])
99 -
    names(dosage.p)<-mrk.names
98 +
    dosage.p1 <- as.integer(pardose[mrk.names,"parent1"])
99 +
    names(dosage.p1) <- mrk.names
100 100
    ## get dosage in parent Q ----------------
101 -
    dosage.q <- as.integer(pardose[mrk.names,"parent2"])
102 -
    names(dosage.q)<-mrk.names
101 +
    dosage.p2 <- as.integer(pardose[mrk.names,"parent2"])
102 +
    names(dosage.p2) <- mrk.names
103 103
    ## monomorphic markers
104 -
    dp<-abs(abs(dosage.p-(ploidy/2))-(ploidy/2))
105 -
    dq<-abs(abs(dosage.q-(ploidy/2))-(ploidy/2))
106 -
    mrk.names<-names(which(dp+dq!=0))
107 -
    dosage.p <- dosage.p[mrk.names]
108 -
    dosage.q <- dosage.q[mrk.names]
104 +
    dp <- abs(abs(dosage.p1-(ploidy/2))-(ploidy/2))
105 +
    dq <- abs(abs(dosage.p2-(ploidy/2))-(ploidy/2))
106 +
    mrk.names <- names(which(dp+dq != 0))
107 +
    dosage.p1 <- dosage.p1[mrk.names]
108 +
    dosage.p2 <- dosage.p2[mrk.names]
109 109
    nphen <- 0
110 110
    phen <- NULL
111 111
    if (verbose){
@@ -126,15 +126,15 @@
Loading
126 126
    colnames(geno) <- c("mrk", "ind", as.character(0:ploidy))
127 127
    ind.names <- unique(geno$ind)
128 128
    mrk.names <- unique(geno$mrk)
129 -
    dosage.p <- dosage.p[mrk.names]
130 -
    dosage.q <- dosage.q[mrk.names]
129 +
    dosage.p1 <- dosage.p1[mrk.names]
130 +
    dosage.p2 <- dosage.p2[mrk.names]
131 131
    
132 132
    ## transforming na's in expected genotypes using Mendelian segregation
133 133
    i.na <- which(apply(geno, 1, function(x) any(is.na(x))))
134 134
    if (length(i.na) > 0) {
135 135
      m.na <- match(geno[i.na, 1], mrk.names)
136 -
      dp.na <- dosage.p[m.na]
137 -
      dq.na <- dosage.q[m.na]
136 +
      dp.na <- dosage.p1[m.na]
137 +
      dq.na <- dosage.p2[m.na]
138 138
      for (i in 1:length(m.na)) geno[i.na[i], -c(1, 2)] <- segreg_poly(ploidy, dp.na[i], dq.na[i])
139 139
    }
140 140
    ## dosage info
@@ -155,15 +155,15 @@
Loading
155 155
    }
156 156
    ## returning the 'mappoly.data' object
157 157
    if (verbose) cat("\n    Done with reading.\n")
158 -
    mappoly.data <- structure(list(m = ploidy,
158 +
    mappoly.data <- structure(list(ploidy = ploidy,
159 159
                                   n.ind = n.ind,
160 160
                                   n.mrk = length(mrk.names),
161 161
                                   ind.names = ind.names,
162 162
                                   mrk.names = mrk.names,
163 -
                                   dosage.p = dosage.p,
164 -
                                   dosage.q = dosage.q,
165 -
                                   sequence = rep(NA, length(mrk.names)),
166 -
                                   sequence.pos = rep(NA, length(mrk.names)),
163 +
                                   dosage.p1 = dosage.p1,
164 +
                                   dosage.p2 = dosage.p2,
165 +
                                   chrom = rep(NA, length(mrk.names)),
166 +
                                   genome.pos = rep(NA, length(mrk.names)),
167 167
                                   seq.ref = NULL,
168 168
                                   seq.alt = NULL,
169 169
                                   all.mrk.depth = NULL,
@@ -178,16 +178,16 @@
Loading
178 178
                              class = "mappoly.data")
179 179
  }
180 180
  if(filter.non.conforming){
181 -
    mappoly.data<-filter_non_conforming_classes(mappoly.data)
181 +
    mappoly.data <- filter_non_conforming_classes(mappoly.data)
182 182
    Ds <- array(NA, dim = c(ploidy+1, ploidy+1, ploidy+1))
183 183
    for(i in 0:ploidy)
184 184
      for(j in 0:ploidy)
185 -
        Ds[i+1,j+1,] <- segreg_poly(m = ploidy, dP = i, dQ = j)
186 -
    Dpop<-cbind(mappoly.data$dosage.p, mappoly.data$dosage.q)
187 -
    M<-t(apply(Dpop, 1, function(x) Ds[x[1]+1, x[2]+1,]))
188 -
    dimnames(M)<-list(mappoly.data$mrk.names, c(0:ploidy))
189 -
    M<-cbind(M, mappoly.data$geno.dose)
190 -
    mappoly.data$chisq.pval<-apply(M, 1, mrk_chisq_test, m = ploidy)
185 +
        Ds[i+1,j+1,] <- segreg_poly(ploidy = ploidy, dP = i, dQ = j)
186 +
    Dpop <- cbind(mappoly.data$dosage.p1, mappoly.data$dosage.p2)
187 +
    M <- t(apply(Dpop, 1, function(x) Ds[x[1]+1, x[2]+1,]))
188 +
    dimnames(M) <- list(mappoly.data$mrk.names, c(0:ploidy))
189 +
    M <- cbind(M, mappoly.data$geno.dose)
190 +
    mappoly.data$chisq.pval <- apply(M, 1, mrk_chisq_test, ploidy = ploidy)
191 191
  }
192 192
  mappoly.data
193 193
}
@@ -229,27 +229,27 @@
Loading
229 229
  }
230 230
  X <- maplist[[1]]
231 231
  if(is.null(ploidy))
232 -
    m <- (ncol(X)-2)/2
232 +
    ploidy <- (ncol(X)-2)/2
233 233
  MAPs <- vector("list", length(maplist))
234 234
  for(i in 1:length(MAPs)){
235 235
    X <- maplist[[i]]
236 236
    seq.num <- match(X$marker, mappoly.data$mrk.names)
237 237
    seq.rf <- mf_h(diff(X$position)) ## Using haldane
238 238
    seq.rf[seq.rf <= 1e-05] <- 1e-4
239 -
    P = ph_matrix_to_list(X[,3:(m+2)])
240 -
    Q = ph_matrix_to_list(X[,3:(m+2) + m])
239 +
    P = ph_matrix_to_list(X[,3:(ploidy+2)])
240 +
    Q = ph_matrix_to_list(X[,3:(ploidy+2) + ploidy])
241 241
    names(P) <- names(Q) <- seq.num
242 242
    seq.ph <- list(P = P, Q = Q)
243 243
    maps <- vector("list", 1)
244 244
    maps[[1]] <- list(seq.num = seq.num, seq.rf = seq.rf, seq.ph = seq.ph, loglike = 0)
245 -
    MAPs[[i]] <- structure(list(info = list(m = (ncol(X)-2)/2,
245 +
    MAPs[[i]] <- structure(list(info = list(ploidy = (ncol(X)-2)/2,
246 246
                                            n.mrk = nrow(X),
247 247
                                            seq.num = seq.num,
248 248
                                            mrk.names = as.character(X$marker),
249 -
                                            seq.dose.p = mappoly.data$dosage.p[as.character(X$marker)],
250 -
                                            seq.dose.q = mappoly.data$dosage.q[as.character(X$marker)],
251 -
                                            sequence = rep(i, nrow(X)),
252 -
                                            sequence.pos = NULL,
249 +
                                            seq.dose.p1 = mappoly.data$dosage.p1[as.character(X$marker)],
250 +
                                            seq.dose.p2 = mappoly.data$dosage.p2[as.character(X$marker)],
251 +
                                            chrom = rep(i, nrow(X)),
252 +
                                            genome.pos = NULL,
253 253
                                            seq.ref = NULL,
254 254
                                            seq.alt = NULL,
255 255
                                            chisq.pval = mappoly.data$chisq.pval[as.character(X$marker)],

@@ -9,7 +9,7 @@
Loading
9 9
#'     s.go <- make_seq_mappoly(id)
10 10
#'     ## Using the 5 contiguous markers
11 11
#'     seq5 <- make_seq_mappoly(hexafake, s.go$seq.mrk.names[6:10])
12 -
#'     twopt<-est_pairwise_rf(seq5)
12 +
#'     twopt <- est_pairwise_rf(seq5)
13 13
#'     l5 <- ls_linkage_phases(input.seq = seq5, thres = 2, twopt = twopt)
14 14
#'     plot(l5)
15 15
#'     
@@ -19,18 +19,18 @@
Loading
19 19
#'       maps1[[i]] <- est_rf_hmm_single(seq5, l5$config.to.test[[i]], 
20 20
#'                                       tol = 10e-3,
21 21
#'                                       high.prec = FALSE)
22 -
#'    (best<-which.max(sapply(maps1, function(x) x$loglike)))
23 -
#'    dist1<-round(cumsum(c(0, imf_h(maps1[[best]]$seq.rf))),2)
22 +
#'    (best <- which.max(sapply(maps1, function(x) x$loglike)))
23 +
#'    dist1 <- round(cumsum(c(0, imf_h(maps1[[best]]$seq.rf))),2)
24 24
#'    
25 25
#'    ## Same thing using automatic search
26 -
#'    maps2<-est_rf_hmm(input.seq = seq5, twopt = twopt, thres = 2, 
26 +
#'    maps2 <- est_rf_hmm(input.seq = seq5, twopt = twopt, thres = 2, 
27 27
#'                      verbose = TRUE, tol = 10e-3, high.prec = FALSE)
28 28
#'    plot(maps2)
29 29
#'    dist1
30 30
#'  }
31 31
#' @author Marcelo Mollinari, \email{mmollin@ncsu.edu}
32 32
#' @export est_rf_hmm_single
33 -
est_rf_hmm_single<-function(input.seq,
33 +
est_rf_hmm_single <- function(input.seq,
34 34
                            input.ph.single,
35 35
                            rf.temp = NULL,
36 36
                            tol,
@@ -43,49 +43,49 @@
Loading
43 43
  if (!inherits(input.seq, input_classes[1])) {
44 44
    stop(deparse(substitute(input.seq)), " is not an object of class 'mappoly.sequence'")
45 45
  }
46 -
  if(length(input.seq$seq.num) == 1)
46 +
  if(length(input.seq$seq.num)  ==  1)
47 47
    stop("Input sequence contains only one marker.", call. = FALSE)
48 48
  if (verbose && !capabilities("long.double") && high.prec){
49 49
    cat("You've requested high precision calculations ('high.prec = TRUE'), but your system's architecture doesn't support long double allocation ('capabilities('long.double') = FALSE'). Running in low precision mode.\n")
50 50
  }
51 51
  if(is.null(rf.temp))
52 -
    rf.temp<-rep(0.001, length(input.seq$seq.num)-1)
52 +
    rf.temp <- rep(0.001, length(input.seq$seq.num)-1)
53 53
  if(!ret.map.no.rf.estimation)
54 54
  {
55 -
    D<-get(input.seq$data.name, pos=1)$geno.dose[input.seq$seq.num,]
56 -
    dp <- get(input.seq$data.name)$dosage.p[input.seq$seq.num]
57 -
    dq <- get(input.seq$data.name)$dosage.q[input.seq$seq.num]
55 +
    D <- get(input.seq$data.name, pos = 1)$geno.dose[input.seq$seq.num,]
56 +
    dp <- get(input.seq$data.name)$dosage.p1[input.seq$seq.num]
57 +
    dq <- get(input.seq$data.name)$dosage.p2[input.seq$seq.num]
58 58
    for (j in 1:nrow(D))
59 -
      D[j, D[j, ] == input.seq$m + 1] <- dp[j] + dq[j] + 1 + as.numeric(dp[j]==0 || dq[j]==0)
59 +
      D[j, D[j, ]  ==  input.seq$ploidy + 1] <- dp[j] + dq[j] + 1 + as.numeric(dp[j] == 0 || dq[j] == 0)
60 60
61 61
    if(high.prec)
62 62
    {
63 -
      res.temp<-.Call("est_map_hmm_highprec",
64 -
                      input.seq$m,
63 +
      res.temp <- .Call("est_map_hmm_highprec",
64 +
                      input.seq$ploidy,
65 65
                      t(D),
66 66
                      lapply(input.ph.single$P, function(x) x-1),
67 67
                      lapply(input.ph.single$Q, function(x) x-1),
68 68
                      rf.temp,
69 -
                      verbose=verbose,
69 +
                      verbose = verbose,
70 70
                      max.rf.to.break.EM,
71 71
                      tol,
72 72
                      PACKAGE = "mappoly")
73 73
    } else{
74 -
      res.temp<-.Call("est_map_hmm",
75 -
                      input.seq$m,
74 +
      res.temp <- .Call("est_map_hmm",
75 +
                      input.seq$ploidy,
76 76
                      t(D),
77 77
                      lapply(input.ph.single$P, function(x) x-1),
78 78
                      lapply(input.ph.single$Q, function(x) x-1),
79 79
                      rf.temp,
80 -
                      verbose=verbose,
80 +
                      verbose = verbose,
81 81
                      max.rf.to.break.EM,
82 82
                      tol,
83 83
                      PACKAGE = "mappoly")
84 84
    }
85 -
  } else res.temp<-list(1, rf.temp)
86 -
  map<-list(seq.num=input.seq$seq.num,
87 -
            seq.rf=res.temp[[2]],
88 -
            seq.ph=input.ph.single,
89 -
            loglike=res.temp[[1]])
85 +
  } else res.temp <- list(1, rf.temp)
86 +
  map <- list(seq.num = input.seq$seq.num,
87 +
            seq.rf = res.temp[[2]],
88 +
            seq.ph = input.ph.single,
89 +
            loglike = res.temp[[1]])
90 90
  map
91 91
}

@@ -23,7 +23,7 @@
Loading
23 23
#' 
24 24
#' @examples
25 25
#'  ## tetraploid example
26 -
#'  probs.t<-calc_genoprob_dist(input.map = solcap.prior.map[[1]],
26 +
#'  probs.t <- calc_genoprob_dist(input.map = solcap.prior.map[[1]],
27 27
#'                            dat.prob = tetra.solcap.geno.dist,
28 28
#'                            verbose = TRUE)
29 29
#'  probs.t
@@ -41,7 +41,7 @@
Loading
41 41
#'     \doi{10.1534/g3.119.400378} 
42 42
#'
43 43
#' @export calc_genoprob_dist
44 -
calc_genoprob_dist<-function(input.map, dat.prob = NULL, phase.config = "best", verbose = TRUE)
44 +
calc_genoprob_dist <- function(input.map, dat.prob = NULL, phase.config = "best", verbose = TRUE)
45 45
{
46 46
  if (!inherits(input.map, "mappoly.map")) {
47 47
    stop(deparse(substitute(input.map)), " is not an object of class 'mappoly.map'")
@@ -51,52 +51,52 @@
Loading
51 51
  }
52 52
  ## choosing the linkage phase configuration
53 53
  LOD.conf <- get_LOD(input.map, sorted = FALSE)
54 -
  if(phase.config == "best") {
54 +
  if(phase.config  ==  "best") {
55 55
    i.lpc <- which.min(LOD.conf)
56 56
  } else if (phase.config > length(LOD.conf)) {
57 57
    stop("invalid linkage phase configuration")
58 58
  } else i.lpc <- phase.config
59 59
  if(is.null(dat.prob)){
60 -
    if(nrow(get(input.map$info$data.name, pos=1)$geno)==get(input.map$info$data.name, pos=1)$n.mrk){
60 +
    if(nrow(get(input.map$info$data.name, pos = 1)$geno) == get(input.map$info$data.name, pos = 1)$n.mrk){
61 61
      stop("
62 62
            The dataset associated to 'input.map'
63 63
            contains no genotypic probability distribution.
64 64
            Please provide provide dataset in argument 
65 65
           'dat.prob'.
66 66
           ")
67 67
    } else {
68 -
      dat.prob <- get(input.map$info$data.name, pos=1)
68 +
      dat.prob <- get(input.map$info$data.name, pos = 1)
69 69
    }
70 70
  }
71 -
  mrk<-NULL
71 +
  mrk <- NULL
72 72
  original.map.mrk <- input.map$info$mrk.names
73 73
  dat.prob.pos <- match(original.map.mrk, dat.prob$mrk.names)
74 -
  which.is.na<-which(is.na(dat.prob.pos))
74 +
  which.is.na <- which(is.na(dat.prob.pos))
75 75
  if(length(which.is.na) > 0)
76 76
    stop("Markers", original.map.mrk[which.is.na], "are not present in the 'dat.prob' object")
77 -
  temp.map<-input.map
78 -
  temp.map$info$seq.num <- temp.map$maps[[i.lpc]]$seq.num<-dat.prob.pos
79 -
  names(temp.map$maps[[i.lpc]]$seq.ph$P)<-names(temp.map$maps[[i.lpc]]$seq.ph$Q)<-dat.prob.pos
77 +
  temp.map <- input.map
78 +
  temp.map$info$seq.num <- temp.map$maps[[i.lpc]]$seq.num <- dat.prob.pos
79 +
  names(temp.map$maps[[i.lpc]]$seq.ph$P) <- names(temp.map$maps[[i.lpc]]$seq.ph$Q) <- dat.prob.pos
80 80
  if(!all(sort(get(temp.map$info$data.name, pos = 1)$ind.names) %in% sort(get(input.map$info$data.name, pos = 1)$ind.names)))
81 81
    stop("The individuals are different in the new and original datasets")
82 -
  geno<-subset(dat.prob$geno, mrk%in%original.map.mrk)
83 -
  geno.new<-NULL
82 +
  geno <- subset(dat.prob$geno, mrk%in%original.map.mrk)
83 +
  geno.new <- NULL
84 84
  for(i in unique(geno$ind))
85 -
    geno.new<-rbind(geno.new, geno[geno[,"ind"] == i, ][match(original.map.mrk, geno[,"mrk"]),])
85 +
    geno.new <- rbind(geno.new, geno[geno[,"ind"]  ==  i, ][match(original.map.mrk, geno[,"mrk"]),])
86 86
  g <- as.double(t(geno.new[, -c(1:2)]))
87 -
  m = as.numeric(temp.map$info$m)
87 +
  ploidy = as.numeric(temp.map$info$ploidy)
88 88
  n.mrk = as.numeric(temp.map$info$n.mrk)
89 89
  n.ind = get(temp.map$info$data.name, pos = 1)$n.ind
90 90
  p = as.numeric(unlist(temp.map$maps[[1]]$seq.ph$P))
91 91
  dp = as.numeric(cumsum(c(0, sapply(temp.map$maps[[1]]$seq.ph$P, function(x) sum(length(x))))))
92 92
  q = as.numeric(unlist(temp.map$maps[[1]]$seq.ph$Q))
93 93
  dq = as.numeric(cumsum(c(0, sapply(temp.map$maps[[1]]$seq.ph$Q, function(x) sum(length(x))))))
94 94
  rf = temp.map$maps[[1]]$seq.rf
95 -
  indnames<-dat.prob$ind.names
96 -
  res.temp <-
95 +
  ind.names <- dat.prob$ind.names
96 +
  res.temp  <- 
97 97
    .Call(
98 98
      "calc_genoprob_prior",
99 -
      as.numeric(m),
99 +
      as.numeric(ploidy),
100 100
      as.numeric(n.mrk),
101 101
      as.numeric(n.ind),
102 102
      as.numeric(p),
@@ -105,15 +105,15 @@
Loading
105 105
      as.numeric(dq),
106 106
      as.double(g),
107 107
      as.double(rf),
108 -
      as.numeric(rep(0, choose(m, m/2)^2 * n.mrk * n.ind)),
108 +
      as.numeric(rep(0, choose(ploidy, ploidy/2)^2 * n.mrk * n.ind)),
109 109
      as.double(0),
110 110
      as.numeric(verbose),
111 111
      PACKAGE = "mappoly"
112 112
    )
113 113
  if (verbose) cat("\n")
114 -
  dim(res.temp[[1]])<-c(choose(m,m/2)^2,n.mrk,n.ind)
115 -
  dimnames(res.temp[[1]])<-list(kronecker(apply(combn(letters[1:m],m/2),2, paste, collapse=""),
116 -
                                          apply(combn(letters[(m+1):(2*m)],m/2),2, paste, collapse=""), paste, sep=":"),
117 -
                                original.map.mrk, indnames)
118 -
  structure(list(probs = res.temp[[1]], map = create_map(input.map)), class="mappoly.genoprob")
114 +
  dim(res.temp[[1]]) <- c(choose(ploidy,ploidy/2)^2,n.mrk,n.ind)
115 +
  dimnames(res.temp[[1]]) <- list(kronecker(apply(combn(letters[1:ploidy],ploidy/2),2, paste, collapse = ""),
116 +
                                          apply(combn(letters[(ploidy+1):(2*ploidy)],ploidy/2),2, paste, collapse = ""), paste, sep = ":"),
117 +
                                original.map.mrk, ind.names)
118 +
  structure(list(probs = res.temp[[1]], map = create_map(input.map)), class = "mappoly.genoprob")
119 119
}

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Click to load this diff.
Loading diff...

Everything is accounted for!

No changes detected that need to be reviewed.
What changes does Codecov check for?
Lines, not adjusted in diff, that have changed coverage data.
Files that introduced coverage data that had none before.
Files that have missing coverage data that once were tracked.
Files Coverage
R -0.60% 69.38%
src 76.71%
Project Totals (59 files) 71.03%
Loading