##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
getProtocolData.Affy <- function(filenames){
scanDates <- data.frame(ScanDate=sapply(filenames, celfileDate))
rownames(scanDates) <- basename(rownames(scanDates))
protocoldata <- new("AnnotatedDataFrame",
data=scanDates,
varMetadata=data.frame(labelDescription=colnames(scanDates),
row.names=colnames(scanDates)))
return(protocoldata)
}
getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
pkgname <- getCrlmmAnnotationName(cdfName)
if(!require(pkgname, character.only=TRUE)){
suggCall <- paste("library(", pkgname, ", lib.loc='/Altern/Lib/Loc')", sep="")
msg <- paste("If", pkgname, "is installed on an alternative location, please load it manually by using", suggCall)
message(strwrap(msg))
stop("Package ", pkgname, " could not be found.")
rm(suggCall, msg)
}
loader("preprocStuff.rda", .crlmmPkgEnv, pkgname)
loader("genotypeStuff.rda", .crlmmPkgEnv, pkgname)
loader("mixtureStuff.rda", .crlmmPkgEnv, pkgname)
gns <- getVarInEnv("gns")
path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))
load(file.path(path, "snpProbes.rda"))
if(copynumber){
load(file.path(path, "cnProbes.rda"))
cnProbes <- get("cnProbes")
snpIndex <- seq(along=gns)
npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex)
featurenames <- c(gns, rownames(cnProbes))
} else featurenames <- gns
fvarlabels=c("chromosome", "position", "isSnp")
M <- matrix(NA, length(featurenames), 3, dimnames=list(featurenames, fvarlabels))
index <- match(rownames(snpProbes), rownames(M)) #only snp probes in M get assigned position
M[index, "position"] <- snpProbes[, grep("pos", colnames(snpProbes))]
M[index, "chromosome"] <- snpProbes[, grep("chr", colnames(snpProbes))]
M[index, "isSnp"] <- 1L
index <- which(is.na(M[, "isSnp"]))
M[index, "isSnp"] <- 1L
if(copynumber){
index <- match(rownames(cnProbes), rownames(M)) #only snp probes in M get assigned position
M[index, "position"] <- cnProbes[, grep("pos", colnames(cnProbes))]
M[index, "chromosome"] <- cnProbes[, grep("chr", colnames(cnProbes))]
M[index, "isSnp"] <- 0L
}
##A few of the snpProbes do not match -- I think it is chromosome Y.
M[is.na(M[, "isSnp"]), "isSnp"] <- 1L
return(new("AnnotatedDataFrame", data=data.frame(M)))
##list(snpIndex, npIndex, fns)
##crlmmOpts$snpRange <- range(snpIndex)
##crlmmOpts$npRange <- range(npIndex)
}
construct <- function(filenames, cdfName, copynumber=FALSE,
sns, verbose=TRUE, batch){
if(verbose) message("Initializing container for assay data elements alleleA, alleleB, call, callProbability, CA, CB")
if(missing(sns)) sns <- basename(filenames)
protocolData <- getProtocolData.Affy(filenames)
if(missing(batch)){
protocolData$batch <- as.numeric(as.factor(protocolData$ScanDate))
} else {
if(length(batch) != length(filenames))
stop("batch variable must be the same length as the filenames")
protocolData$batch <- batch
}
rownames(pData(protocolData)) <- sns
featureData <- getFeatureData.Affy(cdfName, copynumber=copynumber)
nr <- nrow(featureData); nc <- length(filenames)
pd <- data.frame(matrix(NA, nc, 3), row.names=sns)
colnames(pd)=c("SKW", "SNR", "gender")
pD <- new("AnnotatedDataFrame", data=pd)
alleleA <- initializeBigMatrix(name="A", nr, nc)
colnames(alleleA) <- sns
alleleB <- initializeBigMatrix(name="B", nr, nc)
colnames(alleleB) <- sns
call=initializeBigMatrix(name="call", nr, nc)
colnames(call) <- sns
callProbability=initializeBigMatrix(name="callPr", nr,nc)
colnames(callProbability) <- sns
CA=initializeBigMatrix(name="CA", nr, nc)
colnames(CA) <- sns
CB=initializeBigMatrix(name="CB", nr, nc)
colnames(CB) <- sns
callSet <- new("CNSetLM",
alleleA=alleleA,
alleleB=alleleB,
call=call,
callProbability=callProbability,
CA=CA,
CB=CB,
protocolData=protocolData,
featureData=featureData,
phenoData=pD,
annotation=cdfName)
## phenoData(callSet) <- new("AnnotatedDataFrame", data=pd)
## rownames(pData(callSet)) <- NULL
phenoData(callSet)$sampleNames <- sns
return(callSet)
}
genotype <- function(filenames,
cdfName,
batch,
mixtureSampleSize=10^5,
eps=0.1,
verbose=TRUE,
seed=1,
sns,
copynumber=FALSE,
probs=rep(1/3, 3),
DF=6,
SNRMin=5,
recallMin=10,
recallRegMin=1000,
gender=NULL,
returnParams=TRUE,
badSNP=0.7){
if(missing(cdfName)) stop("must specify cdfName")
if(!isValidCdfName(cdfName)) stop("cdfName not valid. see validCdfNames")
if(missing(sns)) sns <- basename(filenames)
## callSet contains potentially very big matrices
## More big matrices are created within snprma, that will then be removed.
callSet <- construct(filenames=filenames,
cdfName=cdfName,
copynumber=copynumber,
sns=sns,
verbose=verbose,
batch=batch)
mixtureParams <- matrix(NA, 4, length(filenames))
snp.index <- which(isSnp(callSet)==1)
batches <- splitIndicesByLength(1:ncol(callSet), ocSamples())
iter <- 1
## for(j in batches){
j <- unlist(batches)
if(verbose) message("Batch ", iter, " of ", length(batches))
snprmaRes <- snprma(filenames=filenames[j],
mixtureSampleSize=mixtureSampleSize,
fitMixture=TRUE,
eps=eps,
verbose=verbose,
seed=seed,
cdfName=cdfName,
sns=sns[j])
stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
pData(callSet)$SKW[j] <- snprmaRes$SKW
pData(callSet)$SNR[j] <- snprmaRes$SNR
suppressWarnings(A(callSet)[snp.index, j] <- snprmaRes[["A"]])
suppressWarnings(B(callSet)[snp.index, j] <- snprmaRes[["B"]])
mixtureParams[, j] <- snprmaRes$mixtureParams
rm(snprmaRes); gc()
if(copynumber){
np.index <- which(isSnp(callSet) == 0)
cnrmaRes <- cnrma(filenames=filenames[j],
cdfName=cdfName,
row.names=featureNames(callSet)[np.index],
sns=sns[j],
seed=seed,
verbose=verbose)
stopifnot(identical(featureNames(callSet)[np.index], rownames(cnrmaRes)))
A(callSet)[np.index, j] <- cnrmaRes
rm(cnrmaRes); gc()
}
## as.matrix needed when ffdf is used
tmp <- crlmmGT(A=as.matrix(A(callSet)[snp.index, j]),
B=as.matrix(B(callSet)[snp.index, j]),
SNR=callSet$SNR[j],
mixtureParams=mixtureParams[, j],
cdfName=annotation(callSet),
row.names=featureNames(callSet)[snp.index],
col.names=sampleNames(callSet)[j],
probs=probs,
DF=DF,
SNRMin=SNRMin,
recallMin=recallMin,
recallRegMin=recallRegMin,
gender=gender,
verbose=verbose,
returnParams=returnParams,
badSNP=badSNP)
snpCall(callSet)[snp.index, j] <- tmp[["calls"]]
snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]]
callSet$gender[j] <- tmp$gender
iter <- iter+1
## }
return(callSet)
}
genotype2 <- function(filenames,
cdfName,
batch,
mixtureSampleSize=10^5,
eps=0.1,
verbose=TRUE,
seed=1,
sns,
copynumber=FALSE,
probs=rep(1/3, 3),
DF=6,
SNRMin=5,
recallMin=10,
recallRegMin=1000,
gender=NULL,
returnParams=TRUE,
badSNP=0.7){
if(!isPackageLoaded("ff")) stop("Must load package 'ff'")
if(!copynumber){
callSet <- crlmm2(filenames=filenames,
cdfName=cdfName,
mixtureSampleSize=mixtureSampleSize,
eps=eps,
verbose=verbose,
sns=sns,
probs=probs,
DF=DF,
SNRMin=SNRMin,
recallMin=recallMin,
recallRegMin=recallRegMin,
gender=gender,
returnParams=returnParams,
badSNP=badSNP)
return(callSet)
}
if(missing(cdfName)) stop("must specify cdfName")
if(!isValidCdfName(cdfName)) stop("cdfName not valid. see validCdfNames")
if(missing(sns)) sns <- basename(filenames)
## callSet contains potentially very big matrices
callSet <- construct(filenames=filenames,
cdfName=cdfName,
copynumber=TRUE,
sns=sns,
verbose=verbose,
batch=batch)
##lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))
mixtureParams <- matrix(NA, 4, length(filenames))
snp.index <- which(isSnp(callSet)==1)
snprmaRes <- snprma2(filenames=filenames,
mixtureSampleSize=mixtureSampleSize,
fitMixture=TRUE,
eps=eps,
verbose=verbose,
seed=seed,
cdfName=cdfName,
sns=sns)
if(verbose) message("Finished preprocessing.")
open(snprmaRes[["A"]])
open(snprmaRes[["B"]])
open(snprmaRes[["SNR"]])
open(snprmaRes[["SKW"]])
open(snprmaRes[["mixtureParams"]])
if(verbose) message("Updating elements of callSet")
## bb = ocProbesets()*ncol(A)*8
## ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]], BATCHBYTES=bb)
## ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]], BATCHBYTES=bb)
##batches <- splitIndicesByLength(1:nrow(snprmaRes[["A"]]), ocProbesets())
AA <- A(callSet)
BB <- B(callSet)
for(j in 1:ncol(callSet)){
AA[snp.index, j] <- snprmaRes[["A"]][, j]
BB[snp.index, j] <- snprmaRes[["B"]][, j]
}
if(verbose) message("Finished updating elements of callSet")
stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
pData(callSet)$SKW <- snprmaRes$SKW
pData(callSet)$SNR <- snprmaRes$SNR
mixtureParams <- snprmaRes$mixtureParams
np.index <- which(isSnp(callSet) == 0)
cnrmaRes <- cnrma2(A=AA,
filenames=filenames,
row.names=featureNames(callSet)[np.index],
cdfName=cdfName,
sns=sns,
seed=seed,
verbose=verbose)
rm(cnrmaRes); gc()
## as.matrix needed when ffdf is used
tmp <- crlmmGT2(A=snprmaRes[["A"]],
B=snprmaRes[["B"]],
SNR=snprmaRes[["SNR"]],
mixtureParams=snprmaRes[["mixtureParams"]],
cdfName=cdfName,
row.names=NULL, ##featureNames(callSet),##[snp.index],
col.names=sampleNames(callSet),
probs=probs,
DF=DF,
SNRMin=SNRMin,
recallMin=recallMin,
recallRegMin=recallRegMin,
gender=gender,
verbose=verbose,
returnParams=returnParams,
badSNP=badSNP)
if(verbose) message("Leaving crlmmGT2")
open(tmp[["calls"]])
open(tmp[["confs"]])
##bb = ocProbesets()*ncol(A)*8
GT <- snpCall(callSet)
GTP <- snpCallProbability(callSet)
for(j in 1:ncol(callSet)){
GT[snp.index, j] <- tmp[["calls"]][, j]
GTP[snp.index, j] <- tmp[["confs"]][, j]
}
## ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb)
## ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb)
callSet$gender <- tmp$gender
return(callSet)
}
genotype3 <- function(filenames,
cdfName,
batch,
mixtureSampleSize=10^5,
eps=0.1,
verbose=TRUE,
seed=1,
sns,
copynumber=FALSE,
probs=rep(1/3, 3),
DF=6,
SNRMin=5,
recallMin=10,
recallRegMin=1000,
gender=NULL,
returnParams=TRUE,
badSNP=0.7){
if(!isPackageLoaded("ff")) stop("Must load package 'ff'")
if(!copynumber){
callSet <- crlmm2(filenames=filenames,
cdfName=cdfName,
mixtureSampleSize=mixtureSampleSize,
eps=eps,
verbose=verbose,
sns=sns,
probs=probs,
DF=DF,
SNRMin=SNRMin,
recallMin=recallMin,
recallRegMin=recallRegMin,
gender=gender,
returnParams=returnParams,
badSNP=badSNP)
return(callSet)
}
if(missing(cdfName)) stop("must specify cdfName")
if(!isValidCdfName(cdfName)) stop("cdfName not valid. see validCdfNames")
if(missing(batch)){
warning("The batch variable is not specified. The scan date of the array will be used as a surrogate for batch. The batch variable does not affect the preprocessing or genotyping, but is important for copy number estimation.")
} else {
if(length(batch) != length(filenames))
stop("batch variable must be the same length as the filenames")
}
if(missing(sns)) sns <- basename(filenames)
## callSet contains potentially very big matrices
callSet <- construct(filenames=filenames,
cdfName=cdfName,
copynumber=TRUE,
sns=sns,
verbose=verbose)
if(missing(batch)){
protocolData(callSet)$batch <- as.numeric(as.factor(protocolData(callSet)$ScanDate))
}
if(!missing(batch)) protocolData(callSet)$batch <- batch
##lM(callSet) <- initializeParamObject(list(featureNames(callSet), unique(protocolData(callSet)$batch)))
mixtureParams <- matrix(NA, 4, length(filenames))
snp.index <- which(isSnp(callSet)==1)
## snprmaRes <- snprma2(filenames=filenames,
## mixtureSampleSize=mixtureSampleSize,
## fitMixture=TRUE,
## eps=eps,
## verbose=verbose,
## seed=seed,
## cdfName=cdfName,
## sns=sns)
## if(verbose) message("Finished preprocessing.")
## open(snprmaRes[["A"]])
## open(snprmaRes[["B"]])
## open(snprmaRes[["SNR"]])
## open(snprmaRes[["SKW"]])
## open(snprmaRes[["mixtureParams"]])
## if(verbose) message("Updating elements of callSet")
#### bb = ocProbesets()*ncol(A)*8
#### ffrowapply(A(callSet)[i1:i2, ] <- snprmaRes[["A"]][i1:i2, ], X=snprmaRes[["A"]], BATCHBYTES=bb)
#### ffrowapply(B(callSet)[i1:i2, ] <- snprmaRes[["B"]][i1:i2, ], X=snprmaRes[["B"]], BATCHBYTES=bb)
## ##batches <- splitIndicesByLength(1:nrow(snprmaRes[["A"]]), ocProbesets())
## for(j in 1:ncol(callSet)){
## A(callSet)[snp.index, j] <- snprmaRes[["A"]][, j]
## B(callSet)[snp.index, j] <- snprmaRes[["B"]][, j]
## }
## if(verbose) message("Finished updating elements of callSet")
## stopifnot(identical(featureNames(callSet)[snp.index], snprmaRes$gns))
## pData(callSet)$SKW <- snprmaRes$SKW
## pData(callSet)$SNR <- snprmaRes$SNR
## mixtureParams <- snprmaRes$mixtureParams
np.index <- which(isSnp(callSet) == 0)
cnrmaRes <- cnrma2(A=A(callSet),
filenames=filenames,
row.names=featureNames(callSet)[np.index],
cdfName=cdfName,
sns=sns,
seed=seed,
verbose=verbose)
## if(verbose) message("Entering crlmmGT2...")
## rm(cnrmaRes); gc()
## ## as.matrix needed when ffdf is used
## tmp <- crlmmGT2(A=snprmaRes[["A"]],
## B=snprmaRes[["B"]],
## SNR=snprmaRes[["SNR"]],
## mixtureParams=snprmaRes[["mixtureParams"]],
## cdfName=cdfName,
## row.names=NULL, ##featureNames(callSet),##[snp.index],
## col.names=sampleNames(callSet),
## probs=probs,
## DF=DF,
## SNRMin=SNRMin,
## recallMin=recallMin,
## recallRegMin=recallRegMin,
## gender=gender,
## verbose=verbose,
## returnParams=returnParams,
## badSNP=badSNP)
## if(verbose) message("Leaving crlmmGT2")
## open(tmp[["calls"]])
## open(tmp[["confs"]])
## ##bb = ocProbesets()*ncol(A)*8
## for(j in 1:ncol(callSet)){
## snpCall(callSet)[snp.index, j] <- tmp[["calls"]][, j]
## snpCallProbability(callSet)[snp.index, j] <- tmp[["confs"]][, j]
## }
#### ffrowapply(snpCall(callSet)[i1:i2, ] <- tmp[["calls"]][i1:i2, ], X=tmp[["calls"]], BATCHBYTES=bb)
#### ffrowapply(snpCallProbability(callSet)[i1:i2, ] <- tmp[["confs"]][i1:i2, ], X=tmp[["confs"]], BATCHBYTES=bb)
## callSet$gender <- tmp$gender
#### cnSet <- as(callSet, "CNSetLM")
## return(callSet)
return(cnrmaRes)
}
##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
## For Illumina
##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
getPhenoData <- function(sampleSheet=NULL, arrayNames=NULL,
arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A")){
if(!is.null(arrayNames)) {
pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))
}
if(!is.null(sampleSheet)) { # get array info from Illumina's sample sheet
if(is.null(arrayNames)){
##arrayNames=NULL
if(!is.null(arrayInfoColNames$barcode) && (arrayInfoColNames$barcode %in% colnames(sampleSheet))) {
barcode = sampleSheet[,arrayInfoColNames$barcode]
arrayNames=barcode
}
if(!is.null(arrayInfoColNames$position) && (arrayInfoColNames$position %in% colnames(sampleSheet))) {
position = sampleSheet[,arrayInfoColNames$position]
if(is.null(arrayNames))
arrayNames=position
else
arrayNames = paste(arrayNames, position, sep=sep)
if(highDensity) {
hdExt = list(A="R01C01", B="R01C02", C="R02C01", D="R02C02")
for(i in names(hdExt))
arrayNames = sub(paste(sep, i, sep=""), paste(sep, hdExt[[i]], sep=""), arrayNames)
}
}
}
pd = new("AnnotatedDataFrame", data = sampleSheet)
sampleNames(pd) <- basename(arrayNames)
}
if(is.null(arrayNames)) {
arrayNames = gsub(paste(sep, fileExt$green, sep=""), "", dir(pattern=fileExt$green, path=path))
if(!is.null(sampleSheet)) {
sampleSheet=NULL
cat("Could not find required info in \'sampleSheet\' - ignoring. Check \'sampleSheet\' and/or \'arrayInfoColNames\'\n")
}
pd = new("AnnotatedDataFrame", data = data.frame(Sample_ID=arrayNames))
}
return(pd)
}
constructRG <- function(filenames, cdfName, sns, verbose, fileExt, sep, sampleSheet, arrayInfoColNames){
if(verbose) message("reading first idat file to extract feature data")
grnfile <- paste(filenames[1], fileExt$green, sep=sep)
if(!file.exists(grnfile)){
stop(paste(grnfile, " does not exist. Check fileExt argument"))
}
G <- readIDAT(grnfile)
idsG = rownames(G$Quants)
nr <- length(idsG)
fD <- new("AnnotatedDataFrame", data=data.frame(row.names=idsG))##, varMetadata=data.frame(labelDescript
nr <- nrow(fD)
dns <- list(featureNames(fD), basename(filenames))
RG <- new("NChannelSet",
R=initializeBigMatrix(name="R", nr=nr, nc=length(filenames)),
G=initializeBigMatrix(name="G", nr=nr, nc=length(filenames)),
zero=initializeBigMatrix(name="zero", nr=nr, nc=length(filenames)),
featureData=fD,
annotation=cdfName)
phenoData(RG) <- getPhenoData(sampleSheet=sampleSheet, arrayNames=filenames,
arrayInfoColNames=arrayInfoColNames)
## pD <- data.frame(matrix(NA, length(sampleNames(RG)), 12), row.names=sampleNames(RG))
## colnames(pD) <- c("Index","HapMap.Name","Name","ID",
## "Gender", "Plate", "Well", "Group", "Parent1",
## "Parent2","Replicate","SentrixPosition")
##phenoData(RG) <- new("AnnotatedDataFrame", data=pD)
pD <- data.frame(matrix(NA, length(sampleNames(RG)), 1), row.names=sampleNames(RG))
colnames(pD) <- "ScanDate"
protocolData(RG) <- new("AnnotatedDataFrame", data=pD)
sampleNames(RG) <- basename(filenames)
storageMode(RG) <- "environment"
RG##featureData=ops$illuminaOpts[["featureData"]])
}
crlmmIlluminaRS <- function(sampleSheet=NULL,
arrayNames=NULL,
batch,
ids=NULL,
path=".",
arrayInfoColNames=list(barcode="SentrixBarcode_A", position="SentrixPosition_A"),
highDensity=FALSE,
sep="_",
fileExt=list(green="Grn.idat", red="Red.idat"),
stripNorm=TRUE,
useTarget=TRUE,
row.names=TRUE,
col.names=TRUE,
probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,
seed=1, save.ab=FALSE, snpFile, cnFile,
mixtureSampleSize=10^5, eps=0.1, verbose=TRUE,
cdfName, sns, recallMin=10, recallRegMin=1000,
returnParams=FALSE, badSNP=.7,
copynumber=FALSE,
load.it=TRUE) {
if(missing(cdfName)) stop("must specify cdfName")
if(!isValidCdfName(cdfName)) stop("cdfName not valid. see validCdfNames")
if(missing(sns)) sns <- basename(arrayNames)
if(missing(batch)){
warning("The batch variable is not specified. The scan date of the array will be used as a surrogate for batch. The batch variable does not affect the preprocessing or genotyping, but is important for copy number estimation.")
} else {
if(length(batch) != length(sns))
stop("batch variable must be the same length as the filenames")
}
batches <- splitIndicesByLength(seq(along=arrayNames), ocSamples())
k <- 1
for(j in batches){
if(verbose) message("Batch ", k, " of ", length(batches))
RG <- readIdatFiles(sampleSheet=sampleSheet[j, ],
arrayNames=arrayNames[j],
ids=ids,
path=path,
arrayInfoColNames=arrayInfoColNames,
highDensity=highDensity,
sep=sep,
fileExt=fileExt,
saveDate=TRUE)
RG <- RGtoXY(RG, chipType=cdfName)
protocolData <- protocolData(RG)
res <- preprocessInfinium2(RG,
mixtureSampleSize=mixtureSampleSize,
fitMixture=TRUE,
verbose=verbose,
seed=seed,
eps=eps,
cdfName=cdfName,
sns=sns[j],
stripNorm=stripNorm,
useTarget=useTarget)
rm(RG); gc()
## MR: number of rows should be number of SNPs + number of nonpolymorphic markers.
## Here, I'm just using the # of rows returned from the above function
if(k == 1){
if(verbose) message("Initializing container for alleleA, alleleB, call, callProbability")
callSet <- new("SnpSuperSet",
alleleA=initializeBigMatrix(name="A", nr=nrow(res[[1]]), nc=length(sns)),
alleleB=initializeBigMatrix(name="B", nr=nrow(res[[1]]), nc=length(sns)),
call=initializeBigMatrix(name="call", nr=nrow(res[[1]]), nc=length(sns)),
callProbability=initializeBigMatrix(name="callPr", nr=nrow(res[[1]]), nc=length(sns)),
annotation=cdfName)
sampleNames(callSet) <- sns
phenoData(callSet) <- getPhenoData(sampleSheet=sampleSheet,
arrayNames=sns,
arrayInfoColNames=arrayInfoColNames)
pD <- data.frame(matrix(NA, length(sns), 1), row.names=sns)
colnames(pD) <- "ScanDate"
protocolData(callSet) <- new("AnnotatedDataFrame", data=pD)
pData(protocolData(callSet))[j, ] <- pData(protocolData)
featureNames(callSet) <- res[["gns"]]
pData(callSet)$SNR <- initializeBigVector("crlmmSNR-", length(sns), "double")
pData(callSet)$SKW <- initializeBigVector("crlmmSKW-", length(sns), "double")
##pData(callSet)$SKW <- rep(NA, length(sns))
##pData(callSet)$SNR <- rep(NA, length(sns))
pData(callSet)$gender <- rep(NA, length(sns))
mixtureParams <- initializeBigMatrix("crlmmMixt-", nr=4, nc=ncol(callSet), vmode="double")
save(mixtureParams, file=file.path(ldPath(), "mixtureParams.rda"))
if(missing(batch)){
protocolData(callSet)$batch <- rep(NA, length(sns))
} else{
protocolData(callSet)$batch <- batch
}
featureData(callSet) <- addFeatureAnnotation(callSet)
open(mixtureParams)
open(callSet$SNR)
open(callSet$SKW)
}
if(k > 1 & nrow(res[[1]]) != nrow(callSet)){
##RS: I don't understand why the IDATS for the
##same platform potentially have different lengths
res[["A"]] <- res[["A"]][res$gns %in% featureNames(callSet), ]
res[["B"]] <- res[["B"]][res$gns %in% featureNames(callSet), ]
}
if(missing(batch)){
protocolData(callSet)$batch[j] <- as.numeric(as.factor(protocolData$ScanDate))
}
## MR: we need to define a snp.index vs np.index
snp.index <- match(res$gns, featureNames(callSet))
A(callSet)[snp.index, j] <- res[["A"]]
B(callSet)[snp.index, j] <- res[["B"]]
pData(callSet)$SKW[j] <- res$SKW
pData(callSet)$SNR[j] <- res$SNR
mixtureParams[, j] <- res$mixtureParams
rm(res); gc()
k <- k+1
}
save(callSet, file=file.path(ldPath(), "callSet.rda"))
##otherwise, A and B get overwritten
##AA <- initializeBigMatrix("crlmmA", nrow(callSet), ncol(callSet), "integer")
##BB <- initializeBigMatrix("crlmmB", nrow(callSet), ncol(callSet), "integer")
##bb = ocProbesets()*ncol(A)*8
AA <- clone(A(callSet))
BB <- clone(B(callSet))
##ffrowapply(AA[i1:i2, ] <- A(callSet)[i1:i2, ], X=A(callSet), BATCHBYTES=bb)
##ffrowapply(BB[i1:i2, ] <- B(callSet)[i1:i2, ], X=B(callSet), BATCHBYTES=bb)
##crlmmGT2 overwrites A and B.
tmp <- crlmmGT2(A=A(callSet),
B=B(callSet),
SNR=callSet$SNR,
mixtureParams=mixtureParams,
cdfName=annotation(callSet),
row.names=featureNames(callSet),
col.names=sampleNames(callSet),
probs=probs,
DF=DF,
SNRMin=SNRMin,
recallMin=recallMin,
recallRegMin=recallRegMin,
gender=gender,
verbose=verbose,
returnParams=returnParams,
badSNP=badSNP)
open(tmp[["calls"]])
open(tmp[["confs"]])
A(callSet) <- AA
B(callSet) <- BB
snpCall(callSet) <- tmp[["calls"]]
## MR: many zeros in the conf. scores (?)
snpCallProbability(callSet) <- tmp[["confs"]]
callSet$gender <- tmp$gender
if(copynumber) cnSet <- as(callSet, "CNSetLM")
close(mixtureParams)
rm(tmp); gc()
return(cnSet)
}
##---------------------------------------------------------------------------
##---------------------------------------------------------------------------
rowCovs <- function(x, y, ...){
notna <- !is.na(x)
N <- rowSums(notna)
x <- suppressWarnings(log2(x))
if(!missing(y)) y <- suppressWarnings(log2(y))
return(rowSums((x - rowMeans(x, ...)) * (y - rowMeans(y, ...)), ...)/(N-1))
}
rowMAD <- function(x, y, ...){
notna <- !is.na(x)
mad <- 1.4826*rowMedians(abs(x-rowMedians(x, ...)), ...)
return(mad)
}
shrink <- function(x, Ns, DF.PRIOR){
DF <- Ns-1
DF[DF < 1] <- 1
x.0 <- apply(x, 2, median, na.rm=TRUE)
x <- (x*DF + x.0*DF.PRIOR)/(DF.PRIOR + DF)
for(j in 1:ncol(x)) x[is.na(x[, j]), j] <- x.0[j]
return(x)
}
applyByGenotype <- function(x, FUN, G){
FUN <- match.fun(FUN)
tmp <- matrix(NA, nrow(x), 3)
for(j in 1:3){
GT <- G==j
GT[GT == FALSE] <- NA
gt.x <- GT*x
tmp[, j] <- FUN(gt.x, na.rm=TRUE)
}
tmp
}
rowCors <- function(x, y, ...){
N <- rowSums(!is.na(x))
x <- suppressWarnings(log2(x))
y <- suppressWarnings(log2(y))
sd.x <- rowSds(x, ...)
sd.y <- rowSds(y, ...)
covar <- rowSums((x - rowMeans(x, ...)) * (y - rowMeans(y, ...)), ...)/(N-1)
return(covar/(sd.x*sd.y))
}
corByGenotype <- function(A, B, G, Ns, which.cluster=c(1,2,3)[1], DF.PRIOR){
x <- A * (G == which.cluster)
x[x==0] <- NA
y <- B * (G == which.cluster)
res <- as.matrix(rowCors(x, y, na.rm=TRUE))
cors <- shrink(res, Ns[, which.cluster], DF.PRIOR)
cors
}
generateX <- function(w, X) as.numeric(diag(w) %*% X)
generateIXTX <- function(x, nrow=3) {
X <- matrix(x, nrow=nrow)
XTX <- crossprod(X)
solve(XTX)
}
nuphiAllele2 <- function(allele, Ystar, W, Ns){
complete <- rowSums(is.na(W)) == 0
if(any(!is.finite(W))){## | any(!is.finite(V))){
i <- which(rowSums(!is.finite(W)) > 0)
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
}
NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
if(missing(allele)) stop("must specify allele")
if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
##How to quickly generate Xstar, Xstar = diag(W) %*% X
Xstar <- apply(W, 1, generateX, X)
IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
betahat <- matrix(NA, 2, nrow(Ystar))
ses <- matrix(NA, 2, nrow(Ystar))
for(i in 1:nrow(Ystar)){
betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
}
nu <- betahat[1, ]
phi <- betahat[2, ]
return(list(nu, phi))
}
nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
complete <- rowSums(is.na(W)) == 0
if(any(!is.finite(W))){## | any(!is.finite(V))){
i <- which(rowSums(!is.finite(W)) > 0)
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
}
##NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
if(missing(allele)) stop("must specify allele")
if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
##How to quickly generate Xstar, Xstar = diag(W) %*% X
Xstar <- apply(W, 1, generateX, X)
IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
betahat <- matrix(NA, 3, nrow(Ystar))
ses <- matrix(NA, 3, nrow(Ystar))
##betahat <- matrix(NA, 2, nrow(Ystar))
##ses <- matrix(NA, 2, nrow(Ystar))
for(i in 1:nrow(Ystar)){
betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
}
nu <- betahat[1, ]
phi <- betahat[2, ]
phi2 <- betahat[3, ]
return(list(nu, phi, phi2))
}
nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
I <- isSnp(object)
Ystar <- Ystar[I, ]
rownames(Ystar) <- featureNames(object)[isSnp(object)]
complete <- rowSums(is.na(W)) == 0 & I
W <- W[I, ]
if(any(!is.finite(W))){## | any(!is.finite(V))){
i <- which(rowSums(!is.finite(W)) > 0)
##browser()
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
}
Ns <- tmp.objects[["Ns"]]
Ns <- Ns[I, ]
CHR <- unique(chromosome(object))
##batch <- unique(object$batch)
batch <- unique(batch(object))
nuA <- getParam(object, "nuA", batch)
nuB <- getParam(object, "nuB", batch)
nuA.se <- getParam(object, "nuA.se", batch)
nuB.se <- getParam(object, "nuB.se", batch)
phiA <- getParam(object, "phiA", batch)
phiB <- getParam(object, "phiB", batch)
phiA.se <- getParam(object, "phiA.se", batch)
phiB.se <- getParam(object, "phiB.se", batch)
if(CHR==23){
phiAX <- getParam(object, "phiAX", batch)
phiBX <- getParam(object, "phiBX", batch)
}
NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
if(missing(allele)) stop("must specify allele")
if(CHR == 23){
##Design matrix for X chromosome depends on whether there was a sufficient number of AB genotypes
if(length(grep("AB", colnames(W))) > 0){
if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
} else{
if(allele == "A") X <- cbind(1, c(1, 0, 2, 0), c(0, 1, 0, 2))
if(allele == "B") X <- cbind(1, c(0, 1, 0, 2), c(1, 0, 2, 0))
}
} else {##autosome
if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
}
##How to quickly generate Xstar, Xstar = diag(W) %*% X
Xstar <- apply(W, 1, generateX, X)
IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
##as.numeric(diag(W[1, ]) %*% X)
if(CHR == 23){
betahat <- matrix(NA, 3, nrow(Ystar))
ses <- matrix(NA, 3, nrow(Ystar))
} else{
betahat <- matrix(NA, 2, nrow(Ystar))
ses <- matrix(NA, 2, nrow(Ystar))
}
for(i in 1:nrow(Ystar)){
betahat[, i] <- crossprod(matrix(IXTX[, i], ncol(X), ncol(X)), crossprod(matrix(Xstar[, i], nrow=nrow(X)), Ystar[i, ]))
ssr <- sum((Ystar[i, ] - matrix(Xstar[, i], nrow(X), ncol(X)) %*% matrix(betahat[, i], ncol(X), 1))^2)
ses[, i] <- sqrt(diag(matrix(IXTX[, i], ncol(X), ncol(X)) * ssr))
}
if(allele == "A"){
nuA[complete] <- betahat[1, ]
phiA[complete] <- betahat[2, ]
nuA.se[complete] <- ses[1, ]
phiA.se[complete] <- ses[2, ]
object <- pr(object, "nuA", batch, nuA)
object <- pr(object, "nuA.se", batch, nuA.se)
object <- pr(object, "phiA", batch, phiA)
object <- pr(object, "phiA.se", batch, phiA.se)
if(CHR == 23){
phiAX[complete] <- betahat[3, ]
object <- pr(object, "phiAX", batch, phiAX)
}
}
if(allele=="B"){
nuB[complete] <- betahat[1, ]
phiB[complete] <- betahat[2, ]
nuB.se[complete] <- ses[1, ]
phiB.se[complete] <- ses[2, ]
object <- pr(object, "nuB", batch, nuB)
object <- pr(object, "nuB.se", batch, nuB.se)
object <- pr(object, "phiB", batch, phiB)
object <- pr(object, "phiB.se", batch, phiB.se)
if(CHR == 23){
phiBX[complete] <- betahat[3, ]
object <- pr(object, "phiBX", batch, phiBX)
}
}
## THR.NU.PHI <- cnOptions$THR.NU.PHI
## if(THR.NU.PHI){
## verbose <- cnOptions$verbose
## if(verbose) message("Thresholding nu and phi")
## object <- thresholdModelParams(object, cnOptions)
## }
return(object)
}
predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){
pkgname <- paste(cdfName, "Crlmm", sep="")
path <- system.file("extdata", package=pkgname)
loader("23index.rda", .crlmmPkgEnv, pkgname)
index <- getVarInEnv("index")
##load(file.path(path, "23index.rda"))
XIndex <- index[[1]]
if(class(res) != "ABset"){
A <- res$A
B <- res$B
} else{
A <- res@assayData[["A"]]
B <- res@assayData[["B"]]
}
tmp <- which(rowSums(is.na(A)) > 0)
XMedian <- apply(log2(A[XIndex,, drop=FALSE])+log2(B[XIndex,, drop=FALSE]), 2, median, na.rm=TRUE)/2
SNR <- res$SNR
if(sum(SNR>SNRMin)==1){
gender <- which.min(c(abs(XMedian-8.9), abs(XMedian-9.5)))
}else{
gender <- kmeans(XMedian, c(min(XMedian[SNR>SNRMin]), max(XMedian[SNR>SNRMin])))[["cluster"]]
}
return(gender)
}
##replace with release/crlmm/R/cnrma-functions
crlmmCopynumber <- function(object,
chromosome=1:23,
which.batches,
MIN.SAMPLES=10,
SNRMin=5,
MIN.OBS=3,
DF.PRIOR=50,
bias.adj=FALSE,
prior.prob=rep(1/4,4),
seed=1,
verbose=TRUE,
GT.CONF.THR=0.99,
PHI.THR=2^6,
nHOM.THR=5,
MIN.NU=2^3,
MIN.PHI=2^3,
THR.NU.PHI=TRUE,
thresholdCopynumber=TRUE){
ffIsLoaded <- class(calls(object))[[1]] == "ff"
if(ffIsLoaded){
open(object)
open(object$SKW)
open(object$SNR)
}
stopifnot("batch" %in% varLabels(protocolData(object)))
stopifnot("chromosome" %in% fvarLabels(object))
stopifnot("position" %in% fvarLabels(object))
stopifnot("isSnp" %in% fvarLabels(object))
##batch <- object$batch
batch <- batch(object)
if(ffIsLoaded) {
open(object$SNR)
SNR <- object$SNR[, ]
} else SNR <- object$SNR
batches <- split((1:ncol(object))[SNR > SNRMin], batch[SNR > SNRMin])
if(any(sapply(batches, length) < MIN.SAMPLES)) message("Excluding batches with fewer than ", MIN.SAMPLES, " samples")
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
if(missing(which.batches)) which.batches <- seq(along=batches)
for(i in chromosome){
if(verbose) cat("Chromosome ", i, "\n")
if(i >= 24) next()
ii <- which.batches[1]
for(j in batches[which.batches]){
if(verbose) message("Batch ", ii, " of ", length(which.batches))
row.index <- which(chromosome(object) == i)
##Note that ffdf assayDataElements are data.frames after subsetting(not matrices)
ca <- as.matrix(CA(object)[row.index, j])
cb <- as.matrix(CB(object)[row.index, j])
dimnames(ca) <- dimnames(cb) <- list(featureNames(object)[row.index], sampleNames(object)[j])
tmp <- new("CNSet",
call=as.matrix(calls(object)[row.index, j]),
callProbability=as.matrix(snpCallProbability(object)[row.index, j]),
alleleA=as.matrix(A(object)[row.index, j]),
alleleB=as.matrix(B(object)[row.index, j]),
CA=ca, CB=cb,
phenoData=phenoData(object)[j, ],
annotation=annotation(object))
featureData(tmp) <- addFeatureAnnotation(tmp)
featureData(tmp) <- lm.parameters(tmp, batch=unique(batch[j]))
tmp$batch <- batch(object)[j]
tmp <- computeCopynumber(tmp,
MIN.OBS=MIN.OBS,
DF.PRIOR=DF.PRIOR,
bias.adj=bias.adj,
prior.prob=prior.prob,
seed=seed,
verbose=verbose,
GT.CONF.THR=GT.CONF.THR,
PHI.THR=PHI.THR,
nHOM.THR=nHOM.THR,
MIN.NU=MIN.NU,
MIN.PHI=MIN.PHI,
THR.NU.PHI=THR.NU.PHI,
thresholdCopynumber=thresholdCopynumber)
fData(tmp) <- fData(tmp)[, -(1:3)]
CA(tmp) <- matrix(as.integer(CA(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
dimnames=list(featureNames(tmp), sampleNames(tmp)))
CB(tmp) <- matrix(as.integer(CB(tmp)*100), nrow=nrow(tmp), ncol=ncol(tmp),
dimnames=list(featureNames(tmp), sampleNames(tmp)))
CA(object)[row.index, j] <- CA(tmp)
CB(object)[row.index, j] <- CB(tmp)
##ad-hocery
batchName <- unique(batch(object)[j])
fvarLabels(tmp)[15:17] <- paste(c("corrAB", "corrBB", "corrAA"), batchName, sep=".")
fvarLabels(tmp)[13:14] <- paste(c("phiPrimeA", "phiPrimeB"), batchName, sep=".")
fvarLabels(tmp) <- gsub("_", ".", fvarLabels(tmp))
##fvarLabels(tmp) <- gsub("\\.[1-9]", "", fvarLabels(tmp))
if(ffIsLoaded){
fData(tmp) <- fData(tmp)[, fvarLabels(tmp) %in% names(physical(lM(object)))]
jj <- match(fvarLabels(tmp), names(lM(object)))
lM(object)[row.index, jj] <- fData(tmp)
} else {
nms <- paste(names(lM(object)), batchName, sep=".")
fData(tmp) <- fData(tmp)[, fvarLabels(tmp) %in% nms]
for(k in seq_along(fvarLabels(tmp))){
kk <- match(fvarLabels(tmp)[k], paste(names(lM(object)), batchName, sep="."))
column <- match(batchName, colnames(lM(object)[[k]]))
lM(object)[[k]][row.index, column] <- fData(tmp)[, k]
}
}
rm(tmp); gc()
ii <- ii+1
}
}
return(object)
}
crlmmCopynumber2 <- function(object,
which.batches,
MIN.SAMPLES=10,
SNRMin=5,
MIN.OBS=1,
DF.PRIOR=50,
bias.adj=FALSE,
prior.prob=rep(1/4,4),
seed=1,
verbose=TRUE,
GT.CONF.THR=0.99,
PHI.THR=2^6,
nHOM.THR=5,
MIN.NU=2^3,
MIN.PHI=2^3,
THR.NU.PHI=TRUE,
thresholdCopynumber=TRUE){
stopifnot("batch" %in% varLabels(protocolData(object)))
stopifnot("chromosome" %in% fvarLabels(object))
stopifnot("position" %in% fvarLabels(object))
stopifnot("isSnp" %in% fvarLabels(object))
##tmp <- initializeParamObject(list(featureNames(object), unique(protocolData(object)$batch)))
##lM(object) <- tmp
object@lM <- initializeParamObject(list(featureNames(object), unique(protocolData(object)$batch)))
lM(object) <- initializeParamObject(list(featureNames(object), unique(protocolData(object)$batch)))
batch <- batch(object)
XIndex.snps <- which(chromosome(object) == 23 & isSnp(object) & !is.na(chromosome(object)))
##YIndex.snps <- (1:nrow(object))[chromosome(object) == 24 & isSnp(object)]
XIndex.nps <- which(chromosome(object) == 23 & !isSnp(object) & !is.na(chromosome(object)))
##autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]
autosomeIndex.snps <- which(chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object)))
autosomeIndex.nps <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))
##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]
## Do chromosome X in batches
Ns <- initializeBigMatrix("Ns", nrow(object), 5)
colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
if(!file.exists(file.path(ldPath(), "normal.rda"))){
normal <- initializeBigMatrix("normal", nrow(object), ncol(object), vmode="integer")
normal[,] <- 1L
save(normal, file=file.path(ldPath(), "normal.rda"))
} else load(file.path(ldPath(), "normal.rda"))
if(!file.exists(file.path(ldPath(), "snpflags.rda"))){
snpflags <- initializeBigMatrix("snpflags", nrow(object), length(unique(batch(object))), vmode="integer")
snpflags[,] <- 0L
save(snpflags, file=file.path(ldPath(), "snpflags.rda"))
} else{
load(file.path(ldPath(), "snpflags.rda"))
}
if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
snpBatches <- splitIndicesByLength(autosomeIndex.snps, ocProbesets())
ocLapply(seq(along=snpBatches),
fit.lm1,
autosomeIndex=autosomeIndex.snps,
object=object,
Ns=Ns,
normal=normal,
snpflags=snpflags,
snpBatches=snpBatches,
batchSize=ocProbesets(),
SNRMin=SNRMin,
MIN.SAMPLES=MIN.SAMPLES,
MIN.OBS=MIN.OBS,
DF=DF.PRIOR,
GT.CONF.THR=GT.CONF.THR,
THR.NU.PHI=THR.NU.PHI,
MIN.NU=MIN.NU,
MIN.PHI=MIN.PHI,
verbose=verbose,
neededPkgs="crlmm")
## autosomal NPs
snpBatches <- splitIndicesByLength(autosomeIndex.nps, ocProbesets())
if(verbose) message("Estimating total copy number at nonpolymorphic loci")
ocLapply(seq(along=snpBatches),
fit.lm2,
autosomeIndex=autosomeIndex.nps,
object=object,
Ns=Ns,
normal=normal,
snpflags=snpflags,
snpBatches=snpBatches,
batchSize=ocProbesets(),
SNRMin=SNRMin,
MIN.SAMPLES=MIN.SAMPLES,
MIN.OBS=MIN.OBS,
DF=DF.PRIOR,
GT.CONF.THR=GT.CONF.THR,
THR.NU.PHI=THR.NU.PHI,
MIN.NU=MIN.NU,
MIN.PHI=MIN.PHI,
verbose=verbose,
neededPkgs="crlmm")
snpBatches <- splitIndicesByLength(XIndex.snps, ocProbesets())
if(verbose) message("Estimating allele-specific copy number at polymorphic loci on chromosome X")
ocLapply(seq(along=snpBatches),
fit.lm3,
autosomeIndex=XIndex.snps,
object=object,
Ns=Ns,
normal=normal,
snpflags=snpflags,
snpBatches=snpBatches,
batchSize=ocProbesets(),
SNRMin=SNRMin,
MIN.SAMPLES=MIN.SAMPLES,
MIN.OBS=MIN.OBS,
DF=DF.PRIOR,
GT.CONF.THR=GT.CONF.THR,
THR.NU.PHI=THR.NU.PHI,
MIN.NU=MIN.NU,
MIN.PHI=MIN.PHI,
verbose=verbose,
neededPkgs="crlmm")
if(verbose) message("Estimating total copy number for nonpolymorphic loci on chromosome X")
snpBatches <- splitIndicesByLength(XIndex.nps, ocProbesets())
tmp <- ocLapply(seq(along=snpBatches),
fit.lm4,
XIndex=XIndex.nps,
object=object,
Ns=Ns,
normal=normal,
snpflags=snpflags,
snpBatches=snpBatches,
batchSize=ocProbesets(),
SNRMin=SNRMin,
MIN.SAMPLES=MIN.SAMPLES,
MIN.OBS=MIN.OBS,
DF=DF.PRIOR,
GT.CONF.THR=GT.CONF.THR,
THR.NU.PHI=THR.NU.PHI,
MIN.NU=MIN.NU,
MIN.PHI=MIN.PHI,
verbose=verbose,
neededPkgs="crlmm")
return(object)
}
fit.lm1 <- function(idxBatch,
snpBatches,
autosomeIndex,
object,
Ns,
normal,
snpflags,
batchSize,
SNRMin,
MIN.SAMPLES,
MIN.OBS,
DF.PRIOR,
GT.CONF.THR,
THR.NU.PHI,
MIN.NU,
MIN.PHI,
verbose, ...){
if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
snps <- snpBatches[[idxBatch]]
batches <- split(seq(along=batch(object)), batch(object))
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
open(object)
open(snpflags)
open(normal)
corrAB <- corrBB <- corrAA <- sig2B <- sig2A <- tau2B <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))
flags <- nuA <- nuB <- phiA <- phiB <- corrAB
normal.snps <- normal[snps, ]
cB <- cA <- matrix(NA, length(snps), ncol(object))
GG <- as.matrix(calls(object)[snps, ])
CP <- as.matrix(snpCallProbability(object)[snps, ])
AA <- as.matrix(A(object)[snps, ])
BB <- as.matrix(B(object)[snps, ])
for(k in batches){
G <- GG[, k]
NORM <- normal.snps[, k]
xx <- CP[, k]
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
G <- G*highConf*NORM
A <- AA[, k]
B <- BB[, k]
##index <- GT.B <- GT.A <- vector("list", 3)
##names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
Ns <- applyByGenotype(matrix(1, nrow(G), ncol(G)), rowSums, G)
muA <- applyByGenotype(A, rowMedians, G)
muB <- applyByGenotype(B, rowMedians, G)
vA <- applyByGenotype(A, rowMAD, G)
vB <- applyByGenotype(B, rowMAD, G)
vA <- shrink(vA, Ns, DF.PRIOR)
vB <- shrink(vB, Ns, DF.PRIOR)
##location and scale
J <- match(unique(batch(object)[k]), unique(batch(object)))
##background variance for alleleA
taus <- applyByGenotype(log2(A), rowMAD, G)^2
tau2A[, J] <- shrink(taus[, 3, drop=FALSE], Ns[, 3], DF.PRIOR)
sig2A[, J] <- shrink(taus[, 1, drop=FALSE], Ns[, 1], DF.PRIOR)
taus <- applyByGenotype(log2(B), rowMAD, G)^2
tau2B[, J] <- shrink(taus[, 3, drop=FALSE], Ns[, 1], DF.PRIOR)
sig2B[, J] <- shrink(taus[, 1, drop=FALSE], Ns[, 3], DF.PRIOR)
corrAB[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=2, DF.PRIOR)
corrAA[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=1, DF.PRIOR)
corrBB[, J] <- corByGenotype(A=A, B=B, G=G, Ns=Ns, which.cluster=3, DF.PRIOR)
##---------------------------------------------------------------------------
## Impute sufficient statistics for unobserved genotypes (plate-specific)
##---------------------------------------------------------------------------
index <- apply(Ns, 2, function(x, MIN.OBS) which(x > MIN.OBS), MIN.OBS)
correct.orderA <- muA[, 1] > muA[, 3]
correct.orderB <- muB[, 3] > muB[, 1]
index.complete <- intersect(which(correct.orderA & correct.orderB), intersect(index[[1]], intersect(index[[2]], index[[3]])))
size <- min(5000, length(index.complete))
if(size == 5000) index.complete <- sample(index.complete, 5000, replace=TRUE)
if(length(index.complete) < 200){
warning("fewer than 200 snps pass criteria for predicting the sufficient statistics")
return()
}
index <- vector("list", 3)
index[[1]] <- which(Ns[, 1] == 0 & (Ns[, 2] >= MIN.OBS & Ns[, 3] >= MIN.OBS))
index[[2]] <- which(Ns[, 2] == 0 & (Ns[, 1] >= MIN.OBS & Ns[, 3] >= MIN.OBS))
index[[3]] <- which(Ns[, 3] == 0 & (Ns[, 2] >= MIN.OBS & Ns[, 1] >= MIN.OBS))
res <- imputeCenter(muA, muB, index.complete, index)
muA <- res[[1]]
muB <- res[[2]]
## Monomorphic SNPs. Mixture model may be better
## Improve estimation by borrowing strength across batch
noAA <- Ns[, 1] < MIN.OBS
noAB <- Ns[, 2] < MIN.OBS
noBB <- Ns[, 3] < MIN.OBS
index[[1]] <- noAA & noAB
index[[2]] <- noBB & noAB
index[[3]] <- noAA & noBB
cols <- c(3, 1, 2)
for(j in 1:3){
if(sum(index[[j]]) == 0) next()
kk <- cols[j]
X <- cbind(1, muA[index.complete, kk], muB[index.complete, kk])
Y <- cbind(muA[index.complete, -kk],
muB[index.complete, -kk])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, muA[index[[j]], kk], muB[index[[j]], kk])
mus <- X %*% betahat
muA[index[[j]], -kk] <- mus[, 1:2]
muB[index[[j]], -kk] <- mus[, 3:4]
}
rm(betahat, X, Y, mus, index, noAA, noAB, noBB, res)
gc()
negA <- rowSums(muA < 0) > 0
negB <- rowSums(muB < 0) > 0
flags[, J] <- rowSums(Ns == 0) > 0
##flags[, J] <- index[[1]] | index[[2]] | index[[3]] | rowSums(
## we're regressing on the medians using the standard errors (hence the division by N) as weights
##formerly coefs()
Np <- Ns
Np[Np < 1] <- 1
vA2 <- vA^2/Np
vB2 <- vB^2/Np
wA <- sqrt(1/vA2)
wB <- sqrt(1/vB2)
YA <- muA*wA
YB <- muB*wB
res <- nuphiAllele2(allele="A", Ystar=YA, W=wA, Ns=Ns)
nuA[, J] <- res[[1]]
phiA[, J] <- res[[2]]
res <- nuphiAllele2(allele="B", Ystar=YB, W=wB, Ns=Ns)
nuB[, J] <- res[[1]]
phiB[, J] <- res[[2]]
if(THR.NU.PHI){
nuA[nuA[, J] < MIN.NU, J] <- MIN.NU
nuB[nuB[, J] < MIN.NU, J] <- MIN.NU
phiA[phiA[, J] < MIN.PHI, J] <- MIN.PHI
phiB[phiB[, J] < MIN.PHI, J] <- MIN.PHI
}
## formerly polymorphic(): calculate copy number
cA[, k] <- matrix((1/phiA[, J]*(A-nuA[, J])), nrow(A), ncol(A))
cB[, k] <- matrix((1/phiB[, J]*(B-nuB[, J])), nrow(B), ncol(B))
rm(G, A, B, NORM, wA, wB, YA,YB, res, negA, negB, Np, Ns)
gc()
}
cA[cA < 0.05] <- 0.05
cB[cB < 0.05] <- 0.05
cA[cA > 5] <- 5
cB[cB > 5] <- 5
cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
cB <- matrix(as.integer(cB*100), nrow(cB), ncol(cB))
CA(object)[snps, ] <- cA
CB(object)[snps, ] <- cB
snpflags[snps, ] <- flags
lapply(lM(object), open)
tmp <- physical(lM(object))$tau2A
tmp[snps, ] <- tau2A
lM(object)$tau2A <- tmp
tmp <- physical(lM(object))$tau2B
tmp[snps, ] <- tau2B
lM(object)$tau2B <- tmp
tmp <- physical(lM(object))$tau2B
tmp[snps, ] <- tau2B
lM(object)$tau2B <- tmp
tmp <- physical(lM(object))$sig2A
tmp[snps, ] <- sig2A
lM(object)$sig2A <- tmp
tmp <- physical(lM(object))$sig2B
tmp[snps, ] <- sig2B
lM(object)$sig2B <- tmp
tmp <- physical(lM(object))$nuA
tmp[snps, ] <- nuA
lM(object)$nuA <- tmp
tmp <- physical(lM(object))$nuB
tmp[snps, ] <- nuB
lM(object)$nuB <- tmp
tmp <- physical(lM(object))$phiA
tmp[snps, ] <- phiA
lM(object)$phiA <- tmp
tmp <- physical(lM(object))$phiB
tmp[snps, ] <- phiB
lM(object)$phiB <- tmp
tmp <- physical(lM(object))$corrAB
tmp[snps, ] <- corrAB
lM(object)$corrAB <- tmp
tmp <- physical(lM(object))$corrAA
tmp[snps, ] <- corrAA
lM(object)$corrAA <- tmp
tmp <- physical(lM(object))$corrBB
tmp[snps, ] <- corrBB
lM(object)$corrAB <- tmp
lapply(assayData(object), close)
lapply(lM(object), close)
TRUE
}
fit.lm2 <- function(idxBatch,
snpBatches,
autosomeIndex,
object,
Ns,
normal,
snpflags,
batchSize,
SNRMin,
MIN.SAMPLES,
MIN.OBS,
DF.PRIOR,
GT.CONF.THR,
THR.NU.PHI,
MIN.NU,
MIN.PHI,
verbose, ...){
## which.batches, ...){
if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
snps <- snpBatches[[idxBatch]]
batches <- split(seq(along=batch(object)), batch(object))
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
open(object)
open(snpflags)
open(normal)
cA <- matrix(NA, length(snps), ncol(object))
ii <- isSnp(object) & chromosome(object) < 23 & !is.na(chromosome(object))
flags <- as.matrix(snpflags[,])
noflags <- rowSums(flags, na.rm=TRUE) == 0 ##NA's for unevaluated batches
## We do not want to write to discuss for each batch. More efficient to
## write to disk after estimating these parameters for all batches.
nuA.np <- phiA.np <- sig2A.np <- matrix(NA, length(snps), length(unique(batch(object))))
## for imputation, we need the corresponding parameters of the snps
NN <- min(10e3, length(which(ii & noflags)))
snp.ind <- sample(which(ii & noflags), NN)
nnuA.snp <- as.matrix(physical(lM(object))$nuA[snp.ind,])
pphiA.snp <- as.matrix(physical(lM(object))$phiA[snp.ind,])
nnuB.snp <- as.matrix(physical(lM(object))$nuB[snp.ind,])
pphiB.snp <- as.matrix(physical(lM(object))$phiB[snp.ind,])
AA.snp <- as.matrix(A(object)[snp.ind, ])
BB.snp <- as.matrix(B(object)[snp.ind, ])
NNORM.snp <- as.matrix(normal[snp.ind, ])
NORM.np <- as.matrix(normal[snps, ])
AA.np <- as.matrix(A(object)[snps, ])
GG <- as.matrix(calls(object)[snp.ind, ])
CP <- as.matrix(snpCallProbability(object)[snp.ind, ])
for(k in batches){
##if(verbose) message("SNP batch ", ii, " of ", length(batches))
J <- match(unique(batch(object)[k]), unique(batch(object)))
## snp.index <- snp.ind & nuA[, J] > 20 & nuB[, J] > 20 & phiA[, J] > 20 & phiB[, J] > 20
## if(sum(snp.index) >= 5000){
## snp.index <- sample(which(snp.index), 5000)
## } else snp.index <- which(snp.index)
phiA.snp <- pphiA.snp[, J]
phiB.snp <- pphiB.snp[, J]
A.snp <- AA.snp[, k]
B.snp <- BB.snp[, k]
NORM.snp <- NNORM.snp[, k]
G <- GG[, k]
xx <- CP[, k]
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
G <- G*highConf*NORM.snp
G[G==0] <- NA
##nonpolymorphic
A.np <- AA.np[, k]
Ns <- applyByGenotype(matrix(1, nrow(G), ncol(G)), rowSums, G)
muA <- applyByGenotype(A.snp, rowMedians, G)
muB <- applyByGenotype(B.snp, rowMedians, G)
muA <- muA[, 1]
muB <- muB[, 3]
X <- cbind(1, log2(c(muA, muB)))
Y <- log2(c(phiA.snp, phiB.snp))
betahat <- solve(crossprod(X), crossprod(X, Y))
##
mus <- rowMedians(A.np * NORM.np[, k], na.rm=TRUE)
crosshyb <- max(median(muA) - median(mus), 0)
X <- cbind(1, log2(mus+crosshyb))
logPhiT <- X %*% betahat
phiA.np[, J] <- 2^(logPhiT)
nuA.np[, J] <- mus-2*phiA.np[, J]
if(THR.NU.PHI){
nuA.np[nuA.np[, J] < MIN.NU, J] <- MIN.NU
phiA.np[phiA.np[, J] < MIN.PHI, J] <- MIN.PHI
}
cA[, k] <- 1/phiA.np[, J] * (A.np - nuA.np[, J])
sig2A.np[, J] <- rowMAD(log2(A.np*NORM.np[, k]), na.rm=TRUE)
rm(NORM.snp, highConf, xx, G, Ns, A.np, X, Y, betahat, mus, logPhiT)
gc()
}
cA[cA < 0.05] <- 0.05
cA[cA > 5] <- 5
cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
CA(object)[snps, ] <- cA
##open(lM(object))
tmp <- physical(lM(object))$nuA
tmp[snps, ] <- nuA.np
lM(object)$nuA <- tmp
tmp <- physical(lM(object))$sig2A
tmp[snps, ] <- sig2A.np
lM(object)$sig2A <- tmp
tmp <- physical(lM(object))$phiA
tmp[snps, ] <- phiA.np
lM(object)$sig2A <- tmp
## lM(object)$sig2A[snps, ] <- sig2A.np
## lM(object)$nuA[snps, ] <- nuA.np
## lM(object)$phiA[snps, ] <- phiA.np
lapply(assayData(object), close)
lapply(lM(object), close)
TRUE
}
fit.lm3 <- function(idxBatch,
snpBatches,
XIndex,
object,
Ns,
normal,
snpflags,
batchSize,
SNRMin,
MIN.SAMPLES,
MIN.OBS,
DF.PRIOR,
GT.CONF.THR,
THR.NU.PHI,
MIN.NU,
MIN.PHI,
verbose, ...){
## which.batches, ...){
if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
snps <- snpBatches[[idxBatch]]
batches <- split(seq(along=batch(object)), batch(object))
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
open(snpflags)
open(normal)
open(object)
corrAB <- corrBB <- corrAA <- sig2B <- sig2A <- tau2B <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))
phiA2 <- phiB2 <- tau2A
flags <- nuA <- nuB <- phiA <- phiB <- corrAB
cB <- cA <- matrix(NA, length(snps), ncol(object))
gender <- object$gender
IX <- matrix(gender, length(snps), ncol(object))
NORM <- normal[snps,]
IX <- IX==2
GG <- as.matrix(calls(object)[snps, ])
CP <- as.matrix(snpCallProbability(object)[snps,])
AA <- as.matrix(A(object)[snps, ])
BB <- as.matrix(B(object)[snps, ])
for(k in batches){
##if(verbose) message("SNP batch ", ii, " of ", length(batches))
## within-genotype moments
gender <- object$gender[k]
G <- GG[, k]
xx <- CP[, k]
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
G <- G*highConf*NORM[, k]
A <- AA[, k]
B <- BB[, k]
##index <- GT.B <- GT.A <- vector("list", 3)
##names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
Ns.F <- applyByGenotype(matrix(1, nrow(G), sum(gender==2)), rowSums, G[, gender==2])
Ns.M <- applyByGenotype(matrix(1, nrow(G), sum(gender==1)), rowSums, G[, gender==1])
Ns <- cbind(Ns.M[, 1], Ns.M[, 3], Ns.F)
muA.F <- applyByGenotype(A[, gender==2], rowMedians, G[, gender==2])
muA.M <- applyByGenotype(A[, gender==1], rowMedians, G[, gender==1])
muB.F <- applyByGenotype(B[, gender==2], rowMedians, G[, gender==2])
muB.M <- applyByGenotype(B[, gender==1], rowMedians, G[, gender==1])
vA.F <- applyByGenotype(A[, gender==2], rowMAD, G[, gender==2])
vB.F <- applyByGenotype(B[, gender==2], rowMAD, G[, gender==2])
vA.M <- applyByGenotype(A[, gender==1], rowMAD, G[, gender==1])
vB.M <- applyByGenotype(B[, gender==1], rowMAD, G[, gender==1])
vA.F <- shrink(vA.F, Ns.F, DF.PRIOR)
vA.M <- shrink(vA.M, Ns.M, DF.PRIOR)
vB.F <- shrink(vB.F, Ns.F, DF.PRIOR)
vB.M <- shrink(vB.M, Ns.M, DF.PRIOR)
##location and scale
J <- match(unique(batch(object)[k]), unique(batch(object)))
##background variance for alleleA
taus <- applyByGenotype(log2(A[, gender==2]), rowMAD, G[, gender==2])^2
tau2A[, J] <- shrink(taus[, 3, drop=FALSE], Ns.F[, 3], DF.PRIOR)
sig2A[, J] <- shrink(taus[, 1, drop=FALSE], Ns.F[, 1], DF.PRIOR)
taus <- applyByGenotype(log2(B[, gender==2]), rowMAD, G[, gender==2])^2
tau2B[, J] <- shrink(taus[, 3, drop=FALSE], Ns.F[, 1], DF.PRIOR)
sig2B[, J] <- shrink(taus[, 1, drop=FALSE], Ns.F[, 3], DF.PRIOR)
corrAB[, J] <- corByGenotype(A=A[, gender==2], B=B[, gender==2], G=G[, gender==2], Ns=Ns.F, which.cluster=2, DF.PRIOR)
corrAA[, J] <- corByGenotype(A=A[, gender==2], B=B[, gender==2], G=G[, gender==2], Ns=Ns.F, which.cluster=1, DF.PRIOR)
corrBB[, J] <- corByGenotype(A=A[, gender==2], B=B[, gender==2], G=G[, gender==2], Ns=Ns.F, which.cluster=3, DF.PRIOR)
##formerly oneBatch()...
##---------------------------------------------------------------------------
## Impute sufficient statistics for unobserved genotypes (plate-specific)
##---------------------------------------------------------------------------
index <- apply(Ns.F, 2, function(x, MIN.OBS) which(x >= MIN.OBS), MIN.OBS)
correct.orderA <- muA.F[, 1] > muA.F[, 3]
correct.orderB <- muB.F[, 3] > muB.F[, 1]
index.complete <- intersect(which(correct.orderA & correct.orderB), intersect(index[[1]], intersect(index[[2]], index[[3]])))
size <- min(5000, length(index.complete))
if(length(index.complete) < 200){
warning("fewer than 200 snps pass criteria for predicting the sufficient statistics")
return()
}
if(size==5000) index.complete <- sample(index.complete, size)
index <- vector("list", 3)
index[[1]] <- which(Ns.F[, 1] == 0 & (Ns.F[, 2] >= MIN.OBS & Ns.F[, 3] >= MIN.OBS))
index[[2]] <- which(Ns.F[, 2] == 0 & (Ns.F[, 1] >= MIN.OBS & Ns.F[, 3] >= MIN.OBS))
index[[3]] <- which(Ns.F[, 3] == 0 & (Ns.F[, 2] >= MIN.OBS & Ns.F[, 1] >= MIN.OBS))
res <- imputeCenter(muA.F, muB.F, index.complete, index)
muA.F <- res[[1]]
muB.F <- res[[2]]
nobsA <- Ns.M[, 1] > MIN.OBS
nobsB <- Ns.M[, 3] > MIN.OBS
notMissing <- !(is.na(muA.M[, 1]) | is.na(muA.M[, 3]) | is.na(muB.M[, 1]) | is.na(muB.M[, 3]))
complete <- list()
complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
size <- min(5000, length(complete[[1]]))
if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
##
res <- imputeCenterX(muA.M, muB.M, Ns.M, complete, MIN.OBS)
muA.M <- res[[1]]
muB.M <- res[[2]]
##
## Monomorphic SNPs. Mixture model may be better
## Improve estimation by borrowing strength across batch
noAA <- Ns.F[, 1] < MIN.OBS
noAB <- Ns.F[, 2] < MIN.OBS
noBB <- Ns.F[, 3] < MIN.OBS
index[[1]] <- noAA & noAB
index[[2]] <- noBB & noAB
index[[3]] <- noAA & noBB
cols <- c(3, 1, 2)
for(j in 1:3){
if(sum(index[[j]]) == 0) next()
kk <- cols[j]
X <- cbind(1, muA.F[index.complete, kk], muB.F[index.complete, kk])
Y <- cbind(muA.F[index.complete, -kk],
muB.F[index.complete, -kk])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, muA.F[index[[j]], kk], muB.F[index[[j]], kk])
mus <- X %*% betahat
muA.F[index[[j]], -kk] <- mus[, 1:2]
muB.F[index[[j]], -kk] <- mus[, 3:4]
}
negA <- rowSums(muA.F < 0) > 0
negB <- rowSums(muB.F < 0) > 0
flags[, J] <- rowSums(Ns.F == 0) > 0 | negA | negB
##flags[, J] <- index[[1]] | index[[2]] | index[[3]] | rowSums(
##formerly coefs()
Np <- cbind(Ns.M[, c(1,3)], Ns.F)
Np[Np < 1] <- 1
vA <- cbind(vA.M[, c(1, 3)], vA.F)
vB <- cbind(vB.M[, c(1, 3)], vB.F)
muA <- cbind(muA.M[, c(1,3)], muA.F)
muB <- cbind(muB.M[, c(1,3)], muB.F)
vA2 <- vA^2/Np
vB2 <- vB^2/Np
wA <- sqrt(1/vA2)
wB <- sqrt(1/vB2)
YA <- muA*wA
YB <- muB*wB
res <- nuphiAlleleX(allele="A", Ystar=YA, W=wA)
nuA[, J] <- res[[1]]
phiA[, J] <- res[[2]]
phiA2[, J] <- res[[3]]
res <- nuphiAlleleX(allele="B", Ystar=YB, W=wB)
nuB[, J] <- res[[1]]
phiB[, J] <- res[[2]]
phiB2[, J] <- res[[3]]
if(THR.NU.PHI){
nuA[nuA[, J] < MIN.NU, J] <- MIN.NU
nuB[nuB[, J] < MIN.NU, J] <- MIN.NU
phiA[phiA[, J] < MIN.PHI, J] <- MIN.PHI
phiA2[phiA2[, J] < MIN.PHI, J] <- MIN.PHI
phiB[phiB[, J] < MIN.PHI, J] <- MIN.PHI
phiB2[phiB2[, J] < MIN.PHI, J] <- MIN.PHI
}
phistar <- phiB2[, J]/phiA[, J]
tmp <- (B-nuB[, J] - phistar*A + phistar*nuA[, J])/phiB[, J]
cB[, k] <- tmp/(1-phistar*phiA2[, J]/phiB[, J])
cA[, k] <- (A-nuA[, J]-phiA2[, J]*cB[, k])/phiA[, J]
##some of the snps are called for the men, but not the women
rm(YA, YB, wA, wB, res, tmp, phistar, A, B, G, index)
gc()
}
cA[cA < 0.05] <- 0.05
cB[cB < 0.05] <- 0.05
cA[cA > 5] <- 5
cB[cB > 5] <- 5
##--------------------------------------------------
##RS: need to fix. why are there NA's by coercion
cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
##--------------------------------------------------
##ii <- rowSums(is.na(cA)) > 0
##these often arise at SNPs with low confidence scores
cB <- matrix(as.integer(cB*100), nrow(cB), ncol(cB))
CA(object)[snps, ] <- cA
CB(object)[snps, ] <- cB
snpflags[snps, ] <- flags
tmp <- physical(lM(object))$tau2A
tmp[snps, ] <- tau2A
lM(object)$tau2A <- tmp
tmp <- physical(lM(object))$tau2B
tmp[snps, ] <- tau2B
lM(object)$tau2B <- tmp
tmp <- physical(lM(object))$tau2B
tmp[snps, ] <- tau2B
lM(object)$tau2B <- tmp
tmp <- physical(lM(object))$sig2A
tmp[snps, ] <- sig2A
lM(object)$sig2A <- tmp
tmp <- physical(lM(object))$sig2B
tmp[snps, ] <- sig2B
lM(object)$sig2B <- tmp
tmp <- physical(lM(object))$nuA
tmp[snps, ] <- nuA
lM(object)$nuA <- tmp
tmp <- physical(lM(object))$nuB
tmp[snps, ] <- nuB
lM(object)$nuB <- tmp
tmp <- physical(lM(object))$phiA
tmp[snps, ] <- phiA
lM(object)$phiA <- tmp
tmp <- physical(lM(object))$phiB
tmp[snps, ] <- phiB
lM(object)$phiB <- tmp
tmp <- physical(lM(object))$phiPrimeA
tmp[snps, ] <- phiA2
lM(object)$phiPrimeA <- tmp
tmp <- physical(lM(object))$phiPrimeB
tmp[snps, ] <- phiB2
lM(object)$phiPrimeB <- tmp
tmp <- physical(lM(object))$corrAB
tmp[snps, ] <- corrAB
lM(object)$corrAB <- tmp
tmp <- physical(lM(object))$corrAA
tmp[snps, ] <- corrAA
lM(object)$corrAA <- tmp
tmp <- physical(lM(object))$corrBB
tmp[snps, ] <- corrBB
lM(object)$corrAB <- tmp
## lM(object)$tau2A[snps, ] <- tau2A
## lM(object)$tau2B[snps, ] <- tau2B
## lM(object)$sig2A[snps, ] <- sig2A
## lM(object)$sig2B[snps, ] <- sig2B
## lM(object)$nuA[snps, ] <- nuA
## lM(object)$nuB[snps, ] <- nuB
## lM(object)$phiA[snps, ] <- phiA
## lM(object)$phiB[snps, ] <- phiB
## lM(object)$phiPrimeA[snps, ] <- phiA2
## lM(object)$phiPrimeB[snps, ] <- phiB2
lapply(assayData(object), close)
lapply(lM(object), close)
TRUE
}
fit.lm4 <- function(idxBatch,
snpBatches,
XIndex,
object,
Ns,
normal,
snpflags,
batchSize,
SNRMin,
MIN.SAMPLES,
MIN.OBS,
DF.PRIOR,
GT.CONF.THR,
THR.NU.PHI,
MIN.NU,
MIN.PHI,
verbose, ...){
if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
open(object)
open(normal)
open(snpflags)
snps <- snpBatches[[idxBatch]]
batches <- split(seq(along=batch(object)), batch(object))
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
nuA <- phiA <- sig2A <- tau2A <- matrix(NA, length(snps), length(unique(batch(object))))
cA <- matrix(NA, length(snps), ncol(object))
ii <- isSnp(object) & chromosome(object) < 23 & !is.na(chromosome(object))
flags <- snpflags[ii, , drop=FALSE]
noflags <- rowSums(flags, na.rm=TRUE) == 0
lapply(lM(object), open)
nuIA <- physical(lM(object))$nuA[ii, ]
nuIB <- physical(lM(object))$nuB[ii, ]
phiIA <- physical(lM(object))$phiA[ii,]
phiIB <- physical(lM(object))$phiB[ii,]
i1 <- rowSums(nuIA < 20, na.rm=TRUE) == 0
i2 <- rowSums(nuIB < 20, na.rm=TRUE) == 0
i3 <- rowSums(phiIA < 20, na.rm=TRUE) == 0
i4 <- rowSums(phiIB < 20, na.rm=TRUE) == 0
snp.index <- which(i1 & i2 & i3 & i4 & noflags)
if(length(snp.index) == 0){
warning("No snps meet the following criteria: (1) nu and phi > 20 and (2) at least MIN.OBS in each genotype cluster. CN not estimated for nonpolymorphic loci on X")
return(TRUE)
}
if(length(snp.index) >= 5000){
snp.index <- sample(snp.index, 5000)
}
phiA.snp <- physical(lM(object))$phiA[snp.index, , drop=FALSE]
phiB.snp <- physical(lM(object))$phiB[snp.index, , drop=FALSE]
A.snp <- as.matrix(A(object)[snp.index, ])
B.snp <- as.matrix(B(object)[snp.index, ])
NORM.snp <- as.matrix(normal[snp.index, ])
NORM.np <- as.matrix(normal[snps, ])
gender <- object$gender
pseudoAR <- position(object)[snps] < 2709520 | (position(object)[snps] > 154584237 & position(object)[snps] < 154913754)
pseudoAR[is.na(pseudoAR)] <- FALSE
GG <- as.matrix(calls(object)[snp.index, ])
CP <- as.matrix(snpCallProbability(object)[snp.index, ])
AA.np <- as.matrix(A(object)[snps, ])
##if(missing(which.batches)) which.batches <- seq(along=batches)
##batches <- batches[which.batches]
for(k in batches){
##if(verbose) message("SNP batch ", ii, " of ", length(batches))
G <- GG[, k]
xx <- CP[, k]
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
G <- G*highConf*NORM.snp[, k]
##snps
AA <- A.snp[, k]
BB <- B.snp[, k]
##index <- GT.B <- GT.A <- vector("list", 3)
##names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
Ns <- applyByGenotype(matrix(1, nrow(G), ncol(G)), rowSums, G)
muA <- applyByGenotype(AA, rowMedians, G)
muB <- applyByGenotype(BB, rowMedians, G)
muA <- muA[, 1]
muB <- muB[, 3]
X <- cbind(1, log2(c(muA, muB)))
J <- match(unique(batch(object)[k]), unique(batch(object)))
Y <- log2(c(phiA.snp[, J], phiB.snp[, J]))
##--------------------------------------------------
##RS: need to fix
remove <- is.na(X[, 2]) | !is.finite(Y)
Y <- Y[!remove]
X <- X[!remove, ]
##--------------------------------------------------
betahat <- solve(crossprod(X), crossprod(X, Y))
##nonpolymorphic
A <- AA.np[, k]
gend <- gender[k]
A.M <- A[, gend==1]
mu1 <- rowMedians(A.M, na.rm=TRUE)
A.F <- A[, gend==2]
mu2 <- rowMedians(A.F, na.rm=TRUE)
mus <- log2(cbind(mu1, mu2))
X.men <- cbind(1, mus[, 1])
X.fem <- cbind(1, mus[, 2])
Yhat1 <- as.numeric(X.men %*% betahat)
Yhat2 <- as.numeric(X.fem %*% betahat)
phi1 <- 2^(Yhat1)
phi2 <- 2^(Yhat2)
nu1 <- 2^(mus[, 1]) - phi1
nu2 <- 2^(mus[, 2]) - 2*phi2
if(any(pseudoAR)){
nu1[pseudoAR] <- 2^(mus[pseudoAR, 1]) - 2*phi1[pseudoAR]
}
normal.f <- NORM.np[, k]
A.F <- A.F*normal.f[, gend==2]
A.F[A.F==0] <- NA
nuA[, J] <- nu2
phiA[, J] <- phi2
sig2A[, J] <- rowMAD(log2(A.F), na.rm=TRUE)^2
if(THR.NU.PHI){
nuA[nuA[, J] < MIN.NU, J] <- MIN.NU
phiA[phiA[, J] < MIN.PHI, J] <- MIN.PHI
}
CT1 <- 1/phi1*(A.M-nu1)
CT2 <- 1/phi2*(A.F-nu2)
tmp <- cA[, k]
tmp[, gend==1] <- CT1
tmp[, gend==2] <- CT2
cA[, k] <- tmp
rm(tmp, CT1, CT2, A.F, normal.f, G, AA, BB, Y, X, Ns)
gc()
}
cA[cA < 0.05] <- 0.05
cA[cA > 5] <- 5
cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
CA(object)[snps, ] <- cA
open(lM(object))
tmp <- physical(lM(object))$nuA
tmp[snps, ] <- nuA
lM(object)$nuA <- tmp
tmp <- physical(lM(object))$sig2A
tmp[snps, ] <- sig2A
lM(object)$sig2A <- tmp
tmp <- physical(lM(object))$phiA
tmp[snps, ] <- phiA
lM(object)$sig2A <- tmp
## lM(object)$sig2A[snps, ] <- sig2A
## lM(object)$nuA[snps, ] <- nuA
## lM(object)$phiA[snps, ] <- phiA
lapply(assayData(object), close)
lapply(lM(object), close)
TRUE
}
whichPlatform <- function(cdfName){
index <- grep("genomewidesnp", cdfName)
if(length(index) > 0){
platform <- "affymetrix"
} else{
index <- grep("human", cdfName)
platform <- "illumina"
}
return(platform)
}
cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
if(missing(cdfName)) stop("must specify cdfName")
pkgname <- getCrlmmAnnotationName(cdfName)
require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")
if (missing(sns)) sns <- basename(filenames)
if(verbose) message("Loading annotations for nonpolymorphic probes")
loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)
fid <- getVarInEnv("npProbesFid")
SKW <- initializeBigVector("crlmmSKW.NP-", length(filenames), "double")
##A <- initializeBigMatrix("crlmmNP-", length(fid), length(filenames), "integer")
sampleBatches <- splitIndicesByNode(seq(along=filenames))
if(verbose) message("Processing nonpolymorphic probes for ", length(filenames), " files.")
ocLapply(sampleBatches,
processCEL2,
row.names=row.names,
filenames=filenames,
A=A,
SKW=SKW,
seed=seed,
pkgname=pkgname,
cdfName=cdfName,
neededPkgs=c("crlmm", pkgname))
list(sns=sns, gns=row.names, SKW=SKW, cdfName=cdfName)
}
processCEL2 <- function(i, filenames, row.names, A, SKW, seed, cdfName, pkgname){
if(cdfName=="genomewidesnp6"){
loader("1m_reference_cn.rda", .crlmmPkgEnv, pkgname)
}
if(cdfName=="genomewidesnp5"){
loader("5.0_reference_cn.rda", .crlmmPkgEnv, pkgname)
}
reference <- getVarInEnv("reference")
loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)
fid <- getVarInEnv("npProbesFid")
fid <- fid[match(row.names, names(fid))]
np.index <- match(row.names, rownames(A))
gns <- names(fid)
set.seed(seed)
idx2 <- sample(length(fid), 10^5)
open(A)
open(SKW)
for (k in i){
y <- as.matrix(read.celfile(filenames[k], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])
x <- log2(y[idx2])
SKW[k] <- mean((x-mean(x))^3)/(sd(x)^3)
rm(x)
A[np.index, k] <- as.integer(normalize.quantiles.use.target(y, target=reference))
rm(y)
}
close(A)
close(SKW)
TRUE
}
# steps: quantile normalize hapmap: create 1m_reference_cn.rda object
cnrma <- function(filenames, cdfName, row.names, sns, seed=1, verbose=FALSE){
if(missing(cdfName)) stop("must specify cdfName")
pkgname <- getCrlmmAnnotationName(cdfName)
require(pkgname, character.only=TRUE) || stop("Package ", pkgname, " not available")
if (missing(sns)) sns <- basename(filenames)
loader("npProbesFid.rda", .crlmmPkgEnv, pkgname)
fid <- getVarInEnv("npProbesFid")
fid <- fid[match(row.names, names(fid))]
set.seed(seed)
idx2 <- sample(length(fid), 10^5) ##for skewness. no need to do everything
SKW <- vector("numeric", length(filenames))
NP <- matrix(NA, length(fid), length(filenames))
verbose <- TRUE
if(verbose){
message("Processing ", length(filenames), " files.")
if (getRversion() > '2.7.0') pb <- txtProgressBar(min=0, max=length(filenames), style=3)
}
if(cdfName=="genomewidesnp6"){
loader("1m_reference_cn.rda", .crlmmPkgEnv, pkgname)
}
if(cdfName=="genomewidesnp5"){
loader("5.0_reference_cn.rda", .crlmmPkgEnv, pkgname)
}
reference <- getVarInEnv("reference")
##if(!is.matrix(reference)) stop("target distribution for quantile normalization not available.")
for(i in seq(along=filenames)){
y <- as.matrix(read.celfile(filenames[i], intensity.means.only=TRUE)[["INTENSITY"]][["MEAN"]][fid])
x <- log2(y[idx2])
SKW[i] <- mean((x-mean(x))^3)/(sd(x)^3)
rm(x)
NP[, i] <- as.integer(normalize.quantiles.use.target(y, target=reference))
if (verbose)
if (getRversion() > '2.7.0') setTxtProgressBar(pb, i)
else cat(".")
}
dimnames(NP) <- list(names(fid), sns)
##dimnames(NP) <- list(map[, "man_fsetid"], sns)
##res3 <- list(NP=NP, SKW=SKW)
cat("\n")
return(NP)
}
getFlags <- function(object, PHI.THR){
batch <- unique(object$batch)
nuA <- getParam(object, "nuA", batch)
nuB <- getParam(object, "nuB", batch)
phiA <- getParam(object, "phiA", batch)
phiB <- getParam(object, "phiB", batch)
negativeNus <- nuA < 1 | nuB < 1
negativePhis <- phiA < PHI.THR | phiB < PHI.THR
negativeCoef <- negativeNus | negativePhis
notfinitePhi <- !is.finite(phiA) | !is.finite(phiB)
flags <- negativeCoef | notfinitePhi
return(flags)
}
instantiateObjects <- function(object, cnOptions){
Ns <- matrix(NA, nrow(object), 5)
colnames(Ns) <- c("A", "B", "AA", "AB", "BB")
vA <- vB <- muB <- muA <- Ns
normal <- matrix(TRUE, nrow(object), ncol(object))
dimnames(normal) <- list(featureNames(object), sampleNames(object))
tmp.objects <- list(vA=vA,
vB=vB,
muB=muB,
muA=muA,
Ns=Ns,
normal=normal)
return(tmp.objects)
}
thresholdCopynumber <- function(object){
ca <- CA(object)
cb <- CB(object)
ca[ca < 0.05] <- 0.05
ca[ca > 5] <- 5
cb[cb < 0.05] <- 0.05
cb[cb > 5] <- 5
CA(object) <- ca
CB(object) <- cb
return(object)
}
##linear model parameters
lm.parameters <- function(object, batch){##cnOptions){
fD <- fData(object)
##batch <- object$batch
uplate <- unique(batch)
parameterNames <- c(paste("tau2A", uplate, sep="_"),
paste("tau2B", uplate, sep="_"),
paste("sig2A", uplate, sep="_"),
paste("sig2B", uplate, sep="_"),
paste("nuA", uplate, sep="_"),
paste("nuA.se", uplate, sep="_"),
paste("nuB", uplate, sep="_"),
paste("nuB.se", uplate, sep="_"),
paste("phiA", uplate, sep="_"),
paste("phiA.se", uplate, sep="_"),
paste("phiB", uplate, sep="_"),
paste("phiB.se", uplate, sep="_"),
paste("phiAX", uplate, sep="_"),
paste("phiBX", uplate, sep="_"),
paste("corr", uplate, sep="_"),
paste("corrA.BB", uplate, sep="_"),
paste("corrB.AA", uplate, sep="_"))
pMatrix <- data.frame(matrix(numeric(0),
nrow(object),
length(parameterNames)),
row.names=featureNames(object))
colnames(pMatrix) <- parameterNames
fD2 <- cbind(fD, pMatrix)
new("AnnotatedDataFrame", data=fD2,
varMetadata=data.frame(labelDescription=colnames(fD2),
row.names=colnames(fD2)))
}
nonpolymorphic <- function(object, cnOptions, tmp.objects){
batch <- unique(object$batch)
CHR <- unique(chromosome(object))
goodSnps <- function(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR){
Ns <- tmp.objects[["Ns"]]
##Ns <- get("Ns", envir)
flags <- getFlags(object, PHI.THR)
fewAA <- Ns[, "AA"] < nAA.THR
fewBB <- Ns[, "BB"] < nBB.THR
flagsA <- flags | fewAA
flagsB <- flags | fewBB
flags <- list(A=flagsA, B=flagsB)
return(flags)
}
nAA.THR <- cnOptions$nHOM.THR
nBB.THR <- cnOptions$nHOM.THR
PHI.THR <- cnOptions$PHI.THR
snpflags <- goodSnps(object, PHI.THR, tmp.objects, nAA.THR, nBB.THR)
flagsA <- snpflags$A
flagsB <- snpflags$B
## if(all(flagsA) | all(flagsB)) stop("all snps are flagged")
nuA <- getParam(object, "nuA", batch)
nuB <- getParam(object, "nuB", batch)
phiA <- getParam(object, "phiA", batch)
phiB <- getParam(object, "phiB", batch)
sns <- sampleNames(object)
muA <- tmp.objects[["muA"]]
muB <- tmp.objects[["muB"]]
A <- A(object)
B <- B(object)
CA <- CA(object)
CB <- CB(object)
if(CHR == 23){
phiAX <- getParam(object, "phiAX", batch)
phiBX <- getParam(object, "phiBX", batch)
}
##---------------------------------------------------------------------------
## Train on unflagged SNPs
##---------------------------------------------------------------------------
##Might be best to train using the X chromosome, since for the
##X phi and nu have been adjusted for cross-hybridization
##plateInd <- plate == uplate[p]
##muA <- muA[!flagsA, p, c("A", "AA")]
##muB <- muB[!flagsB, p, c("B", "BB")]
muA <- muA[!flagsA, "AA"]
muB <- muB[!flagsB, "BB"]
X <- cbind(1, log2(c(muA, muB)))
Y <- log2(c(phiA[!flagsA], phiB[!flagsB]))
if(nrow(X) > 5000){
ix <- sample(1:nrow(X), 5000)
} else {
ix <- 1:nrow(X)
}
betahat <- solve(crossprod(X[ix, ]), crossprod(X[ix, ], Y[ix]))
normal <- tmp.objects[["normal"]][!isSnp(object), ]
if(CHR == 23){
##normalNP <- envir[["normalNP"]]
##normalNP <- normalNP[, plate==uplate[p]]
##nuT <- envir[["nuT"]]
##phiT <- envir[["phiT"]]
##cnvs <- envir[["cnvs"]]
##loader("cnProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
##cnProbes <- get("cnProbes", envir=.crlmmPkgEnv)
##cnProbes <- cnProbes[match(cnvs, rownames(cnProbes)), ]
##For build Hg18
##http://genome.ucsc.edu/cgi-bin/hgGateway
##pseudo-autosomal regions on X
##chrX:1-2,709,520 and chrX:154584237-154913754, respectively
##par:pseudo-autosomal regions
pseudoAR <- position(object) < 2709520 | (position(object) > 154584237 & position(object) < 154913754)
##pseudoAR <- cnProbes[, "position"] < 2709520 | (cnProbes[, "position"] > 154584237 & cnProbes[, "position"] < 154913754)
##in case some of the cnProbes are not annotated
pseudoAR[is.na(pseudoAR)] <- FALSE
pseudoAR <- pseudoAR[!isSnp(object)]
##gender <- envir[["gender"]]
gender <- object$gender
obj1 <- object[!isSnp(object), ]
A.male <- A(obj1[, gender==1])
mu1 <- rowMedians(A.male, na.rm=TRUE)
##mu1 <- rowMedians(NP[, gender=="male"], na.rm=TRUE)
##mu2 <- rowMedians(NP[, gender=="female"], na.rm=TRUE)
A.female <- A(obj1[, gender==2])
mu2 <- rowMedians(A.female, na.rm=TRUE)
mus <- log2(cbind(mu1, mu2))
X.men <- cbind(1, mus[, 1])
X.fem <- cbind(1, mus[, 2])
Yhat1 <- as.numeric(X.men %*% betahat)
Yhat2 <- as.numeric(X.fem %*% betahat)
phi1 <- 2^(Yhat1)
phi2 <- 2^(Yhat2)
nu1 <- 2^(mus[, 1]) - phi1
nu2 <- 2^(mus[, 2]) - 2*phi2
if(any(pseudoAR)){
nu1[pseudoAR] <- 2^(mus[pseudoAR, 1]) - 2*phi1[pseudoAR]
}
CT1 <- 1/phi1*(A.male-nu1)
CT2 <- 1/phi2*(A.female-nu2)
##CT2 <- 1/phi2*(NP[, gender=="female"]-nu2)
##CT1 <- matrix(as.integer(100*CT1), nrow(CT1), ncol(CT1))
##CT2 <- matrix(as.integer(100*CT2), nrow(CT2), ncol(CT2))
##CT <- envir[["CT"]]
CA <- CA(obj1)
CA[, gender==1] <- CT1
CA[, gender==2] <- CT2
CA(object)[!isSnp(object), ] <- CA
##CT[, plate==uplate[p] & gender=="male"] <- CT1
##CT[, plate==uplate[p] & gender=="female"] <- CT2
##envir[["CT"]] <- CT
##only using females to compute the variance
##normalNP[, gender=="male"] <- NA
normal[, gender==1] <- NA
sig2A <- getParam(object, "sig2A", batch)
normal.f <- normal[, object$gender==2]
sig2A[!isSnp(object)] <- rowMAD(log2(A.female*normal.f), na.rm=TRUE)^2
sig2A[!isSnp(object) & is.na(sig2A)] <- median(sig2A[!isSnp(object)], na.rm=TRUE)
##sig2T[, p] <- rowMAD(log2(NP*normalNP), na.rm=TRUE)^2
object <- pr(object, "sig2A", batch, sig2A)
nuA[!isSnp(object)] <- nu2
phiA[!isSnp(object)] <- phi2
THR.NU.PHI <- cnOptions$THR.NU.PHI
if(THR.NU.PHI){
verbose <- cnOptions$verbose
##Assign values to object
object <- pr(object, "nuA", batch, nuA)
object <- pr(object, "phiA", batch, phiA)
##if(verbose) message("Thresholding nu and phi")
object <- thresholdModelParams(object, cnOptions)
} else {
object <- pr(object, "nuA", batch, nuA)
object <- pr(object, "phiA", batch, phiA)
}
} else {
A <- A(object)[!isSnp(object), ]
mus <- rowMedians(A * normal, na.rm=TRUE)
crosshyb <- max(median(muA) - median(mus), 0)
X <- cbind(1, log2(mus+crosshyb))
logPhiT <- X %*% betahat
phiA[!isSnp(object)] <- 2^(logPhiT)
nuA[!isSnp(object)] <- mus-2*phiA[!isSnp(object)]
THR.NU.PHI <- cnOptions$THR.NU.PHI
if(THR.NU.PHI){
verbose <- cnOptions$verbose
##Assign values to object
object <- pr(object, "nuA", batch, nuA)
object <- pr(object, "phiA", batch, phiA)
##if(verbose) message("Thresholding nu and phi")
object <- thresholdModelParams(object, cnOptions)
##reassign values (now thresholded at MIN.NU and MIN.PHI
nuA <- getParam(object, "nuA", batch)
phiA <- getParam(object, "phiA", batch)
}
CA(object)[!isSnp(object), ] <- 1/phiA[!isSnp(object)]*(A - nuA[!isSnp(object)])
sig2A <- getParam(object, "sig2A", batch)
sig2A[!isSnp(object)] <- rowMAD(log2(A*normal), na.rm=TRUE)^2
object <- pr(object, "sig2A", batch, sig2A)
##added
object <- pr(object, "nuA", batch, nuA)
object <- pr(object, "phiA", batch, phiA)
}
return(object)
}
##sufficient statistics on the intensity scale
withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
normal <- tmp.objects[["normal"]]
## muA, muB: robust estimates of the within-genotype center (intensity scale)
muA <- tmp.objects[["muA"]]
muB <- tmp.objects[["muB"]]
## vA, vB: robust estimates of the within-genotype variance (intensity scale)
vA <- tmp.objects[["vA"]]
vB <- tmp.objects[["vB"]]
Ns <- tmp.objects[["Ns"]]
G <- snpCall(object)
GT.CONF.THR <- cnOptions$GT.CONF.THR
CHR <- unique(chromosome(object))
A <- A(object)
B <- B(object)
## highConf <- (1-exp(-confs(object)/1000)) > GT.CONF.THR
xx <- snpCallProbability(object)
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
##highConf <- confs(object) > GT.CONF.THR
##highConf <- highConf > GT.CONF.THR
if(CHR == 23){
gender <- object$gender
## gender <- envir[["gender"]]
IX <- matrix(gender, nrow(G), ncol(G), byrow=TRUE)
## IX <- IX == "female"
IX <- IX == 2 ##2=female, 1=male
} else IX <- matrix(TRUE, nrow(G), ncol(G))
index <- GT.B <- GT.A <- vector("list", 3)
names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
##--------------------------------------------------
##within-genotype sufficient statistics
##--------------------------------------------------
##GT.B <- GT.A <- list()
snpIndicator <- matrix(isSnp(object), nrow(object), ncol(object)) ##RS: added
for(j in 1:3){
GT <- G==j & highConf & IX & snpIndicator
GT <- GT * normal
Ns[, j+2] <- rowSums(GT, na.rm=TRUE)
GT[GT == FALSE] <- NA
GT.A[[j]] <- GT*A
GT.B[[j]] <- GT*B
index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added
muA[, j+2] <- rowMedians(GT.A[[j]], na.rm=TRUE)
muB[, j+2] <- rowMedians(GT.B[[j]], na.rm=TRUE)
vA[, j+2] <- rowMAD(GT.A[[j]], na.rm=TRUE)
vB[, j+2] <- rowMAD(GT.B[[j]], na.rm=TRUE)
##Shrink towards the typical variance
DF <- Ns[, j+2]-1
DF[DF < 1] <- 1
v0A <- median(vA[, j+2], na.rm=TRUE)
v0B <- median(vB[, j+2], na.rm=TRUE)
if(v0A == 0) v0A <- NA
if(v0B == 0) v0B <- NA
DF.PRIOR <- cnOptions$DF.PRIOR
vA[, j+2] <- (vA[, j+2]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
vA[is.na(vA[, j+2]), j+2] <- v0A
vB[, j+2] <- (vB[, j+2]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
vB[is.na(vB[, j+2]), j+2] <- v0B
}
if(CHR == 23){
k <- 1
for(j in c(1,3)){
GT <- G==j & highConf & !IX
Ns[, k] <- rowSums(GT)
GT[GT == FALSE] <- NA
muA[, k] <- rowMedians(GT*A, na.rm=TRUE)
muB[, k] <- rowMedians(GT*B, na.rm=TRUE)
vA[, k] <- rowMAD(GT*A, na.rm=TRUE)
vB[, k] <- rowMAD(GT*B, na.rm=TRUE)
DF <- Ns[, k]-1
DF[DF < 1] <- 1
v0A <- median(vA[, k], na.rm=TRUE)
v0B <- median(vB[, k], na.rm=TRUE)
vA[, k] <- (vA[, k]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
vA[is.na(vA[, k]), k] <- v0A
vB[, k] <- (vB[, k]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
vB[is.na(vB[, k]), k] <- v0B
k <- k+1
}
}
tmp.objects[["Ns"]] <- Ns
tmp.objects[["vA"]] <- vA
tmp.objects[["vB"]] <- vB
tmp.objects[["muA"]] <- muA
tmp.objects[["muB"]] <- muB
tmp.objects$index <- index
tmp.objects$GT.A <- GT.A
tmp.objects$GT.B <- GT.B
return(tmp.objects)
}
imputeCenter <- function(muA, muB, index.complete, index.missing){
index <- index.missing
mnA <- muA
mnB <- muB
for(j in 1:3){
if(length(index[[j]]) == 0) next()
X <- cbind(1, mnA[index.complete, -j, drop=FALSE], mnB[index.complete, -j, drop=FALSE])
Y <- cbind(mnA[index.complete, j], mnB[index.complete, j])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, mnA[index[[j]], -j, drop=FALSE], mnB[index[[j]], -j, drop=FALSE])
mus <- X %*% betahat
muA[index[[j]], j] <- mus[, 1]
muB[index[[j]], j] <- mus[, 2]
}
list(muA, muB)
}
imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
index1 <- which(Ns[, 1] == 0 & Ns[, 3] > MIN.OBS)
if(length(index1) > 0){
X <- cbind(1, muA[index.complete[[1]], 3], muB[index.complete[[1]], 3])
Y <- cbind(1, muA[index.complete[[1]], 1], muB[index.complete[[1]], 1])
betahat <- solve(crossprod(X), crossprod(X,Y))
##now with the incomplete SNPs
X <- cbind(1, muA[index1, 3], muB[index1, 3])
mus <- X %*% betahat
muA[index1, 1] <- mus[, 2]
muB[index1, 1] <- mus[, 3]
}
index1 <- which(Ns[, 3] == 0)
if(length(index1) > 0){
X <- cbind(1, muA[index.complete[[2]], 1], muB[index.complete[[2]], 1])
Y <- cbind(1, muA[index.complete[[2]], 3], muB[index.complete[[2]], 3])
betahat <- solve(crossprod(X), crossprod(X,Y))
##now with the incomplete SNPs
X <- cbind(1, muA[index1, 1], muB[index1, 1])
mus <- X %*% betahat
muA[index1, 3] <- mus[, 2]
muB[index1, 3] <- mus[, 3]
}
list(muA, muB)
}
oneBatch <- function(object, cnOptions, tmp.objects){
muA <- tmp.objects[["muA"]]
muB <- tmp.objects[["muB"]]
Ns <- tmp.objects[["Ns"]]
CHR <- unique(chromosome(object))
##---------------------------------------------------------------------------
## Impute sufficient statistics for unobserved genotypes (plate-specific)
##---------------------------------------------------------------------------
MIN.OBS <- cnOptions$MIN.OBS
index.AA <- which(Ns[, "AA"] >= MIN.OBS)
index.AB <- which(Ns[, "AB"] >= MIN.OBS)
index.BB <- which(Ns[, "BB"] >= MIN.OBS)
correct.orderA <- muA[, "AA"] > muA[, "BB"]
correct.orderB <- muB[, "BB"] > muB[, "AA"]
##For chr X, this will ignore the males
nobs <- rowSums(Ns[, 3:5] >= MIN.OBS, na.rm=TRUE) == 3
index.complete <- which(correct.orderA & correct.orderB & nobs) ##be selective here
size <- min(5000, length(index.complete))
if(size == 5000) index.complete <- sample(index.complete, 5000)
if(length(index.complete) < 200){
stop("fewer than 200 snps pass criteria for predicting the sufficient statistics")
}
index <- tmp.objects[["index"]]
index[[1]] <- which(Ns[, "AA"] == 0 & (Ns[, "AB"] >= MIN.OBS & Ns[, "BB"] >= MIN.OBS))
index[[2]] <- which(Ns[, "AB"] == 0 & (Ns[, "AA"] >= MIN.OBS & Ns[, "BB"] >= MIN.OBS))
index[[3]] <- which(Ns[, "BB"] == 0 & (Ns[, "AB"] >= MIN.OBS & Ns[, "AA"] >= MIN.OBS))
mnA <- muA[, 3:5]
mnB <- muB[, 3:5]
for(j in 1:3){
if(length(index[[j]]) == 0) next()
X <- cbind(1, mnA[index.complete, -j, drop=FALSE], mnB[index.complete, -j, drop=FALSE])
Y <- cbind(mnA[index.complete, j], mnB[index.complete, j])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, mnA[index[[j]], -j, drop=FALSE], mnB[index[[j]], -j, drop=FALSE])
mus <- X %*% betahat
muA[index[[j]], j+2] <- mus[, 1]
muB[index[[j]], j+2] <- mus[, 2]
}
nobsA <- Ns[, "A"] > MIN.OBS
nobsB <- Ns[, "B"] > MIN.OBS
notMissing <- !(is.na(muA[, "A"]) | is.na(muA[, "B"]) | is.na(muB[, "A"]) | is.na(muB[, "B"]))
complete <- list()
complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
size <- min(5000, length(complete[[1]]))
if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
if(CHR == 23){
index <- list()
index[[1]] <- which(Ns[, "A"] == 0)
index[[2]] <- which(Ns[, "B"] == 0)
cols <- 2:1
for(j in 1:2){
if(length(index[[j]]) == 0) next()
X <- cbind(1, muA[complete[[j]], cols[j]], muB[complete[[j]], cols[j]])
Y <- cbind(muA[complete[[j]], j], muB[complete[[j]], j])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, muA[index[[j]], cols[j]], muB[index[[j]], cols[j]])
mus <- X %*% betahat
muA[index[[j]], j] <- mus[, 1]
muB[index[[j]], j] <- mus[, 2]
}
}
##missing two genotypes
noAA <- Ns[, "AA"] < MIN.OBS
noAB <- Ns[, "AB"] < MIN.OBS
noBB <- Ns[, "BB"] < MIN.OBS
index[[1]] <- noAA & noAB
index[[2]] <- noBB & noAB
index[[3]] <- noAA & noBB
## snpflags <- envir[["snpflags"]]
snpflags <- index[[1]] | index[[2]] | index[[3]]
## snpflags[, p] <- index[[1]] | index[[2]] | index[[3]]
##---------------------------------------------------------------------------
## Two genotype clusters not observed -- would sequence help? (didn't seem to)
## 1. extract index of complete data
## 2. Regress mu_missing ~ sequence + mu_observed
## 3. solve for nu assuming the median is 2
##---------------------------------------------------------------------------
cols <- c(3, 1, 2)
for(j in 1:3){
if(sum(index[[j]]) == 0) next()
k <- cols[j]
X <- cbind(1, mnA[index.complete, k], mnB[index.complete, k])
Y <- cbind(mnA[index.complete, -k],
mnB[index.complete, -k])
betahat <- solve(crossprod(X), crossprod(X,Y))
X <- cbind(1, mnA[index[[j]], k], mnB[index[[j]], k])
mus <- X %*% betahat
muA[index[[j]], -c(1, 2, k+2)] <- mus[, 1:2]
muB[index[[j]], -c(1, 2, k+2)] <- mus[, 3:4]
}
negA <- rowSums(muA < 0) > 0
negB <- rowSums(muB < 0) > 0
snpflags <- snpflags| negA | negB | rowSums(is.na(muA[, 3:5]), na.rm=TRUE) > 0
tmp.objects$snpflags <- snpflags
tmp.objects[["muA"]] <- muA
tmp.objects[["muB"]] <- muB
tmp.objects[["snpflags"]] <- snpflags
return(tmp.objects)
}
##Estimate tau2, sigma2, and correlation (updates the object)
locationAndScale <- function(object, cnOptions, tmp.objects){
DF.PRIOR <- cnOptions$DF.PRIOR
Ns <- tmp.objects[["Ns"]]
index <- tmp.objects[["index"]]
index.AA <- index[[1]]
index.AB <- index[[2]]
index.BB <- index[[3]]
rm(index); gc()
GT.A <- tmp.objects[["GT.A"]]
GT.B <- tmp.objects[["GT.B"]]
AA.A <- GT.A[[1]]
AB.A <- GT.A[[2]]
BB.A <- GT.A[[3]]
AA.B <- GT.B[[1]]
AB.B <- GT.B[[2]]
BB.B <- GT.B[[3]]
x <- BB.A[index.BB, ]
##batch <- unique(object$batch)
batch <- unique(batch(object))
tau2A <- getParam(object, "tau2A", batch)
tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
DF <- Ns[, "BB"]-1
DF[DF < 1] <- 1
med <- median(tau2A, na.rm=TRUE)
tau2A <- (tau2A * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
tau2A[is.na(tau2A) & isSnp(object)] <- med
object <- pr(object, "tau2A", batch, tau2A)
sig2B <- getParam(object, "sig2B", batch)
x <- BB.B[index.BB, ]
sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
med <- median(sig2B, na.rm=TRUE)
sig2B <- (sig2B * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
sig2B[is.na(sig2B) & isSnp(object)] <- med
object <- pr(object, "sig2B", batch, sig2B)
tau2B <- getParam(object, "tau2B", batch)
x <- AA.B[index.AA, ]
tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
DF <- Ns[, "AA"]-1
DF[DF < 1] <- 1
med <- median(tau2B, na.rm=TRUE)
tau2B <- (tau2B * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
tau2B[is.na(tau2B) & isSnp(object)] <- med
object <- pr(object, "tau2B", batch, tau2B)
sig2A <- getParam(object, "sig2A", batch)
x <- AA.A[index.AA, ]
sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)
med <- median(sig2A, na.rm=TRUE)
sig2A <- (sig2A * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
sig2A[is.na(sig2A) & isSnp(object)] <- med
object <- pr(object, "sig2A", batch, sig2A)
if(length(index.AB) > 0){ ##all homozygous is possible
corr <- getParam(object, "corr", batch)
x <- AB.A[index.AB, ]
y <- AB.B[index.AB, ]
corr[index.AB] <- rowCors(x, y, na.rm=TRUE)
corr[corr < 0] <- 0
DF <- Ns[, "AB"]-1
DF[DF<1] <- 1
med <- median(corr, na.rm=TRUE)
corr <- (corr*DF + med * DF.PRIOR)/(DF.PRIOR + DF)
corr[is.na(corr) & isSnp(object)] <- med
object <- pr(object, "corr", batch, corr)
}
corrB.AA <- getParam(object, "corrB.AA", batch)
backgroundB <- AA.B[index.AA, ]
signalA <- AA.A[index.AA, ]
corrB.AA[index.AA] <- rowCors(backgroundB, signalA, na.rm=TRUE)
DF <- Ns[, "AA"]-1
DF[DF < 1] <- 1
med <- median(corrB.AA, na.rm=TRUE)
corrB.AA <- (corrB.AA*DF + med*DF.PRIOR)/(DF.PRIOR + DF)
corrB.AA[is.na(corrB.AA) & isSnp(object)] <- med
object <- pr(object, "corrB.AA", batch, corrB.AA)
corrA.BB <- getParam(object, "corrA.BB", batch)
backgroundA <- BB.A[index.BB, ]
signalB <- BB.B[index.BB, ]
corrA.BB[index.BB] <- rowCors(backgroundA, signalB, na.rm=TRUE)
DF <- Ns[, "BB"]-1
DF[DF < 1] <- 1
med <- median(corrA.BB, na.rm=TRUE)
corrA.BB <- (corrA.BB*DF + med*DF.PRIOR)/(DF.PRIOR + DF)
corrA.BB[is.na(corrA.BB) & isSnp(object)] <- med
object <- pr(object, "corrA.BB", batch, corrA.BB)
return(object)
}
coefs <- function(object, cnOptions, tmp.objects){
##batch <- unique(object$batch)
batch <- unique(batch(object))
CHR <- unique(chromosome(object))
muA <- tmp.objects[["muA"]]
muB <- tmp.objects[["muB"]]
vA <- tmp.objects[["vA"]]
vB <- tmp.objects[["vB"]]
Ns <- tmp.objects[["Ns"]]
if(CHR != 23){
IA <- muA[, 3:5]
IB <- muB[, 3:5]
vA <- vA[, 3:5]
vB <- vB[, 3:5]
Np <- Ns[, 3:5]
} else {
NOHET <- is.na(median(vA[, "AB"], na.rm=TRUE))
if(NOHET){
IA <- muA[, -4]
IB <- muB[, -4]
vA <- vA[, -4]
vB <- vB[, -4]
Np <- Ns[, -4]
} else{
IA <- muA
IB <- muB
vA <- vA
vB <- vB
Np <- Ns
}
}
Np[Np < 1] <- 1
vA2 <- vA^2/Np
vB2 <- vB^2/Np
wA <- sqrt(1/vA2)
wB <- sqrt(1/vB2)
YA <- IA*wA
YB <- IB*wB
##update lm.coefficients stored in object
object <- nuphiAllele(object, allele="A", Ystar=YA, W=wA, tmp.objects=tmp.objects, cnOptions=cnOptions)
object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)
##---------------------------------------------------------------------------
##Estimate crosshyb using X chromosome and sequence information
##---------------------------------------------------------------------------
##browser()
####data(sequences, package="genomewidesnp6Crlmm")
##snpflags <- envir[["snpflags"]]
##muA <- envir[["muA"]][, p, 3:5]
##muB <- envir[["muB"]][, p, 3:5]
##Y <- envir[["phiAx"]]
##load("sequences.rda")
##seqA <- sequences[, "A", ][, 1]
##seqA <- seqA[match(snps, names(seqA))]
##X <- cbind(1, sequenceDesignMatrix(seqA))
##X <- cbind(X, nuA[, p], phiA[, p], nuB[, p], phiB[, p])
##missing <- rowSums(is.na(X)) > 0
##betahat <- solve(crossprod(X[!missing, ]), crossprod(X[!missing, ], Y[!missing]))
return(object)
}
polymorphic <- function(object, cnOptions, tmp.objects){
##batch <- unique(object$batch)
batch <- unique(batch(object))
CHR <- unique(chromosome(object))
vA <- tmp.objects[["vA"]]
vB <- tmp.objects[["vB"]]
Ns <- tmp.objects[["Ns"]]
nuA <- getParam(object, "nuA", batch)
nuB <- getParam(object, "nuB", batch)
nuA.se <- getParam(object, "nuA.se", batch)
nuB.se <- getParam(object, "nuB.se", batch)
phiA <- getParam(object, "phiA", batch)
phiB <- getParam(object, "phiB", batch)
phiA.se <- getParam(object, "phiA.se", batch)
phiB.se <- getParam(object, "phiB.se", batch)
A <- A(object)
B <- B(object)
NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
##---------------------------------------------------------------------------
## Estimate CA, CB
##---------------------------------------------------------------------------
if(CHR == 23){
phiAX <- getParam(object, "phiAX", batch) ##nonspecific hybridization coef
phiBX <- getParam(object, "phiBX", batch) ##nonspecific hybridization coef
phistar <- phiBX/phiA
tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
copyB <- tmp/(1-phistar*phiAX/phiB)
copyA <- (A-nuA-phiAX*copyB)/phiA
CB(object) <- copyB ## multiplies by 100 and converts to integer
CA(object) <- copyA
} else{
CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
CB(object) <- matrix((1/phiB*(B-nuB)), nrow(B), ncol(B))
}
return(object)
}
posteriorProbability.snps <- function(object, cnOptions, tmp.objects=list()){
I <- isSnp(object)
gender <- object$gender
CHR <- unique(chromosome(object))
##batch <- unique(object$batch)
batch <- unique(batch(object))
if(CHR == 23){
phiAX <- getParam(object, "phiAX", batch)[I]
phiBX <- getParam(object, "phiBX", batch)[I]
}
A <- A(object)[I, ]
B <- B(object)[I, ]
sig2A <- getParam(object, "sig2A", batch)[I]
sig2B <- getParam(object, "sig2B", batch)[I]
tau2A <- getParam(object, "tau2A", batch)[I]
tau2B <- getParam(object, "tau2B", batch)[I]
corrA.BB <- getParam(object, "corrA.BB", batch)[I]
corrB.AA <- getParam(object, "corrB.AA", batch)[I]
corr <- getParam(object, "corr", batch)[I]
nuA <- getParam(object, "nuA", batch)[I]
nuB <- getParam(object, "nuB", batch)[I]
phiA <- getParam(object, "phiA", batch)[I]
phiB <- getParam(object, "phiB", batch)[I]
normal <- tmp.objects[["normal"]][I, ]
prior.prob <- cnOptions$prior.prob
emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'
lA <- log2(A)
lB <- log2(B)
X <- cbind(lA, lB)
counter <- 1##state counter
for(CT in 0:3){
for(CA in 0:CT){
cat(".")
CB <- CT-CA
A.scale <- sqrt(tau2A*(CA==0) + sig2A*(CA > 0))
B.scale <- sqrt(tau2B*(CA==0) + sig2B*(CA > 0))
if(CA == 0 & CB == 0) rho <- 0
if(CA == 0 & CB > 0) rho <- corrA.BB
if(CA > 0 & CB == 0) rho <- corrB.AA
if(CA > 0 & CB > 0) rho <- corr
if(CHR == 23){
##means <- cbind(suppressWarnings(log2(nuA[, p]+CA*phiA[, p] + CB*phiAx[, p])), suppressWarnings(log2(nuB[, p]+CB*phiB[, p] + CA*phiBx[, p])))
meanA <- suppressWarnings(log2(nuA+CA*phiA + CB*phiAX))
meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))
} else{
##means <- cbind(suppressWarnings(log2(nuA+CA*phiA)), suppressWarnings(log2(nuB+CB*phiB)))
meanA <- suppressWarnings(log2(nuA+CA*phiA))
meanB <- suppressWarnings(log2(nuB+CB*phiB))
covs <- rho*A.scale*B.scale
A.scale2 <- A.scale^2
B.scale2 <- B.scale^2
}
Q.x.y <- 1/(1-rho^2)*(((lA - meanA)/A.scale)^2 + ((lB - meanB)/B.scale)^2 - 2*rho*((lA - meanA)*(lB - meanB))/(A.scale*B.scale))
f.x.y <- 1/(2*pi*A.scale*B.scale*sqrt(1-rho^2))*exp(-0.5*Q.x.y)
emit[, , counter] <- f.x.y
counter <- counter+1
}
}
priorProb <- cnOptions$prior.prob
homDel <- priorProb[1]*emit[, , 1]
hemDel <- priorProb[2]*emit[, , c(2, 3)] # + priorProb[3]*emit[, c(4, 5, 6)] + priorProb[4]*emit[, c(7:10)]
norm <- priorProb[3]*emit[, , 4:6]
amp <- priorProb[4]*emit[, , 7:10]
##sum over the different combinations within each copy number state
hemDel <- hemDel[, , 1] + hemDel[, , 2]
norm <- norm[, , 1] + norm[, , 2] + norm[, , 3]
amp <- amp[, , 1] + amp[, , 2] + amp[ , , 3] + amp[, , 4]
total <- homDel + hemDel + norm + amp
homDel <- homDel/total
hemDel <- hemDel/total
norm <- norm/total
amp <- amp/total
tmp.objects$posteriorProb <- list(hemDel=hemDel, norm=norm, amp=amp)
##envir[["posteriorProb"]] <- list(hemDel=hemDel, norm=norm, amp=amp)
posteriorProb <- array(NA, dim=c(nrow(A), ncol(A), 4))
posteriorProb[, , 1] <- homDel
posteriorProb[, , 2] <- hemDel
posteriorProb[, , 3] <- norm
posteriorProb[, , 4] <- amp
return(list(tmp.objects, posteriorProb))
}
biasAdj <- function(object, cnOptions, tmp.objects){
gender <- object$gender
CHR <- unique(chromosome(object))
I <- isSnp(object)
A <- A(object)[I, ]
normal <- tmp.objects[["normal"]][I, ]
results <- posteriorProbability.snps(object=object, cnOptions=cnOptions, tmp.objects=tmp.objects)
posteriorProb <- results[[2]]
tmp.objects <- results[[1]]
mostLikelyState <- apply(posteriorProb, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
if(CHR == 23){
##so state index 3 is the most likely state for men and women
mostLikelyState[, gender==1] <- mostLikelyState[, gender==1] + 1
}
proportionSamplesAltered <- rowMeans(mostLikelyState != 3)
ii <- proportionSamplesAltered < 0.8 & proportionSamplesAltered > 0.01
## only exclude observations from one tail, depending on
## whether more are up or down
moreup <- rowSums(mostLikelyState > 3) > rowSums(mostLikelyState < 3) ##3 is normal
## if equal, which points have a high posterior probability of altered copy number. drop those.
NORM <- matrix(FALSE, nrow(A), ncol(A))
##NORM[proportionSamplesAltered > 0.8, ] <- FALSE
ratioUp <- posteriorProb[, , 4]/posteriorProb[, , 3]
NORM[ii & moreup, ] <- ratioUp[moreup & ii] < 1 ##normal more likely
ratioDown <- posteriorProb[, , 2]/posteriorProb[, , 3]
NORM[ii & !moreup, ] <- ratioDown[!moreup & ii] < 1 ##normal more likely
normal <- NORM*normal
tmp <- tmp.objects[["normal"]]
tmp[I, ] <- normal
tmp.objects[["normal"]] <- tmp
return(tmp.objects)
}
bias1 <- function(idxBatch,
snpBatches,
index,
object,
normal,
emit,
prior.prob,
MIN.SAMPLES,
verbose){
}
bias2 <- function(idxBatch,
snpBatches,
index,
object,
normal,
prior.prob,
MIN.SAMPLES,
verbose){
open(object)
open(normal)
nps <- snpBatches[[idxBatch]]
nuA <- lM(object)$nuA[nps, , drop=FALSE]
phiA <- lM(object)$phiA[nps, , drop=FALSE]
sig2A <- lM(object)$sig2A[nps, , drop=FALSE]
AA <- as.matrix(A(object)[nps, ])
batches <- split(seq(along=batch(object)), batch(object))
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
cn.lik <- matrix(NA, length(nps)*ncol(object), 4)
argmax.cn <- emit[nps, ]
norm <- matrix(1L, length(nps), ncol(object))
for(k in batches){
J <- match(unique(batch(object)[k]), unique(batch(object)))
lT <- log2(AA[, k])
counter <- 1 ##state counter
for(CT in c(0, 1.5, 2, 2.5)){
##sds <- sqrt(sig2A[, J]*(CT==0) + sig2A[ , J]*(CT > 0))
sds <- sqrt(sig2A[, J])
means <- suppressWarnings(log2(nuA[, J]+CT*phiA[, J]))
lik <- log(dnorm(lT, mean=means, sd=sds))
##emit[[counter]][nps, ] <- tmp
cn.lik[, counter] <- as.numeric(lik)
counter <- counter+1
}
outlier <- matrix(rowSums(cn.lik < -10) == 4, length(nps), ncol(object))
argmax.cn.lik <- apply(cn.lik, 1, function(x) order(x, decreasing=TRUE)[1])
argmax.cn <- matrix(argmax.cn.lik, length(nps), length(k))
isUp <- argmax.cn > 3
prUp <- rowMeans(isUp)
isDn <- argmax.cn < 3
prDn <- rowMeans(isDn)
index <- which(prUp > 0.05 & prUp > prDn)
##if proportion up greater than 5%, trim the high cn est.
norm[index, k] <- argmax.cn[index, ] > 3
index <- which(prDn > 0.05 & prDn > prUp)
norm[index, k] <- argmax.cn[index, ] < 3
norm[index, k] <- norm[index, k]*!outlier
}
normal[nps, ] <- norm
TRUE
}
biasAdjust <- function(object, prior.prob=rep(1/4, 4), MIN.SAMPLES=10, verbose=TRUE){
load(file.path(ldPath(), "normal.rda"))
autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]
## emit <- initializeBigMatrix("emit",
## nrow(object),
## ncol(object),
## vmode="double")
if(verbose) message("Bias adjustment for nonpolymorphic loci on chromosomes 1-22.")
snpBatches <- splitIndicesByLength(autosomeIndex.nps, ocProbesets())
ocLapply(seq(along=snpBatches),
bias2,
index=autosomeIndex.nps,
snpBatches=snpBatches,
object=object,
normal=normal,
prior.prob=prior.prob,
MIN.SAMPLES=MIN.SAMPLES,
verbose=verbose)
if(verbose) message("Bias adjustment for polymorphic loci on chromosomes 1-22.")
autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]
snpBatches <- splitIndicesByLength(autosomeIndex.snps, ocProbesets())
ocLapply(seq(along=snpBatches),
bias1,
index=autosomeIndex.snps,
snpBatches=snpBatches,
object=object,
normal=normal,
prior.prob=prior.prob,
emit=emit,
MIN.SAMPLES=MIN.SAMPLES,
verbose=verbose)
}
##biasAdjNP <- function(plateIndex, envir, priorProb){
biasAdjNP <- function(object, cnOptions, tmp.objects){
##batch <- unique(object$batch)
batch <- unique(batch(object))
normalNP <- tmp.objects[["normal"]][!isSnp(object), ]
CHR <- unique(chromosome(object))
A <- A(object)[!isSnp(object), ]
sig2A <- getParam(object, "sig2A", batch)
gender <- object$gender
##Assume that on the log-scale, that the background variance is the same...
tau2A <- sig2A
nuA <- getParam(object, "nuA", batch)
phiA <- getParam(object, "phiA", batch)
prior.prob <- cnOptions$prior.prob
emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'
lT <- log2(A)
I <- isSnp(object)
counter <- 1 ##state counter
## for(CT in 0:3){
## sds <- sqrt(tau2A[I]*(CT==0) + sig2A[I]*(CT > 0))
## means <- suppressWarnings(log2(nuA[I]+CT*phiA[I]))
## tmp <- dnorm(lT, mean=means, sd=sds)
## emit[, , counter] <- tmp
## counter <- counter+1
## }
## mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
counter <- 1
for(CT in c(0,1,2,2.5)){
sds <- sqrt(tau2A[I]*(CT==0) + sig2A[I]*(CT > 0))
means <- suppressWarnings(log2(nuA[I]+CT*phiA[I]))
tmp <- dnorm(lT, mean=means, sd=sds)
emit[, , counter] <- tmp
counter <- counter+1
}
mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
if(CHR == 23){
## the state index for male on chromosome 23 is 2
## add 1 so that the state index is 3 for 'normal' state
mostLikelyState[, gender=="male"] <- mostLikelyState[, gender==1] + 1
}
tmp3 <- mostLikelyState != 3
##Those near 1 have NaNs for nu and phi. this occurs by NaNs in the muA[,, "A"] or muA[, , "B"] for X chromosome
proportionSamplesAltered <- rowMeans(tmp3)##prop normal
ii <- proportionSamplesAltered < 0.75
moreup <- rowSums(mostLikelyState > 3) > rowSums(mostLikelyState < 3)
notUp <- mostLikelyState[ii & moreup, ] <= 3
notDown <- mostLikelyState[ii & !moreup, ] >= 3
NORM <- matrix(TRUE, nrow(A), ncol(A))
NORM[ii & moreup, ] <- notUp
NORM[ii & !moreup, ] <- notDown
normalNP <- normalNP*NORM
##flagAltered <- which(proportionSamplesAltered > 0.5)
##envir[["flagAlteredNP"]] <- flagAltered
normal <- tmp.objects[["normal"]]
normal[!isSnp(object), ] <- normalNP
tmp.objects[["normal"]] <- normal
return(tmp.objects)
}
getParams <- function(object, batch){
##batch <- unique(object$batch)
batch <- unique(batch(object))
if(length(batch) > 1) stop("batch variable not unique")
nuA <- as.numeric(fData(object)[, match(paste("nuA", batch, sep="_"), fvarLabels(object))])
nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])
phiA <- as.numeric(fData(object)[, match(paste("phiA", batch, sep="_"), fvarLabels(object))])
phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])
tau2A <- as.numeric(fData(object)[, match(paste("tau2A", batch, sep="_"), fvarLabels(object))])
tau2B <- as.numeric(fData(object)[, match(paste("tau2B", batch, sep="_"), fvarLabels(object))])
sig2A <- as.numeric(fData(object)[, match(paste("sig2A", batch, sep="_"), fvarLabels(object))])
sig2B <- as.numeric(fData(object)[, match(paste("sig2B", batch, sep="_"), fvarLabels(object))])
corrA.BB <- as.numeric(fData(object)[, match(paste("corrA.BB", batch, sep="_"), fvarLabels(object))])
corrB.AA <- as.numeric(fData(object)[, match(paste("corrB.AA", batch, sep="_"), fvarLabels(object))])
corr <- as.numeric(fData(object)[, match(paste("corr", batch, sep="_"), fvarLabels(object))])
params <- list(nuA=nuA,
nuB=nuB,
phiA=phiA,
phiB=phiB,
tau2A=tau2A,
tau2B=tau2B,
sig2A=sig2A,
sig2B=sig2B,
corrA.BB=corrA.BB,
corrB.AA=corrB.AA,
corr=corr)
return(params)
}
## Constrain nu and phi to positive values
thresholdModelParams <- function(object, cnOptions){
MIN.NU <- cnOptions$MIN.NU
MIN.PHI <- cnOptions$MIN.PHI
batch <- unique(object$batch)
nuA <- getParam(object, "nuA", batch)
nuA[nuA < MIN.NU] <- MIN.NU
object <- pr(object, "nuA", batch, nuA)
nuB <- getParam(object, "nuB", batch)
if(!all(is.na(nuB))){
nuB[nuB < MIN.NU] <- MIN.NU
object <- pr(object, "nuB", batch, nuB)
}
phiA <- getParam(object, "phiA", batch)
phiA[phiA < MIN.PHI] <- MIN.PHI
object <- pr(object, "phiA", batch, phiA)
phiB <- getParam(object, "phiB", batch)
if(!all(is.na(phiB))){
phiB[phiB < MIN.PHI] <- MIN.PHI
object <- pr(object, "phiB", batch, phiB)
}
phiAX <- as.numeric(getParam(object, "phiAX", batch))
if(!all(is.na(phiAX))){
phiAX[phiAX < MIN.PHI] <- MIN.PHI
object <- pr(object, "phiAX", batch, phiAX)
}
phiBX <- as.numeric(getParam(object, "phiBX", batch))
if(!all(is.na(phiBX))){
phiBX[phiBX < MIN.PHI] <- MIN.PHI
object <- pr(object, "phiBX", batch, phiBX)
}
return(object)
}
##computeCopynumber.SnpSuperSet <- function(object, cnOptions){
#### use.ff <- cnOptions[["use.ff"]]
#### if(!use.ff){
#### object <- as(object, "CrlmmSet")
#### } else object <- as(object, "CrlmmSetFF")
## bias.adj <- cnOptions[["bias.adj"]]
## ##must be FALSE to initialize parameters
## cnOptions[["bias.adj"]] <- FALSE
## ## Add linear model parameters to the CrlmmSet object
## featureData(object) <- lm.parameters(object, cnOptions)
## if(!isValidCdfName(annotation(object))) stop(annotation(object), " not supported.")
## object <- as(object, "CNSet")
## object <- computeCopynumber.CNSet(object, cnOptions)
## if(bias.adj==TRUE){## run a second time
## object <- computeCopynumber.CNSet(object, cnOptions)
## }
## return(object)
##}
computeCopynumber.CNSet <- function(object, cnOptions){
##PLATE <- unique(object$batch)
PLATE <- unique(batch(object))
verbose <- cnOptions$verbose
tmp.objects <- instantiateObjects(object, cnOptions)
bias.adj <- cnOptions$bias.adj
if(bias.adj & ncol(object) <= 15){
warning(paste("bias.adj is TRUE, but too few samples to perform this step"))
cnOptions$bias.adj <- bias.adj <- FALSE
}
if(bias.adj){
if(verbose) message("Dropping samples with low posterior prob. of normal copy number (samples dropped is locus-specific)")
tmp.objects <- biasAdjNP(object, cnOptions, tmp.objects)
tmp.objects <- biasAdj(object, cnOptions, tmp.objects)
if(verbose) message("Recomputing location and scale parameters")
}
##update tmp.objects
tmp.objects <- withinGenotypeMoments(object,
cnOptions=cnOptions,
tmp.objects=tmp.objects)
object <- locationAndScale(object, cnOptions, tmp.objects)
tmp.objects <- oneBatch(object,
cnOptions=cnOptions,
tmp.objects=tmp.objects)
##coefs calls nuphiAllele.
object <- coefs(object, cnOptions, tmp.objects)
##nuA=getParam(object, "nuA", PLATE)
THR.NU.PHI <- cnOptions$THR.NU.PHI
if(THR.NU.PHI){
verbose <- cnOptions$verbose
##if(verbose) message("Thresholding nu and phi")
object <- thresholdModelParams(object, cnOptions)
}
##if(verbose) message("\nAllele specific copy number")
object <- polymorphic(object, cnOptions, tmp.objects)
if(any(!isSnp(object))){ ## there are nonpolymorphic probes
##if(verbose) message("\nCopy number for nonpolymorphic probes...")
object <- nonpolymorphic(object, cnOptions, tmp.objects)
}
##---------------------------------------------------------------------------
##Note: the replacement method multiples by 100
## CA(object)[, batch==PLATE] <- CA(object)
## CB(object)[, batch==PLATE] <- CB(object)
##---------------------------------------------------------------------------
##update-the plate-specific parameters for copy number
object <- pr(object, "nuA", PLATE, getParam(object, "nuA", PLATE))
object <- pr(object, "nuA.se", PLATE, getParam(object, "nuA.se", PLATE))
object <- pr(object, "nuB", PLATE, getParam(object, "nuB", PLATE))
object <- pr(object, "nuB.se", PLATE, getParam(object, "nuB.se", PLATE))
object <- pr(object, "phiA", PLATE, getParam(object, "phiA", PLATE))
object <- pr(object, "phiA.se", PLATE, getParam(object, "phiA.se", PLATE))
object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))
object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))
object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
object <- pr(object, "corrA.BB", PLATE, getParam(object, "corrA.BB", PLATE))
object <- pr(object, "corrB.AA", PLATE, getParam(object, "corrB.AA", PLATE))
##object <- object[order(chromosome(object), position(object)), ]
if(cnOptions[["thresholdCopynumber"]]){
object <- thresholdCopynumber(object)
}
return(object)
}
constructFeatureData <- function(gns, cdfName){
pkgname <- paste(cdfName, "Crlmm", sep="")
path <- system.file("extdata", package=pkgname)
load(file.path(path, "cnProbes.rda"))
load(file.path(path, "snpProbes.rda"))
cnProbes$chr <- chromosome2integer(cnProbes$chr)
cnProbes <- as.matrix(cnProbes)
snpProbes$chr <- chromosome2integer(snpProbes$chr)
snpProbes <- as.matrix(snpProbes)
mapping <- rbind(snpProbes, cnProbes, deparse.level=0)
mapping <- mapping[match(gns, rownames(mapping)), ]
isSnp <- 1L-as.integer(gns %in% rownames(cnProbes))
mapping <- cbind(mapping, isSnp, deparse.level=0)
stopifnot(identical(rownames(mapping), gns))
colnames(mapping) <- c("chromosome", "position", "isSnp")
new("AnnotatedDataFrame",
data=data.frame(mapping),
varMetadata=data.frame(labelDescription=colnames(mapping)))
}
constructAssayData <- function(np, snp, object, storage.mode="environment", order.index){
stopifnot(identical(snp$gns, featureNames(object)))
gns <- c(featureNames(object), np$gns)
sns <- np$sns
np <- np[1:2]
snp <- snp[1:2]
stripnames <- function(x) {
dimnames(x) <- NULL
x
}
np <- lapply(np, stripnames)
snp <- lapply(snp, stripnames)
A <- rbind(snp[[1]], np[[1]], deparse.level=0)[order.index, ]
B <- rbind(snp[[2]], np[[2]], deparse.level=0)[order.index, ]
gt <- stripnames(calls(object))
emptyMatrix <- matrix(integer(), nrow(np[[1]]), ncol(A))
gt <- rbind(gt, emptyMatrix, deparse.level=0)[order.index,]
pr <- stripnames(snpCallProbability(object))
pr <- rbind(pr, emptyMatrix, deparse.level=0)[order.index, ]
emptyMatrix <- matrix(integer(), nrow(A), ncol(A))
aD <- assayDataNew(storage.mode,
alleleA=A,
alleleB=B,
call=gt,
callProbability=pr,
CA=emptyMatrix,
CB=emptyMatrix)
}
