|
1 | 1 | ###################################################################### |
2 | 2 | # phyloqtl_scan.R |
3 | 3 | # |
4 | | -# copyright (c) 2009-2019, Karl W Broman |
5 | | -# last modified Dec, 2019 |
| 4 | +# copyright (c) 2009-2020, Karl W Broman |
| 5 | +# last modified Dec, 2020 |
6 | 6 | # first written May, 2009 |
7 | 7 | # |
8 | 8 | # This program is free software; you can redistribute it and/or |
@@ -327,16 +327,24 @@ inferredpartitions <- |
327 | 327 | lodthreshold=0 |
328 | 328 | } |
329 | 329 |
|
| 330 | + if(probthreshold >= 1 || probthreshold <= 0) { |
| 331 | + stop("probthreshold should be in (0,1)") |
| 332 | + } |
| 333 | + |
330 | 334 | output <- output[output[,1]==chr,] |
331 | 335 | output[,1] <- as.factor(as.character(output[,1])) |
332 | 336 | output <- summary(output, format="postprob") |
333 | 337 |
|
334 | 338 | if(output$maxlod < lodthreshold) return("null") |
335 | 339 |
|
336 | | - prob <- sort(output[,3:(ncol(output)-2)], decreasing=TRUE) |
| 340 | + prob <- sort(unlist(output[,3:(ncol(output)-2)]), decreasing=TRUE) |
337 | 341 | cs <- cumsum(as.numeric(prob)) |
| 342 | + if(!any(cs >= probthreshold)) { |
| 343 | + warning("No values >= probthreshold") |
| 344 | + return(NULL) |
| 345 | + } |
338 | 346 | wh <- min(which(cs >= probthreshold)) |
339 | | - names(prob)[1:wh] |
| 347 | + names(prob)[seq_len(wh)] |
340 | 348 | } |
341 | 349 |
|
342 | 350 | # end of phyloqtl_scan.R |
0 commit comments