# \references{
# [1] P. Neuvial, \emph{Improved CN estimation and detection for SNP arrays}, private notes, 2009-12-01.\cr
# }
##
## NOMENCLATURE?!?
##
## ascnT [Jx2], ascnN [Jx2] => (thetaA, thetaB)
##
## pscnT [Jx2], pscnN [Jx2] => (theta, beta)
##
## pscnT [Jx2], pscnN [Jx2] => (theta1, theta2)
##
## data <- calculatePairedPSCNByGenotype(pscnT, pscnN)
## fit <- segmentByPairedPSCBS(data, tbn=FALSE)
## fit <- callROH(fit)
##
setMethodS3("calculatePairedPSCNByGenotype", "numeric", function(thetaT, betaT, thetaN, betaN, muN=callNaiveGenotypes(betaN, ...), flavor=c("v1"), ...) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'thetaT':
thetaT <- as.numeric(thetaT)
# Argument 'betaT':
betaT <- as.numeric(betaT)
J <- length(thetaT)
if (length(betaT) != J) {
stop("The length of arguments 'betaT' and 'thetaT' differ: ",
length(betaT), " != ", J)
}
# Argument 'thetaN':
thetaN <- as.numeric(thetaN)
if (length(thetaN) != J) {
stop("The length of arguments 'thetaN' and 'thetaT' differ: ",
length(thetaN), " != ", J)
}
# Argument 'betaN':
betaN <- as.numeric(betaN)
if (length(betaN) != J) {
stop("The length of arguments 'betaN' and 'thetaT' differ: ",
length(betaN), " != ", J)
}
# Argument 'muN':
if (is.null(muN)) {
muN <- callNaiveGenotypes(betaN, ...)
}
if (length(muN) != J) {
stop("The length of arguments 'muN' and 'thetaT' differ: ",
length(muN), " != ", J)
}
knownGenotypes <- c(0, 1/2, 1, NA)
unknown <- which(!is.element(muN, knownGenotypes))
n <- length(unknown)
if (n > 0) {
unknown <- unique(muN[unknown])
unknownStr <- paste(unknown, collapse=", ")
stop("Argument 'muN' contains unknown values: ", unknownStr)
}
# Argument 'flavor':
flavor <- match.arg(flavor)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Calculate CT for all loci [and record (betaT,muN)]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
CTx <- 2 * thetaT/thetaN
betaTx <- betaT
muNx <- muN
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Calculate (CT,betaT) for SNPs
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Identify set to be updated (our best what are SNPs)
isSNP <- which(is.finite(betaN) & is.finite(muN))
# Subset
thetaT <- thetaT[isSNP]
betaT <- betaT[isSNP]
thetaN <- thetaN[isSNP]
betaN <- betaN[isSNP]
muN <- muN[isSNP]
# Transform to (thetaA, thetaB) for tumors and normals (only for SNPs)
thetaTB <- thetaT * betaT
thetaTA <- thetaT - thetaTB
thetaNB <- thetaN * betaN
thetaNA <- thetaN - thetaNB
# Identify genotypes
isAA <- (muN == 0)
isAB <- (muN == 1/2)
isBB <- (muN == 1)
# The true allele-specific CNs in the normal (CNA,CNB)
CNA <- 2*isAA + 1*isAB
CNB <- 2*isBB + 1*isAB
# Calculate CT stratified by genotype
CTA <- thetaTA/thetaNA*CNA
CTB <- thetaTB/thetaNB*CNB
CT <- CTA + CTB
betaT <- CTB / CT
# Sanity checks
.stop_if_not(all(betaT[isAA] == 0))
.stop_if_not(all(betaT[isBB] == 1))
.stop_if_not(all(CTA[isBB] == 0))
.stop_if_not(all(CTB[isAA] == 0))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Expand to all loci
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Update SNPs
CTx[isSNP] <- CT
betaTx[isSNP] <- betaT
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Return calculate data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Sanity check
.stop_if_not(length(CTx) == J)
.stop_if_not(length(betaTx) == J)
.stop_if_not(length(muNx) == J)
list(CT=CTx, betaT=betaTx, muN=muNx)
}) # calculatePairedPSCNByGenotype()
setMethodS3("calculatePairedPSCNByGenotype", "data.frame", function(data, ...) {
res <- calculatePairedPSCNByGenotype(thetaT=data$thetaT, betaT=data$betaT, thetaN=data$thetaN, betaN=data$betaN, muN=data$muN, ...)
data$CT <- res$CT
data$betaT <- res$betaT
data$muN <- res$muN
data
})
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# TO BE RETIRED
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
# \references{
# [1] P. Neuvial, \emph{Improved CN estimation and detection for SNP arrays}, private notes, 2009-12-01.\cr
# }
setMethodS3("calculateTumorPSCNByGenotypeUsingCTCN", "numeric", function(CT, betaT, betaN, muN=callNaiveGenotypes(betaN, ...), flavor=c("v1"), ...) {
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Validate arguments
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Argument 'CT':
CT <- as.numeric(CT)
# Argument 'betaT':
betaT <- as.numeric(betaT)
J <- length(CT)
if (length(betaT) != J) {
stop("The length of arguments 'betaT' and 'CT' differ: ",
length(betaT), " != ", J)
}
# Argument 'betaN':
betaN <- as.numeric(betaN)
if (length(betaN) != J) {
stop("The length of arguments 'betaN' and 'CT' differ: ",
length(betaN), " != ", J)
}
# Argument 'muN':
if (is.null(muN)) {
muN <- callNaiveGenotypes(betaN, ...)
}
if (length(muN) != J) {
stop("The length of arguments 'muN' and 'CT' differ: ",
length(muN), " != ", J)
}
knownGenotypes <- c(0, 1/2, 1, NA)
unknown <- which(!is.element(muN, knownGenotypes))
n <- length(unknown)
if (n > 0) {
unknown <- unique(muN[unknown])
unknownStr <- paste(unknown, collapse=", ")
stop("Argument 'muN' contains unknown values: ", unknownStr)
}
# Argument 'flavor':
flavor <- match.arg(flavor)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Calculate CT for all loci [and record (betaT,muN)]
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
CTx <- CT
betaTx <- betaT
muNx <- muN
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Calculate (CT,betaT) for SNPs
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Identify set to be updated (our best what are SNPs)
isSNP <- which(is.finite(betaN) & is.finite(muN))
# Subset
CT <- CT[isSNP]
betaT <- betaT[isSNP]
betaN <- betaN[isSNP]
muN <- muN[isSNP]
# Identify genotypes
isAA <- (muN == 0)
isAB <- (muN == 1/2)
isBB <- (muN == 1)
# True (CNA,CNB)
CNA <- 2*isAA + 1*isAB
CNB <- 2*isBB + 1*isAB
# Calculate (CTA,CTB) stratified by genotype
alphaT <- 1 - betaT
alphaN <- 1 - betaN
kappaA <- (alphaT/alphaN) * CT/2
kappaB <- ( betaT/ betaN) * CT/2
CTA <- kappaA * CNA
CTB <- kappaB * CNB
CT <- CTA + CTB
betaT <- CTB / CT
# Sanity checks
.stop_if_not(all(betaT[isAA] == 0, na.rm=TRUE))
.stop_if_not(all(betaT[isBB] == 1, na.rm=TRUE))
.stop_if_not(all(CTA[isBB] == 0, na.rm=TRUE))
.stop_if_not(all(CTB[isAA] == 0, na.rm=TRUE))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Expand to all loci
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Update SNPs
CTx[isSNP] <- CT
betaTx[isSNP] <- betaT
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Return calculate data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Sanity check
.stop_if_not(length(CTx) == J)
.stop_if_not(length(betaTx) == J)
.stop_if_not(length(muNx) == J)
list(CT=CTx, betaT=betaTx, muN=muNx)
}) # calculateTumorPSCNByGenotypeUsingCTCN()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.