Nothing
################################################################################
# TODO LIST
# TODO: Re-write the whole function when a preferred model has been found.
# The function has become very complex and messy...
# TODO: Option to perform multiple random selections in one go. Which would mean
# multiple 'ModelX' columns.
################################################################################
# CHANGE LOG (last 20 changes)
# 09.07.2022: Fixed "...URLs which should use \doi (with the DOI name only)".
# 24.08.2018: Removed unused variables.
# 06.08.2017: Added audit trail.
# 18.09.2016: Fixed dataset saved to attributes.
# 15.09.2016: Implemented new filterProfile function to remove sex markers and qs.
# 29.06.2016: Added option to remove sex markers and quality sensor.
# 09.01.2016: Added more attributes to result.
# 07.12.2015: Fixed reference sample name subsetting bug.
# 05.12.2015: More information in 'stop' messages.
# 'warning' for unhandled combinations changed to 'stop'.
# 05.10.2015: Added attributes to result.
# 28.09.2015: Remove rows with missing alleles from the reference dataset.
# 11.09.2015: Handle reference allele is NA.
# 28.08.2015: Added importFrom
# 26.06.2015: More precise warning messages (include sample name and marker).
# 15.12.2014: Changed parameter names to format: lower.case
# 20.01.2014: Changed 'saveImage_gui' for 'ggsave_gui'.
# 16.01.2014: Adde option for selection of one or more scoring methods.
# 16.01.2014: Changed names 'Model[]'/'Method[]'.
# 02.12.2013: Changed name 'Hybrid'/'Hybrid.Ph' to 'ModelL'/'ModelL.Ph'.
# 13.11.2013: Concurrently score both Random, Allele1, and Allele2.
# 07.11.2013: Fixed dropout check for homozygous loci in 'Hybrid' method.
#' @title Calculate Drop-out Events
#'
#' @description
#' Calculate drop-out events (allele and locus) and records the surviving peak height.
#'
#' @details
#' Calculates drop-out events. In case of allele dropout the peak height of the
#' surviving allele is given. Homozygous alleles in the reference set can be
#' either single or double notation (X or X X). Markers present in the
#' reference set but not in the data set will be added to the result.
#' NB! 'Sample.Name' in 'ref' must be unique core name of replicate sample
#' names in 'data'.
#' Use \code{checkSubset} to make sure subsetting works as intended.
#' There are options to remove sex markers and quality sensors from analysis.
#'
#' NB! There are several methods of scoring drop-out events for regression.
#' Currently the 'MethodX', 'Method1', and 'Method2' are endorsed by the DNA
#' commission (see Appendix B in ref 1). However, an alternative method is to
#' consider the whole locus and score drop-out if any allele is missing.
#'
#' Explanation of the methods:
#' Dropout - all alleles are scored according to AT. This is pure observations
#' and is not used for modeling.
#' MethodX - a random reference allele is selected and drop-out is scored in
#' relation to the the partner allele.
#' Method1 - the low molecular weight allele is selected and drop-out is
#' scored in relation to the partner allele.
#' Method2 - the high molecular weight allele is selected and drop-out is
#' scored in relation to the partner allele.
#' MethodL - drop-out is scored per locus i.e. drop-out if any allele has
#' dropped out.
#'
#' Method X/1/2 records the peak height of the partner allele to be used as
#' the explanatory variable in the logistic regression. The locus method L also
#' do this when there has been a drop-out, if not the the mean peak height for
#' the locus is used. Peak heights for the locus method are stored in a
#' separate column.
#'
#' @param data data frame in GeneMapper format containing at least a column 'Allele'.
#' @param ref data frame in GeneMapper format.
#' @param threshold numeric, threshold in RFU defining a dropout event.
#' Default is 'NULL' and dropout is scored purely on the absence of a peak.
#' @param method character vector, specifying which scoring method(s) to use.
#' Method 'X' for random allele, '1' or '2' for the low/high molecular weight allele,
#' and 'L' for the locus method (the option is case insensitive).
#' @param ignore.case logical, default TRUE for case insensitive.
#' @param sex.rm logical, default FALSE to include sex markers in the analysis.
#' @param qs.rm logical, default TRUE to exclude quality sensors from the analysis.
#' @param kit character, required if sex.rm=TRUE or qs.rm=TRUE to define the kit.
#' @param debug logical indicating printing debug information.
#'
#' @return data.frame with columns 'Sample.Name', 'Marker', 'Allele', 'Height', 'Dropout',
#' 'Rfu', 'Heterozygous', and 'Model'.
#' Dropout: 0 indicate no dropout, 1 indicate allele dropout, and 2 indicate locus dropout.
#' Rfu: height of surviving allele.
#' Heterozygous: 1 for heterozygous and 0 for homozygous.
#' And any of the following containing the response (or explanatory) variable used
#' for modeling by logistic regression in function \code{modelDropout}:
#' 'MethodX', 'Method1', 'Method2', 'MethodL' and 'MethodL.Ph'.
#'
#' @export
#'
#' @importFrom utils str
#'
#' @references
#' Peter Gill et.al.,
#' DNA commission of the International Society of Forensic Genetics:
#' Recommendations on the evaluation of STR typing results that may
#' include drop-out and/or drop-in using probabilistic methods,
#' Forensic Science International: Genetics, Volume 6, Issue 6, December 2012,
#' Pages 679-688, ISSN 1872-4973, 10.1016/j.fsigen.2012.06.002.
#' \doi{10.1016/j.fsigen.2012.06.002}
#' @references
#' Peter Gill, Roberto Puch-Solis, James Curran,
#' The low-template-DNA (stochastic) threshold-Its determination relative to
#' risk analysis for national DNA databases,
#' Forensic Science International: Genetics, Volume 3, Issue 2, March 2009,
#' Pages 104-111, ISSN 1872-4973, 10.1016/j.fsigen.2008.11.009.
#' \doi{10.1016/j.fsigen.2008.11.009}
#'
#' @examples
#' data(set4)
#' data(ref4)
#' drop <- calculateDropout(data = set4, ref = ref4, kit = "ESX17", ignore.case = TRUE)
calculateDropout <- function(data, ref, threshold = NULL, method = c("1", "2", "X", "L"),
ignore.case = TRUE, sex.rm = FALSE, qs.rm = TRUE,
kit = NULL, debug = FALSE) {
if (debug) {
print(paste("IN:", match.call()[[1]]))
}
# CHECK DATA ----------------------------------------------------------------
# Check dataset columns.
if (!any(grepl("Sample.Name", names(data)))) {
stop("'data' must contain a column 'Sample.Name'.",
call. = TRUE
)
}
if (!any(grepl("Marker", names(data)))) {
stop("'data' must contain a column 'Marker'.",
call. = TRUE
)
}
if (!any(grepl("Allele", names(data)))) {
stop("'data' must contain a column 'Allele'.",
call. = TRUE
)
}
if (!any(grepl("Height", names(data)))) {
stop("'data' must contain a column 'Height'.",
call. = TRUE
)
}
# Check dataset for factors.
if (is.factor(data$Sample.Name)) {
message("Data 'Sample.Name' is factors. Converting to character!")
data$Sample.Name <- as.character(data$Sample.Name)
}
if (is.factor(data$Marker)) {
message("Data 'Marker' is factors. Converting to character!")
data$Marker <- as.character(data$Marker)
}
# Check reference dataset.
if (!any(grepl("Sample.Name", names(ref)))) {
stop("'ref' must contain a column 'Sample.Name'.",
call. = TRUE
)
}
if (!any(grepl("Marker", names(ref)))) {
stop("'ref' must contain a column 'Marker'.",
call. = TRUE
)
}
if (!any(grepl("Allele", names(ref)))) {
stop("'ref' must contain a column 'Allele'.",
call. = TRUE
)
}
# Check dataset for factors.
if (is.factor(ref$Sample.Name)) {
message("Reference 'Sample.Name' is factors. Converting to character!")
ref$Sample.Name <- as.character(ref$Sample.Name)
}
if (is.factor(ref$Marker)) {
message("Reference 'Marker' is factors. Converting to character!")
ref$Marker <- as.character(ref$Marker)
}
# Check if slim format.
if (sum(grepl("Allele", names(ref)) > 1)) {
stop("'ref' must be in 'slim' format.",
call. = TRUE
)
}
if (sum(grepl("Allele", names(data)) > 1)) {
stop("'data' must be in 'slim' format.",
call. = TRUE
)
}
# Check if character data.
if (!is.character(ref$Allele)) {
message("'Allele' must be character. 'ref' converted!")
data$Allele <- as.character(data$Allele)
}
if (!is.character(data$Allele)) {
message("'Allele' must be character. 'data' converted!")
data$Allele <- as.character(data$Allele)
}
# Check if numeric data.
if (!is.numeric(data$Height)) {
message("'Height' must be numeric. 'data' converted!")
data$Height <- as.numeric(data$Height)
}
# Check threshold parameter.
if (!is.null(threshold)) {
if (!is.numeric(threshold)) {
message("'threshold' must be numeric. 'threshold' converted!")
threshold <- as.numeric(threshold)
}
if (is.na(threshold)) {
stop("'threshold' must be numeric.",
call. = TRUE
)
}
}
# Check logical arguments.
if (!is.logical(ignore.case)) {
stop("'ignore.case' must be logical.")
}
if (!is.logical(sex.rm)) {
stop("'sex.rm' must be logical.")
}
if (!is.logical(qs.rm)) {
stop("'qs.rm' must be logical.")
}
# Check kit.
if (!is.null(kit)) {
if (is.null(nrow(getKit(kit = kit)))) {
stop("'kit' was not found in the kit definition file.")
}
}
# PREPARE -------------------------------------------------------------------
# Remove sex markers and quality sensors from dataset.
if (sex.rm || qs.rm) {
message("Removing gender markers and/or quality sensors from dataset...")
data <- filterProfile(
data = data, filter.allele = FALSE,
sex.rm = sex.rm, qs.rm = qs.rm, kit = kit,
debug = debug
)
message("Removing gender markers and/or quality sensors from reference dataset...")
ref <- filterProfile(
data = ref, filter.allele = FALSE,
sex.rm = sex.rm, qs.rm = qs.rm, kit = kit,
debug = debug
)
}
# Check for NA alleles (this must be after 'remove sex markers').
if (any(is.na(ref$Allele))) {
# Remove markers with no allele in reference dataset.
tmp1 <- nrow(ref)
ref <- ref[!is.na(ref$Allele), ]
tmp2 <- nrow(ref)
message(paste("Removed", tmp1 - tmp2, "NA rows in 'ref'. This may be DYS markers in female profiles."))
}
# ANALYZE -------------------------------------------------------------------
# Create vectors for temporary.
samplesVec <- character()
markersVec <- character()
allelesVec <- character()
heightsVec <- numeric()
dropoutVec <- numeric()
rfuVec <- numeric()
hetVec <- numeric()
methodXVec <- numeric()
method1Vec <- numeric()
method2Vec <- numeric()
methodLVec <- numeric()
methodLPhVec <- numeric()
samplesTmp <- NULL
markersTmp <- NULL
allelesTmp <- NULL
heightsTmp <- NULL
dropoutTmp <- NULL
methodLTmp <- NULL
methodLPh <- NULL
rfuTmp <- NULL
hetTmp <- NULL
pass <- TRUE
# Get threshold if not given.
if (is.null(threshold)) {
threshold <- min(data$Height, na.rm = TRUE)
}
# Get sample names.
sampleNamesRef <- unique(ref$Sample.Name)
# Loop through all reference samples.
for (r in seq(along = sampleNamesRef)) {
# Select current subsets.
if (ignore.case) {
selectedSamples <- grepl(
toupper(sampleNamesRef[r]),
toupper(data$Sample.Name)
)
selectedRefs <- grepl(
paste("\\b",
toupper(sampleNamesRef[r]),
"\\b",
sep = ""
),
toupper(ref$Sample.Name)
)
} else {
selectedSamples <- grepl(
sampleNamesRef[r],
data$Sample.Name
)
selectedRefs <- grepl(
paste("\\b",
sampleNamesRef[r],
"\\b",
sep = ""
),
ref$Sample.Name
)
}
# Extract the results from the current reference sample.
dataSubset <- data[selectedSamples, ]
refSubset <- ref[selectedRefs, ]
# Get markers for current ref sample.
# NB! Needed for handling mixed typing kits.
markers <- unique(refSubset$Marker)
# Get sample names in subset.
sampleNames <- unique(dataSubset$Sample.Name)
# Loop through all individual samples.
# NB! Needed for detection of locus dropout.
for (s in seq(along = sampleNames)) {
# Extract the result for the current sample.
selectedReplicate <- data$Sample.Name == sampleNames[s]
dataSample <- data[selectedReplicate, ]
# Loop through all markers.
for (m in seq(along = markers)) {
# Get reference alleles and calculate expected number of alleles.
refAlleles <- unique(refSubset$Allele[refSubset$Marker == markers[m]])
expected <- length(refAlleles)
# Get sample alleles and peak heights.
selectMarker <- dataSample$Marker == markers[m]
dataAlleles <- dataSample$Allele[selectMarker]
dataHeight <- dataSample$Height[selectMarker]
# Get data mathcing ref alleles and their heights.
matching <- dataAlleles %in% refAlleles
matchedAlleles <- dataAlleles[matching]
peakHeight <- dataHeight[matching]
if (debug) {
print("refAlleles")
print(refAlleles)
print("expected")
print(expected)
print("matching")
print(matching)
print("matchedAlleles")
print(matchedAlleles)
print("peakHeight")
print(peakHeight)
}
# Count the number of observed alleles.
observed <- length(matchedAlleles)
# Score dropout for modeling -----------------------------------------
# Reset variables.
methodXTmp <- NULL
method1Tmp <- NULL
method2Tmp <- NULL
methodLTmp <- NULL
methodLPh <- NULL
if (observed == 1 || observed == 2) {
# Check expected number of alleles.
if (expected == 1) {
# Expected homozygous.
methodXTmp <- NA
method1Tmp <- NA
method2Tmp <- NA
} else if (expected == 2) {
# Expected heterozygous.
# Randomly choose an allele.
random <- sample(c(1, 2), 2, replace = FALSE)
selected0 <- match(refAlleles[random[1]], matchedAlleles)
partner0 <- match(refAlleles[random[2]], matchedAlleles)
# Choose Allele 1
selected1 <- match(refAlleles[1], matchedAlleles)
partner1 <- match(refAlleles[2], matchedAlleles)
# Choose Allele 2
selected2 <- match(refAlleles[2], matchedAlleles)
partner2 <- match(refAlleles[1], matchedAlleles)
# SCORE RANDOM ALLELE ---------------------------------------BEGIN-
if ("X" %in% toupper(method)) {
# Dropout if not partner.
if (is.na(partner0)) {
# Score as dropout.
modeldrop <- 1
methodXTmp <- modeldrop
} else if (peakHeight[partner0] < threshold) {
# If both alleles, check if below AT.
# Score as dropout.
modeldrop <- 1
# Save one or two entries.
if (observed == 1) {
methodXTmp <- modeldrop
} else if (observed == 2) {
methodXTmp[selected0] <- modeldrop
methodXTmp[partner0] <- NA
}
} else {
# Partner peak is above AT.
# Score as NO dropout.
modeldrop <- 0
# Save one or two entries.
if (observed == 1) {
methodXTmp <- modeldrop
} else if (observed == 2) {
methodXTmp[selected0] <- modeldrop
methodXTmp[partner0] <- NA
}
}
}
# SCORE RANDOM ALLELE -----------------------------------------END-
# SCORE ALLELE 1 --------------------------------------------BEGIN-
if ("1" %in% toupper(method)) {
# Dropout if not partner.
if (is.na(partner1)) {
# Score as dropout.
modeldrop <- 1
method1Tmp <- modeldrop
} else if (peakHeight[partner1] < threshold) {
# If both alleles, check if below AT.
# Score as dropout.
modeldrop <- 1
# Save one or two entries.
if (observed == 1) {
method1Tmp <- modeldrop
} else if (observed == 2) {
method1Tmp[selected1] <- modeldrop
method1Tmp[partner1] <- NA
}
} else {
# Partner peak is above AT.
# Score as NO dropout.
modeldrop <- 0
# Save one or two entries.
if (observed == 1) {
method1Tmp <- modeldrop
} else if (observed == 2) {
method1Tmp[selected1] <- modeldrop
method1Tmp[partner1] <- NA
}
}
}
# SCORE ALLELE 1 ----------------------------------------------END-
# SCORE ALLELE 2 --------------------------------------------BEGIN-
if ("2" %in% toupper(method)) {
# Dropout if not partner.
if (is.na(partner2)) {
# Score as dropout.
modeldrop <- 1
method2Tmp <- modeldrop
} else if (peakHeight[partner2] < threshold) {
# If both alleles, check if below AT.
# Score as dropout.
modeldrop <- 1
# Save one or two entries.
if (observed == 1) {
method2Tmp <- modeldrop
} else if (observed == 2) {
method2Tmp[selected2] <- modeldrop
method2Tmp[partner2] <- NA
}
} else {
# Partner peak is above AT.
# Score as NO dropout.
modeldrop <- 0
# Save one or two entries.
if (observed == 1) {
method2Tmp <- modeldrop
} else if (observed == 2) {
method2Tmp[selected2] <- modeldrop
method2Tmp[partner2] <- NA
}
}
}
}
# SCORE ALLELE 2 ----------------------------------------------END-
# SCORE LOCUS -----------------------------------------------BEGIN-
if ("L" %in% toupper(method)) {
# Check expected number of alleles.
if (expected == 1) {
# Expected homozygous.
if (peakHeight < threshold) {
# Dropout.
methodLTmp <- 2
methodLPh <- NA
} else {
# No dropout.
methodLTmp <- 0
methodLPh <- NA
}
} else if (expected == 2) {
# Expected heterozygous.
# Marker approach
selA1 <- match(refAlleles[1], matchedAlleles)
selA2 <- match(refAlleles[2], matchedAlleles)
if (debug) {
print("Peak heights:")
print(peakHeight[selA1])
print(peakHeight[selA2])
}
if (is.na(peakHeight[selA1])) {
drop1 <- TRUE
} else {
drop1 <- peakHeight[selA1] < threshold
}
if (is.na(peakHeight[selA2])) {
drop2 <- TRUE
} else {
drop2 <- peakHeight[selA2] < threshold
}
passingAlleles <- sum(!drop1, !drop2, na.rm = TRUE)
if (debug) {
print("drop1")
print(drop1)
print("drop2")
print(drop2)
print("passingAlleles:")
print(passingAlleles)
print("observed:")
print(observed)
}
# Check for any dropouts.
if (any(drop1, drop2, na.rm = TRUE)) {
if (debug) {
print("At least one dropout!")
}
# Save one or two entries.
if (passingAlleles == 1) {
if (observed == 1) {
methodLTmp <- 1
methodLPh <- max(peakHeight[selA1], peakHeight[selA2], na.rm = TRUE)
} else if (observed == 2) {
methodLTmp[1] <- 1
methodLTmp[2] <- NA
methodLPh[1] <- max(peakHeight[selA1], peakHeight[selA2], rm.na = TRUE)
methodLPh[2] <- NA
} else {
stop(
paste("Sample: ", sampleNames[s],
", Marker: ", markers[m],
" - unhandled number of observed alleles",
" (observed = ", observed,
", matchedAlleles = ", paste(matchedAlleles, collapse = "/"),
", passingAlleles = ", paste(passingAlleles, collapse = "/"),
")",
sep = ""
),
call. = TRUE
)
}
} else if (passingAlleles == 2) {
methodLTmp[1] <- 0
methodLTmp[2] <- NA
methodLPh[1] <- mean(c(peakHeight[selA1], peakHeight[selA2]))
methodLPh[2] <- NA
} else if (passingAlleles == 0) {
if (observed == 1) {
methodLTmp <- 2
methodLPh <- NA
} else if (observed == 2) {
methodLTmp[1] <- 2
methodLTmp[2] <- NA
methodLPh[1] <- NA
methodLPh[2] <- NA
} else {
stop(
paste("Sample: ", sampleNames[s],
" Marker: ", markers[m],
" - unhandled number of observed alleles",
" (observed = ", observed,
", matchedAlleles = ", paste(matchedAlleles, collapse = "/"),
", passingAlleles = ", paste(passingAlleles, collapse = "/"),
")",
sep = ""
),
call. = TRUE
)
}
} else {
stop(
paste("Sample: ", sampleNames[s],
" Marker: ", markers[m],
" - unhandled number of observed alleles > AT",
" (passingAlleles = ", paste(passingAlleles, collapse = "/"),
", matchedAlleles = ", paste(matchedAlleles, collapse = "/"),
")",
sep = ""
),
call. = TRUE
)
}
} else { # No dropout or NA.
if (debug) {
print("No dropout:")
print("Observed:")
print(observed)
}
# Save one or two entries.
if (observed == 1) {
methodLTmp <- 1
methodLPh <- max(peakHeight[selA1], peakHeight[selA2], rm.na = TRUE)
} else if (observed == 2) {
methodLTmp[1] <- 0
methodLTmp[2] <- NA
methodLPh[1] <- mean(c(peakHeight[selA1], peakHeight[selA2]))
methodLPh[2] <- NA
}
if (debug) {
print("Pk:")
print(methodLPh)
}
}
}
}
# SCORE LOCUS -------------------------------------------------END-
} else { # Zero or more peaks.
if (observed == 0) {
methodXTmp <- NA
method1Tmp <- NA
method2Tmp <- NA
methodLTmp <- 2
methodLPh <- NA
} else {
methodXTmp <- rep(NA, observed)
method1Tmp <- rep(NA, observed)
method2Tmp <- rep(NA, observed)
methodLTmp <- rep(NA, observed)
methodLPh <- rep(NA, observed)
}
}
# Score all dropouts --------------------------------------------------
# TODO: This is not needed if the 'locus' method is used. With that
# approach both regression and heat-maps can use the same data.
# Compare to threshold
if (!is.null(threshold)) {
# Mark peaks above threhold.
pass <- peakHeight >= threshold
if (length(pass) == 0) {
# This happens when there are no mathing alleles.
pass <- NULL
}
}
# Count the number of observed peaks passing the threshold.
observedPass <- sum(pass)
# Count the number of alleles that have dropped out.
dropCount <- expected - observedPass
if (debug) {
if (length(dataHeight) != length(dataAlleles)) {
print("WARNING! Different length for dataHeight and dataAlleles")
}
print("dataAlleles")
print(dataAlleles)
print("dataHeight")
print(dataHeight)
print("pass")
print(pass)
print("expected")
print(expected)
print("observedPass")
print(observedPass)
print("Sample")
print(sampleNames[s])
print("Marker")
print(markers[m])
}
# Handle locus dropout.
if (length(matchedAlleles) == 0) {
allelesTmp <- NA
heightsTmp <- NA
records <- 1
} else {
# Store all peaks.
allelesTmp <- matchedAlleles
heightsTmp <- peakHeight
records <- length(matchedAlleles)
}
# Uppdate vector.
allelesVec <- c(allelesVec, allelesTmp)
heightsVec <- c(heightsVec, heightsTmp)
methodXVec <- c(methodXVec, methodXTmp)
method1Vec <- c(method1Vec, method1Tmp)
method2Vec <- c(method2Vec, method2Tmp)
methodLVec <- c(methodLVec, methodLTmp)
methodLPhVec <- c(methodLPhVec, methodLPh)
# Samples and markers.
samplesTmp <- rep(sampleNames[s], records)
markersTmp <- rep(markers[m], records)
# Update vectors.
samplesVec <- c(samplesVec, samplesTmp)
markersVec <- c(markersVec, markersTmp)
# Indicate zygosity (1-Heterozygote, 0-Homozygote).
if (expected == 1) {
hetTmp <- rep(0, records)
het <- FALSE
} else if (expected == 2) {
hetTmp <- rep(1, records)
het <- TRUE
} else {
stop(
paste("Sample: ", sampleNames[s],
" Marker: ", markers[m],
" - unhandled number of expected alleles",
" (expected = ", expected,
", refAlleles = ", paste(refAlleles, collapse = "/"),
")",
sep = ""
),
call. = TRUE
)
}
# Update vector.
hetVec <- c(hetVec, hetTmp)
# Indicate dropout:
# 0 - for no dropout, 1 - for allele dropout, 2 - for locus dropout
if (dropCount == 0) {
dropoutTmp <- rep(0, records)
} else if (dropCount == 1 & het) {
dropoutTmp <- rep(1, records)
} else if (dropCount == 1 & !het) {
dropoutTmp <- rep(2, records)
} else if (dropCount == 2 & het) {
dropoutTmp <- rep(2, records)
} else {
stop(
paste("Sample: ", sampleNames[s],
" Marker: ", markers[m],
" - unhandled combination",
" (dropCount = ", dropCount,
", het = ", het,
")",
sep = ""
),
call. = TRUE
)
}
# Replace dropout for heights below threshold with NA or 2 if locus dropout.
if (!is.null(pass)) {
if (all(!pass)) {
dropoutTmp[!pass] <- 2
} else {
dropoutTmp[!pass] <- NA
}
}
# Update vector.
dropoutVec <- c(dropoutVec, dropoutTmp)
# Store peak height of surviving allele, or NA.
if (dropCount == 1 & observed > 0) {
rfuTmp <- peakHeight
} else {
rfuTmp <- rep(NA, records)
}
# Replace rfu for heights below threshold with NA.
if (!is.null(pass)) {
rfuTmp[!pass] <- NA
}
# Update vector.
rfuVec <- c(rfuVec, rfuTmp)
if (debug) {
print("dataAlleles")
print(dataAlleles)
print("dataHeight")
print(dataHeight)
print("pass")
print(pass)
print("expected")
print(expected)
print("observedPass")
print(observedPass)
print("Sample")
print(sampleNames[s])
print("Marker")
print(markers[m])
print("records")
print(records)
print("samplesTmp")
print(samplesTmp)
print("markersTmp")
print(markersTmp)
print("allelesTmp")
print(allelesTmp)
print("heightsTmp")
print(heightsTmp)
print("dropoutTmp")
print(dropoutTmp)
print("rfuTmp")
print(rfuTmp)
print("hetTmp")
print(hetTmp)
print("methodXTmp")
print(methodXTmp)
print("method1Tmp")
print(method1Tmp)
print("method2Tmp")
print(method2Tmp)
}
}
}
}
if (debug) {
print("samplesVec")
print(str(samplesVec))
print("markersVec")
print(str(markersVec))
print("allelesVec")
print(str(allelesVec))
print("heightsVec")
print(str(heightsVec))
print("dropooutVec")
print(str(dropoutVec))
print("rfuVec")
print(str(rfuVec))
print("hetVec")
print(str(hetVec))
print("method1Vec")
print(str(method1Vec))
print("method2Vec")
print(str(method2Vec))
print("methodLVec")
print(str(methodLVec))
print("methodLPhVec")
print(str(methodLPhVec))
}
# Create return dataframe.
dataDrop <- data.frame(
Sample.Name = samplesVec,
Marker = markersVec,
Allele = allelesVec,
Height = heightsVec,
Dropout = dropoutVec,
Rfu = rfuVec,
Heterozygous = hetVec,
stringsAsFactors = FALSE
)
# Add according to scoring method(s):
if ("X" %in% toupper(method)) {
dataDrop$MethodX <- methodXVec
}
if ("1" %in% toupper(method)) {
dataDrop$Method1 <- method1Vec
}
if ("2" %in% toupper(method)) {
dataDrop$Method2 <- method2Vec
}
if ("L" %in% toupper(method)) {
dataDrop$MethodL <- methodLVec
dataDrop$MethodL.Ph <- methodLPhVec
}
# Add attributes to result.
attr(dataDrop, which = "kit") <- kit
# Update audit trail.
dataDrop <- auditTrail(obj = dataDrop, f.call = match.call(), package = "strvalidator")
if (debug) {
print(paste("EXIT:", match.call()[[1]]))
}
return(dataDrop)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.