... |
... |
@@ -28,7 +28,7 @@ getFeatureData.Affy <- function(cdfName, copynumber=FALSE){
|
28 |
28 |
load(file.path(path, "cnProbes.rda"))
|
29 |
29 |
cnProbes <- get("cnProbes")
|
30 |
30 |
snpIndex <- seq(along=gns)
|
31 |
|
- npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex)
|
|
31 |
+ npIndex <- seq(along=rownames(cnProbes)) + max(snpIndex)
|
32 |
32 |
featurenames <- c(gns, rownames(cnProbes))
|
33 |
33 |
} else featurenames <- gns
|
34 |
34 |
fvarlabels=c("chromosome", "position", "isSnp")
|
... |
... |
@@ -84,7 +84,7 @@ construct <- function(filenames, cdfName, copynumber=FALSE,
|
84 |
84 |
colnames(CA) <- sns
|
85 |
85 |
CB=initializeBigMatrix(name="CB", nr, nc)
|
86 |
86 |
colnames(CB) <- sns
|
87 |
|
- callSet <- new("CNSetLM",
|
|
87 |
+ callSet <- new("CNSetLM",
|
88 |
88 |
alleleA=alleleA,
|
89 |
89 |
alleleB=alleleB,
|
90 |
90 |
call=call,
|
... |
... |
@@ -157,7 +157,7 @@ genotype <- function(filenames,
|
157 |
157 |
np.index <- which(isSnp(callSet) == 0)
|
158 |
158 |
cnrmaRes <- cnrma(filenames=filenames[j],
|
159 |
159 |
cdfName=cdfName,
|
160 |
|
- row.names=featureNames(callSet)[np.index],
|
|
160 |
+ row.names=featureNames(callSet)[np.index],
|
161 |
161 |
sns=sns[j],
|
162 |
162 |
seed=seed,
|
163 |
163 |
verbose=verbose)
|
... |
... |
@@ -459,7 +459,7 @@ getPhenoData <- function(sampleSheet=NULL, arrayNames=NULL,
|
459 |
459 |
barcode = sampleSheet[,arrayInfoColNames$barcode]
|
460 |
460 |
arrayNames=barcode
|
461 |
461 |
}
|
462 |
|
- if(!is.null(arrayInfoColNames$position) && (arrayInfoColNames$position %in% colnames(sampleSheet))) {
|
|
462 |
+ if(!is.null(arrayInfoColNames$position) && (arrayInfoColNames$position %in% colnames(sampleSheet))) {
|
463 |
463 |
position = sampleSheet[,arrayInfoColNames$position]
|
464 |
464 |
if(is.null(arrayNames))
|
465 |
465 |
arrayNames=position
|
... |
... |
@@ -473,7 +473,7 @@ getPhenoData <- function(sampleSheet=NULL, arrayNames=NULL,
|
473 |
473 |
}
|
474 |
474 |
}
|
475 |
475 |
pd = new("AnnotatedDataFrame", data = sampleSheet)
|
476 |
|
- sampleNames(pd) <- basename(arrayNames)
|
|
476 |
+ sampleNames(pd) <- basename(arrayNames)
|
477 |
477 |
}
|
478 |
478 |
if(is.null(arrayNames)) {
|
479 |
479 |
arrayNames = gsub(paste(sep, fileExt$green, sep=""), "", dir(pattern=fileExt$green, path=path))
|
... |
... |
@@ -528,7 +528,7 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,
|
528 |
528 |
fileExt=list(green="Grn.idat", red="Red.idat"),
|
529 |
529 |
stripNorm=TRUE,
|
530 |
530 |
useTarget=TRUE,
|
531 |
|
- row.names=TRUE,
|
|
531 |
+ row.names=TRUE,
|
532 |
532 |
col.names=TRUE,
|
533 |
533 |
probs=c(1/3, 1/3, 1/3), DF=6, SNRMin=5, gender=NULL,
|
534 |
534 |
seed=1, save.ab=FALSE, snpFile, cnFile,
|
... |
... |
@@ -545,7 +545,7 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,
|
545 |
545 |
} else {
|
546 |
546 |
if(length(batch) != length(sns))
|
547 |
547 |
stop("batch variable must be the same length as the filenames")
|
548 |
|
- }
|
|
548 |
+ }
|
549 |
549 |
batches <- splitIndicesByLength(seq(along=arrayNames), ocSamples())
|
550 |
550 |
k <- 1
|
551 |
551 |
for(j in batches){
|
... |
... |
@@ -618,7 +618,7 @@ crlmmIlluminaRS <- function(sampleSheet=NULL,
|
618 |
618 |
protocolData(callSet)$batch[j] <- as.numeric(as.factor(protocolData$ScanDate))
|
619 |
619 |
}
|
620 |
620 |
## MR: we need to define a snp.index vs np.index
|
621 |
|
- snp.index <- match(res$gns, featureNames(callSet))
|
|
621 |
+ snp.index <- match(res$gns, featureNames(callSet))
|
622 |
622 |
A(callSet)[snp.index, j] <- res[["A"]]
|
623 |
623 |
B(callSet)[snp.index, j] <- res[["B"]]
|
624 |
624 |
pData(callSet)$SKW[j] <- res$SKW
|
... |
... |
@@ -733,7 +733,7 @@ generateIXTX <- function(x, nrow=3) {
|
733 |
733 |
}
|
734 |
734 |
|
735 |
735 |
nuphiAllele2 <- function(allele, Ystar, W, Ns){
|
736 |
|
- complete <- rowSums(is.na(W)) == 0
|
|
736 |
+ complete <- rowSums(is.na(W)) == 0
|
737 |
737 |
if(any(!is.finite(W))){## | any(!is.finite(V))){
|
738 |
738 |
i <- which(rowSums(!is.finite(W)) > 0)
|
739 |
739 |
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
|
... |
... |
@@ -741,7 +741,7 @@ nuphiAllele2 <- function(allele, Ystar, W, Ns){
|
741 |
741 |
NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
|
742 |
742 |
if(missing(allele)) stop("must specify allele")
|
743 |
743 |
if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
|
744 |
|
- if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
|
744 |
+ if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
745 |
745 |
##How to quickly generate Xstar, Xstar = diag(W) %*% X
|
746 |
746 |
Xstar <- apply(W, 1, generateX, X)
|
747 |
747 |
IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
|
... |
... |
@@ -758,7 +758,7 @@ nuphiAllele2 <- function(allele, Ystar, W, Ns){
|
758 |
758 |
}
|
759 |
759 |
|
760 |
760 |
nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
|
761 |
|
- complete <- rowSums(is.na(W)) == 0
|
|
761 |
+ complete <- rowSums(is.na(W)) == 0
|
762 |
762 |
if(any(!is.finite(W))){## | any(!is.finite(V))){
|
763 |
763 |
i <- which(rowSums(!is.finite(W)) > 0)
|
764 |
764 |
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
|
... |
... |
@@ -766,14 +766,14 @@ nuphiAlleleX <- function(allele, Ystar, W, Ns, chrX=FALSE){
|
766 |
766 |
##NOHET <- mean(Ns[, 2], na.rm=TRUE) < 0.05
|
767 |
767 |
if(missing(allele)) stop("must specify allele")
|
768 |
768 |
if(allele == "A") X <- cbind(1, c(1, 0, 2, 1, 0), c(0, 1, 0, 1, 2))
|
769 |
|
- if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
|
770 |
|
- ##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
|
769 |
+ if(allele == "B") X <- cbind(1, c(0, 1, 0, 1, 2), c(1, 0, 2, 1, 0))
|
|
770 |
+ ##if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
771 |
771 |
##How to quickly generate Xstar, Xstar = diag(W) %*% X
|
772 |
772 |
Xstar <- apply(W, 1, generateX, X)
|
773 |
773 |
IXTX <- apply(Xstar, 2, generateIXTX, nrow=nrow(X))
|
774 |
774 |
|
775 |
775 |
betahat <- matrix(NA, 3, nrow(Ystar))
|
776 |
|
- ses <- matrix(NA, 3, nrow(Ystar))
|
|
776 |
+ ses <- matrix(NA, 3, nrow(Ystar))
|
777 |
777 |
##betahat <- matrix(NA, 2, nrow(Ystar))
|
778 |
778 |
##ses <- matrix(NA, 2, nrow(Ystar))
|
779 |
779 |
for(i in 1:nrow(Ystar)){
|
... |
... |
@@ -798,10 +798,10 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
798 |
798 |
i <- which(rowSums(!is.finite(W)) > 0)
|
799 |
799 |
##browser()
|
800 |
800 |
stop("Possible zeros in the within-genotype estimates of the spread (vA, vB). ")
|
801 |
|
- }
|
802 |
|
- Ns <- tmp.objects[["Ns"]]
|
|
801 |
+ }
|
|
802 |
+ Ns <- tmp.objects[["Ns"]]
|
803 |
803 |
Ns <- Ns[I, ]
|
804 |
|
-
|
|
804 |
+
|
805 |
805 |
CHR <- unique(chromosome(object))
|
806 |
806 |
##batch <- unique(object$batch)
|
807 |
807 |
batch <- unique(batch(object))
|
... |
... |
@@ -812,10 +812,10 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
812 |
812 |
phiA <- getParam(object, "phiA", batch)
|
813 |
813 |
phiB <- getParam(object, "phiB", batch)
|
814 |
814 |
phiA.se <- getParam(object, "phiA.se", batch)
|
815 |
|
- phiB.se <- getParam(object, "phiB.se", batch)
|
|
815 |
+ phiB.se <- getParam(object, "phiB.se", batch)
|
816 |
816 |
if(CHR==23){
|
817 |
817 |
phiAX <- getParam(object, "phiAX", batch)
|
818 |
|
- phiBX <- getParam(object, "phiBX", batch)
|
|
818 |
+ phiBX <- getParam(object, "phiBX", batch)
|
819 |
819 |
}
|
820 |
820 |
NOHET <- mean(Ns[, "AB"], na.rm=TRUE) < 0.05
|
821 |
821 |
if(missing(allele)) stop("must specify allele")
|
... |
... |
@@ -830,7 +830,7 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
830 |
830 |
}
|
831 |
831 |
} else {##autosome
|
832 |
832 |
if(allele == "A") X <- cbind(1, 2:0) else X <- cbind(1, 0:2)
|
833 |
|
- if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
|
833 |
+ if(NOHET) X <- X[-2, ] ##more than 1 X chromosome, but all homozygous
|
834 |
834 |
}
|
835 |
835 |
##How to quickly generate Xstar, Xstar = diag(W) %*% X
|
836 |
836 |
Xstar <- apply(W, 1, generateX, X)
|
... |
... |
@@ -838,7 +838,7 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
838 |
838 |
##as.numeric(diag(W[1, ]) %*% X)
|
839 |
839 |
if(CHR == 23){
|
840 |
840 |
betahat <- matrix(NA, 3, nrow(Ystar))
|
841 |
|
- ses <- matrix(NA, 3, nrow(Ystar))
|
|
841 |
+ ses <- matrix(NA, 3, nrow(Ystar))
|
842 |
842 |
} else{
|
843 |
843 |
betahat <- matrix(NA, 2, nrow(Ystar))
|
844 |
844 |
ses <- matrix(NA, 2, nrow(Ystar))
|
... |
... |
@@ -859,7 +859,7 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
859 |
859 |
object <- pr(object, "phiA.se", batch, phiA.se)
|
860 |
860 |
if(CHR == 23){
|
861 |
861 |
phiAX[complete] <- betahat[3, ]
|
862 |
|
- object <- pr(object, "phiAX", batch, phiAX)
|
|
862 |
+ object <- pr(object, "phiAX", batch, phiAX)
|
863 |
863 |
}
|
864 |
864 |
}
|
865 |
865 |
if(allele=="B"){
|
... |
... |
@@ -869,7 +869,7 @@ nuphiAllele <- function(object, allele, Ystar, W, tmp.objects, cnOptions){
|
869 |
869 |
phiB.se[complete] <- ses[2, ]
|
870 |
870 |
object <- pr(object, "nuB", batch, nuB)
|
871 |
871 |
object <- pr(object, "nuB.se", batch, nuB.se)
|
872 |
|
- object <- pr(object, "phiB", batch, phiB)
|
|
872 |
+ object <- pr(object, "phiB", batch, phiB)
|
873 |
873 |
object <- pr(object, "phiB.se", batch, phiB.se)
|
874 |
874 |
if(CHR == 23){
|
875 |
875 |
phiBX[complete] <- betahat[3, ]
|
... |
... |
@@ -893,7 +893,7 @@ predictGender <- function(res, cdfName="genomewidesnp6", SNRMin=5){
|
893 |
893 |
loader("23index.rda", .crlmmPkgEnv, pkgname)
|
894 |
894 |
index <- getVarInEnv("index")
|
895 |
895 |
##load(file.path(path, "23index.rda"))
|
896 |
|
- XIndex <- index[[1]]
|
|
896 |
+ XIndex <- index[[1]]
|
897 |
897 |
if(class(res) != "ABset"){
|
898 |
898 |
A <- res$A
|
899 |
899 |
B <- res$B
|
... |
... |
@@ -1013,7 +1013,7 @@ crlmmCopynumber <- function(object,
|
1013 |
1013 |
column <- match(batchName, colnames(lM(object)[[k]]))
|
1014 |
1014 |
lM(object)[[k]][row.index, column] <- fData(tmp)[, k]
|
1015 |
1015 |
}
|
1016 |
|
- }
|
|
1016 |
+ }
|
1017 |
1017 |
rm(tmp); gc()
|
1018 |
1018 |
ii <- ii+1
|
1019 |
1019 |
}
|
... |
... |
@@ -1053,7 +1053,7 @@ crlmmCopynumber2 <- function(object,
|
1053 |
1053 |
##autosomeIndex.snps <- (1:nrow(object))[chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object))]
|
1054 |
1054 |
autosomeIndex.snps <- which(chromosome(object) < 23 & isSnp(object) & !is.na(chromosome(object)))
|
1055 |
1055 |
autosomeIndex.nps <- which(chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object)))
|
1056 |
|
- ##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]
|
|
1056 |
+ ##autosomeIndex.nps <- (1:nrow(object))[chromosome(object) < 23 & !isSnp(object) & !is.na(chromosome(object))]
|
1057 |
1057 |
|
1058 |
1058 |
## Do chromosome X in batches
|
1059 |
1059 |
Ns <- initializeBigMatrix("Ns", nrow(object), 5)
|
... |
... |
@@ -1069,8 +1069,8 @@ crlmmCopynumber2 <- function(object,
|
1069 |
1069 |
save(snpflags, file=file.path(ldPath(), "snpflags.rda"))
|
1070 |
1070 |
} else{
|
1071 |
1071 |
load(file.path(ldPath(), "snpflags.rda"))
|
1072 |
|
-
|
1073 |
|
- }
|
|
1072 |
+
|
|
1073 |
+ }
|
1074 |
1074 |
if(verbose) message("Estimating allele-specific copy number at autosomal SNPs")
|
1075 |
1075 |
snpBatches <- splitIndicesByLength(autosomeIndex.snps, ocProbesets())
|
1076 |
1076 |
ocLapply(seq(along=snpBatches),
|
... |
... |
@@ -1181,7 +1181,7 @@ fit.lm1 <- function(idxBatch,
|
1181 |
1181 |
snps <- snpBatches[[idxBatch]]
|
1182 |
1182 |
batches <- split(seq(along=batch(object)), batch(object))
|
1183 |
1183 |
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
|
1184 |
|
-
|
|
1184 |
+
|
1185 |
1185 |
open(object)
|
1186 |
1186 |
open(snpflags)
|
1187 |
1187 |
open(normal)
|
... |
... |
@@ -1312,14 +1312,14 @@ fit.lm1 <- function(idxBatch,
|
1312 |
1312 |
cB <- matrix(as.integer(cB*100), nrow(cB), ncol(cB))
|
1313 |
1313 |
|
1314 |
1314 |
|
1315 |
|
-
|
|
1315 |
+
|
1316 |
1316 |
CA(object)[snps, ] <- cA
|
1317 |
1317 |
CB(object)[snps, ] <- cB
|
1318 |
1318 |
|
1319 |
|
-
|
|
1319 |
+
|
1320 |
1320 |
snpflags[snps, ] <- flags
|
1321 |
1321 |
lapply(lM(object), open)
|
1322 |
|
-
|
|
1322 |
+
|
1323 |
1323 |
tmp <- physical(lM(object))$tau2A
|
1324 |
1324 |
tmp[snps, ] <- tau2A
|
1325 |
1325 |
lM(object)$tau2A <- tmp
|
... |
... |
@@ -1356,7 +1356,7 @@ fit.lm1 <- function(idxBatch,
|
1356 |
1356 |
tmp <- physical(lM(object))$corrBB
|
1357 |
1357 |
tmp[snps, ] <- corrBB
|
1358 |
1358 |
lM(object)$corrAB <- tmp
|
1359 |
|
-
|
|
1359 |
+
|
1360 |
1360 |
lapply(assayData(object), close)
|
1361 |
1361 |
lapply(lM(object), close)
|
1362 |
1362 |
TRUE
|
... |
... |
@@ -1384,15 +1384,15 @@ fit.lm2 <- function(idxBatch,
|
1384 |
1384 |
snps <- snpBatches[[idxBatch]]
|
1385 |
1385 |
batches <- split(seq(along=batch(object)), batch(object))
|
1386 |
1386 |
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
|
1387 |
|
-
|
|
1387 |
+
|
1388 |
1388 |
open(object)
|
1389 |
1389 |
open(snpflags)
|
1390 |
1390 |
open(normal)
|
1391 |
1391 |
|
1392 |
|
-
|
|
1392 |
+
|
1393 |
1393 |
cA <- matrix(NA, length(snps), ncol(object))
|
1394 |
1394 |
ii <- isSnp(object) & chromosome(object) < 23 & !is.na(chromosome(object))
|
1395 |
|
- flags <- as.matrix(snpflags[,])
|
|
1395 |
+ flags <- as.matrix(snpflags[,])
|
1396 |
1396 |
noflags <- rowSums(flags, na.rm=TRUE) == 0 ##NA's for unevaluated batches
|
1397 |
1397 |
## We do not want to write to discuss for each batch. More efficient to
|
1398 |
1398 |
## write to disk after estimating these parameters for all batches.
|
... |
... |
@@ -1501,7 +1501,7 @@ fit.lm3 <- function(idxBatch,
|
1501 |
1501 |
snps <- snpBatches[[idxBatch]]
|
1502 |
1502 |
batches <- split(seq(along=batch(object)), batch(object))
|
1503 |
1503 |
batches <- batches[sapply(batches, length) >= MIN.SAMPLES]
|
1504 |
|
-
|
|
1504 |
+
|
1505 |
1505 |
open(snpflags)
|
1506 |
1506 |
open(normal)
|
1507 |
1507 |
open(object)
|
... |
... |
@@ -1519,7 +1519,7 @@ fit.lm3 <- function(idxBatch,
|
1519 |
1519 |
AA <- as.matrix(A(object)[snps, ])
|
1520 |
1520 |
BB <- as.matrix(B(object)[snps, ])
|
1521 |
1521 |
for(k in batches){
|
1522 |
|
- ##if(verbose) message("SNP batch ", ii, " of ", length(batches))
|
|
1522 |
+ ##if(verbose) message("SNP batch ", ii, " of ", length(batches))
|
1523 |
1523 |
## within-genotype moments
|
1524 |
1524 |
gender <- object$gender[k]
|
1525 |
1525 |
G <- GG[, k]
|
... |
... |
@@ -1540,7 +1540,7 @@ fit.lm3 <- function(idxBatch,
|
1540 |
1540 |
vA.F <- applyByGenotype(A[, gender==2], rowMAD, G[, gender==2])
|
1541 |
1541 |
vB.F <- applyByGenotype(B[, gender==2], rowMAD, G[, gender==2])
|
1542 |
1542 |
vA.M <- applyByGenotype(A[, gender==1], rowMAD, G[, gender==1])
|
1543 |
|
- vB.M <- applyByGenotype(B[, gender==1], rowMAD, G[, gender==1])
|
|
1543 |
+ vB.M <- applyByGenotype(B[, gender==1], rowMAD, G[, gender==1])
|
1544 |
1544 |
vA.F <- shrink(vA.F, Ns.F, DF.PRIOR)
|
1545 |
1545 |
vA.M <- shrink(vA.M, Ns.M, DF.PRIOR)
|
1546 |
1546 |
vB.F <- shrink(vB.F, Ns.F, DF.PRIOR)
|
... |
... |
@@ -1583,7 +1583,7 @@ fit.lm3 <- function(idxBatch,
|
1583 |
1583 |
notMissing <- !(is.na(muA.M[, 1]) | is.na(muA.M[, 3]) | is.na(muB.M[, 1]) | is.na(muB.M[, 3]))
|
1584 |
1584 |
complete <- list()
|
1585 |
1585 |
complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
|
1586 |
|
- complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
|
|
1586 |
+ complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
|
1587 |
1587 |
size <- min(5000, length(complete[[1]]))
|
1588 |
1588 |
if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
|
1589 |
1589 |
##
|
... |
... |
@@ -1658,7 +1658,7 @@ fit.lm3 <- function(idxBatch,
|
1658 |
1658 |
cB[cB < 0.05] <- 0.05
|
1659 |
1659 |
cA[cA > 5] <- 5
|
1660 |
1660 |
cB[cB > 5] <- 5
|
1661 |
|
-
|
|
1661 |
+
|
1662 |
1662 |
##--------------------------------------------------
|
1663 |
1663 |
##RS: need to fix. why are there NA's by coercion
|
1664 |
1664 |
cA <- matrix(as.integer(cA*100), nrow(cA), ncol(cA))
|
... |
... |
@@ -1702,7 +1702,7 @@ fit.lm3 <- function(idxBatch,
|
1702 |
1702 |
tmp <- physical(lM(object))$phiPrimeB
|
1703 |
1703 |
tmp[snps, ] <- phiB2
|
1704 |
1704 |
lM(object)$phiPrimeB <- tmp
|
1705 |
|
-
|
|
1705 |
+
|
1706 |
1706 |
tmp <- physical(lM(object))$corrAB
|
1707 |
1707 |
tmp[snps, ] <- corrAB
|
1708 |
1708 |
lM(object)$corrAB <- tmp
|
... |
... |
@@ -1711,7 +1711,7 @@ fit.lm3 <- function(idxBatch,
|
1711 |
1711 |
lM(object)$corrAA <- tmp
|
1712 |
1712 |
tmp <- physical(lM(object))$corrBB
|
1713 |
1713 |
tmp[snps, ] <- corrBB
|
1714 |
|
- lM(object)$corrAB <- tmp
|
|
1714 |
+ lM(object)$corrAB <- tmp
|
1715 |
1715 |
## lM(object)$tau2A[snps, ] <- tau2A
|
1716 |
1716 |
## lM(object)$tau2B[snps, ] <- tau2B
|
1717 |
1717 |
## lM(object)$sig2A[snps, ] <- sig2A
|
... |
... |
@@ -1744,7 +1744,7 @@ fit.lm4 <- function(idxBatch,
|
1744 |
1744 |
MIN.NU,
|
1745 |
1745 |
MIN.PHI,
|
1746 |
1746 |
verbose, ...){
|
1747 |
|
- if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
|
|
1747 |
+ if(verbose) message("Probe batch ", idxBatch, " of ", length(snpBatches))
|
1748 |
1748 |
open(object)
|
1749 |
1749 |
open(normal)
|
1750 |
1750 |
open(snpflags)
|
... |
... |
@@ -1766,7 +1766,7 @@ fit.lm4 <- function(idxBatch,
|
1766 |
1766 |
i2 <- rowSums(nuIB < 20, na.rm=TRUE) == 0
|
1767 |
1767 |
i3 <- rowSums(phiIA < 20, na.rm=TRUE) == 0
|
1768 |
1768 |
i4 <- rowSums(phiIB < 20, na.rm=TRUE) == 0
|
1769 |
|
-
|
|
1769 |
+
|
1770 |
1770 |
snp.index <- which(i1 & i2 & i3 & i4 & noflags)
|
1771 |
1771 |
if(length(snp.index) == 0){
|
1772 |
1772 |
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")
|
... |
... |
@@ -1793,7 +1793,7 @@ fit.lm4 <- function(idxBatch,
|
1793 |
1793 |
##if(missing(which.batches)) which.batches <- seq(along=batches)
|
1794 |
1794 |
##batches <- batches[which.batches]
|
1795 |
1795 |
for(k in batches){
|
1796 |
|
- ##if(verbose) message("SNP batch ", ii, " of ", length(batches))
|
|
1796 |
+ ##if(verbose) message("SNP batch ", ii, " of ", length(batches))
|
1797 |
1797 |
G <- GG[, k]
|
1798 |
1798 |
xx <- CP[, k]
|
1799 |
1799 |
highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
|
... |
... |
@@ -1802,7 +1802,7 @@ fit.lm4 <- function(idxBatch,
|
1802 |
1802 |
AA <- A.snp[, k]
|
1803 |
1803 |
BB <- B.snp[, k]
|
1804 |
1804 |
|
1805 |
|
-
|
|
1805 |
+
|
1806 |
1806 |
##index <- GT.B <- GT.A <- vector("list", 3)
|
1807 |
1807 |
##names(index) <- names(GT.B) <- names(GT.A) <- c("AA", "AB", "BB")
|
1808 |
1808 |
Ns <- applyByGenotype(matrix(1, nrow(G), ncol(G)), rowSums, G)
|
... |
... |
@@ -1812,7 +1812,7 @@ fit.lm4 <- function(idxBatch,
|
1812 |
1812 |
muB <- muB[, 3]
|
1813 |
1813 |
X <- cbind(1, log2(c(muA, muB)))
|
1814 |
1814 |
J <- match(unique(batch(object)[k]), unique(batch(object)))
|
1815 |
|
-
|
|
1815 |
+
|
1816 |
1816 |
Y <- log2(c(phiA.snp[, J], phiB.snp[, J]))
|
1817 |
1817 |
|
1818 |
1818 |
##--------------------------------------------------
|
... |
... |
@@ -1822,7 +1822,7 @@ fit.lm4 <- function(idxBatch,
|
1822 |
1822 |
X <- X[!remove, ]
|
1823 |
1823 |
##--------------------------------------------------
|
1824 |
1824 |
betahat <- solve(crossprod(X), crossprod(X, Y))
|
1825 |
|
-
|
|
1825 |
+
|
1826 |
1826 |
|
1827 |
1827 |
##nonpolymorphic
|
1828 |
1828 |
A <- AA.np[, k]
|
... |
... |
@@ -1917,7 +1917,7 @@ cnrma2 <- function(A, filenames, row.names, verbose=TRUE, seed=1, cdfName, sns){
|
1917 |
1917 |
row.names=row.names,
|
1918 |
1918 |
filenames=filenames,
|
1919 |
1919 |
A=A,
|
1920 |
|
- SKW=SKW,
|
|
1920 |
+ SKW=SKW,
|
1921 |
1921 |
seed=seed,
|
1922 |
1922 |
pkgname=pkgname,
|
1923 |
1923 |
cdfName=cdfName,
|
... |
... |
@@ -2051,15 +2051,15 @@ lm.parameters <- function(object, batch){##cnOptions){
|
2051 |
2051 |
paste("sig2A", uplate, sep="_"),
|
2052 |
2052 |
paste("sig2B", uplate, sep="_"),
|
2053 |
2053 |
paste("nuA", uplate, sep="_"),
|
2054 |
|
- paste("nuA.se", uplate, sep="_"),
|
|
2054 |
+ paste("nuA.se", uplate, sep="_"),
|
2055 |
2055 |
paste("nuB", uplate, sep="_"),
|
2056 |
|
- paste("nuB.se", uplate, sep="_"),
|
|
2056 |
+ paste("nuB.se", uplate, sep="_"),
|
2057 |
2057 |
paste("phiA", uplate, sep="_"),
|
2058 |
|
- paste("phiA.se", uplate, sep="_"),
|
|
2058 |
+ paste("phiA.se", uplate, sep="_"),
|
2059 |
2059 |
paste("phiB", uplate, sep="_"),
|
2060 |
|
- paste("phiB.se", uplate, sep="_"),
|
|
2060 |
+ paste("phiB.se", uplate, sep="_"),
|
2061 |
2061 |
paste("phiAX", uplate, sep="_"),
|
2062 |
|
- paste("phiBX", uplate, sep="_"),
|
|
2062 |
+ paste("phiBX", uplate, sep="_"),
|
2063 |
2063 |
paste("corr", uplate, sep="_"),
|
2064 |
2064 |
paste("corrA.BB", uplate, sep="_"),
|
2065 |
2065 |
paste("corrB.AA", uplate, sep="_"))
|
... |
... |
@@ -2098,7 +2098,7 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
|
2098 |
2098 |
nuA <- getParam(object, "nuA", batch)
|
2099 |
2099 |
nuB <- getParam(object, "nuB", batch)
|
2100 |
2100 |
phiA <- getParam(object, "phiA", batch)
|
2101 |
|
- phiB <- getParam(object, "phiB", batch)
|
|
2101 |
+ phiB <- getParam(object, "phiB", batch)
|
2102 |
2102 |
sns <- sampleNames(object)
|
2103 |
2103 |
muA <- tmp.objects[["muA"]]
|
2104 |
2104 |
muB <- tmp.objects[["muB"]]
|
... |
... |
@@ -2128,13 +2128,13 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
|
2128 |
2128 |
ix <- 1:nrow(X)
|
2129 |
2129 |
}
|
2130 |
2130 |
betahat <- solve(crossprod(X[ix, ]), crossprod(X[ix, ], Y[ix]))
|
2131 |
|
- normal <- tmp.objects[["normal"]][!isSnp(object), ]
|
|
2131 |
+ normal <- tmp.objects[["normal"]][!isSnp(object), ]
|
2132 |
2132 |
if(CHR == 23){
|
2133 |
2133 |
##normalNP <- envir[["normalNP"]]
|
2134 |
2134 |
##normalNP <- normalNP[, plate==uplate[p]]
|
2135 |
2135 |
##nuT <- envir[["nuT"]]
|
2136 |
2136 |
##phiT <- envir[["phiT"]]
|
2137 |
|
-
|
|
2137 |
+
|
2138 |
2138 |
##cnvs <- envir[["cnvs"]]
|
2139 |
2139 |
##loader("cnProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
|
2140 |
2140 |
##cnProbes <- get("cnProbes", envir=.crlmmPkgEnv)
|
... |
... |
@@ -2162,7 +2162,7 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
|
2162 |
2162 |
mus <- log2(cbind(mu1, mu2))
|
2163 |
2163 |
X.men <- cbind(1, mus[, 1])
|
2164 |
2164 |
X.fem <- cbind(1, mus[, 2])
|
2165 |
|
-
|
|
2165 |
+
|
2166 |
2166 |
Yhat1 <- as.numeric(X.men %*% betahat)
|
2167 |
2167 |
Yhat2 <- as.numeric(X.fem %*% betahat)
|
2168 |
2168 |
phi1 <- 2^(Yhat1)
|
... |
... |
@@ -2199,17 +2199,17 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
|
2199 |
2199 |
|
2200 |
2200 |
nuA[!isSnp(object)] <- nu2
|
2201 |
2201 |
phiA[!isSnp(object)] <- phi2
|
2202 |
|
-
|
|
2202 |
+
|
2203 |
2203 |
THR.NU.PHI <- cnOptions$THR.NU.PHI
|
2204 |
2204 |
if(THR.NU.PHI){
|
2205 |
2205 |
verbose <- cnOptions$verbose
|
2206 |
2206 |
##Assign values to object
|
2207 |
2207 |
object <- pr(object, "nuA", batch, nuA)
|
2208 |
|
- object <- pr(object, "phiA", batch, phiA)
|
|
2208 |
+ object <- pr(object, "phiA", batch, phiA)
|
2209 |
2209 |
##if(verbose) message("Thresholding nu and phi")
|
2210 |
2210 |
object <- thresholdModelParams(object, cnOptions)
|
2211 |
2211 |
} else {
|
2212 |
|
- object <- pr(object, "nuA", batch, nuA)
|
|
2212 |
+ object <- pr(object, "nuA", batch, nuA)
|
2213 |
2213 |
object <- pr(object, "phiA", batch, phiA)
|
2214 |
2214 |
}
|
2215 |
2215 |
} else {
|
... |
... |
@@ -2226,7 +2226,7 @@ nonpolymorphic <- function(object, cnOptions, tmp.objects){
|
2226 |
2226 |
verbose <- cnOptions$verbose
|
2227 |
2227 |
##Assign values to object
|
2228 |
2228 |
object <- pr(object, "nuA", batch, nuA)
|
2229 |
|
- object <- pr(object, "phiA", batch, phiA)
|
|
2229 |
+ object <- pr(object, "phiA", batch, phiA)
|
2230 |
2230 |
##if(verbose) message("Thresholding nu and phi")
|
2231 |
2231 |
object <- thresholdModelParams(object, cnOptions)
|
2232 |
2232 |
##reassign values (now thresholded at MIN.NU and MIN.PHI
|
... |
... |
@@ -2254,7 +2254,7 @@ withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
|
2254 |
2254 |
vA <- tmp.objects[["vA"]]
|
2255 |
2255 |
vB <- tmp.objects[["vB"]]
|
2256 |
2256 |
Ns <- tmp.objects[["Ns"]]
|
2257 |
|
- G <- snpCall(object)
|
|
2257 |
+ G <- snpCall(object)
|
2258 |
2258 |
GT.CONF.THR <- cnOptions$GT.CONF.THR
|
2259 |
2259 |
CHR <- unique(chromosome(object))
|
2260 |
2260 |
A <- A(object)
|
... |
... |
@@ -2281,11 +2281,11 @@ withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
|
2281 |
2281 |
for(j in 1:3){
|
2282 |
2282 |
GT <- G==j & highConf & IX & snpIndicator
|
2283 |
2283 |
GT <- GT * normal
|
2284 |
|
- Ns[, j+2] <- rowSums(GT, na.rm=TRUE)
|
|
2284 |
+ Ns[, j+2] <- rowSums(GT, na.rm=TRUE)
|
2285 |
2285 |
GT[GT == FALSE] <- NA
|
2286 |
2286 |
GT.A[[j]] <- GT*A
|
2287 |
2287 |
GT.B[[j]] <- GT*B
|
2288 |
|
- index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added
|
|
2288 |
+ index[[j]] <- which(Ns[, j+2] > 0 & isSnp(object)) ##RS: added
|
2289 |
2289 |
muA[, j+2] <- rowMedians(GT.A[[j]], na.rm=TRUE)
|
2290 |
2290 |
muB[, j+2] <- rowMedians(GT.B[[j]], na.rm=TRUE)
|
2291 |
2291 |
vA[, j+2] <- rowMAD(GT.A[[j]], na.rm=TRUE)
|
... |
... |
@@ -2307,14 +2307,14 @@ withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
|
2307 |
2307 |
if(CHR == 23){
|
2308 |
2308 |
k <- 1
|
2309 |
2309 |
for(j in c(1,3)){
|
2310 |
|
- GT <- G==j & highConf & !IX
|
|
2310 |
+ GT <- G==j & highConf & !IX
|
2311 |
2311 |
Ns[, k] <- rowSums(GT)
|
2312 |
2312 |
GT[GT == FALSE] <- NA
|
2313 |
2313 |
muA[, k] <- rowMedians(GT*A, na.rm=TRUE)
|
2314 |
2314 |
muB[, k] <- rowMedians(GT*B, na.rm=TRUE)
|
2315 |
2315 |
vA[, k] <- rowMAD(GT*A, na.rm=TRUE)
|
2316 |
2316 |
vB[, k] <- rowMAD(GT*B, na.rm=TRUE)
|
2317 |
|
-
|
|
2317 |
+
|
2318 |
2318 |
DF <- Ns[, k]-1
|
2319 |
2319 |
DF[DF < 1] <- 1
|
2320 |
2320 |
v0A <- median(vA[, k], na.rm=TRUE)
|
... |
... |
@@ -2322,7 +2322,7 @@ withinGenotypeMoments <- function(object, cnOptions, tmp.objects){
|
2322 |
2322 |
vA[, k] <- (vA[, k]*DF + v0A*DF.PRIOR)/(DF.PRIOR+DF)
|
2323 |
2323 |
vA[is.na(vA[, k]), k] <- v0A
|
2324 |
2324 |
vB[, k] <- (vB[, k]*DF + v0B*DF.PRIOR)/(DF.PRIOR+DF)
|
2325 |
|
- vB[is.na(vB[, k]), k] <- v0B
|
|
2325 |
+ vB[is.na(vB[, k]), k] <- v0B
|
2326 |
2326 |
k <- k+1
|
2327 |
2327 |
}
|
2328 |
2328 |
}
|
... |
... |
@@ -2394,7 +2394,7 @@ oneBatch <- function(object, cnOptions, tmp.objects){
|
2394 |
2394 |
index.BB <- which(Ns[, "BB"] >= MIN.OBS)
|
2395 |
2395 |
correct.orderA <- muA[, "AA"] > muA[, "BB"]
|
2396 |
2396 |
correct.orderB <- muB[, "BB"] > muB[, "AA"]
|
2397 |
|
- ##For chr X, this will ignore the males
|
|
2397 |
+ ##For chr X, this will ignore the males
|
2398 |
2398 |
nobs <- rowSums(Ns[, 3:5] >= MIN.OBS, na.rm=TRUE) == 3
|
2399 |
2399 |
index.complete <- which(correct.orderA & correct.orderB & nobs) ##be selective here
|
2400 |
2400 |
size <- min(5000, length(index.complete))
|
... |
... |
@@ -2423,7 +2423,7 @@ oneBatch <- function(object, cnOptions, tmp.objects){
|
2423 |
2423 |
notMissing <- !(is.na(muA[, "A"]) | is.na(muA[, "B"]) | is.na(muB[, "A"]) | is.na(muB[, "B"]))
|
2424 |
2424 |
complete <- list()
|
2425 |
2425 |
complete[[1]] <- which(correct.orderA & correct.orderB & nobsA & notMissing) ##be selective here
|
2426 |
|
- complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
|
|
2426 |
+ complete[[2]] <- which(correct.orderA & correct.orderB & nobsB & notMissing) ##be selective here
|
2427 |
2427 |
size <- min(5000, length(complete[[1]]))
|
2428 |
2428 |
if(size > 5000) complete <- lapply(complete, function(x) sample(x, size))
|
2429 |
2429 |
if(CHR == 23){
|
... |
... |
@@ -2473,7 +2473,7 @@ oneBatch <- function(object, cnOptions, tmp.objects){
|
2473 |
2473 |
muB[index[[j]], -c(1, 2, k+2)] <- mus[, 3:4]
|
2474 |
2474 |
}
|
2475 |
2475 |
negA <- rowSums(muA < 0) > 0
|
2476 |
|
- negB <- rowSums(muB < 0) > 0
|
|
2476 |
+ negB <- rowSums(muB < 0) > 0
|
2477 |
2477 |
snpflags <- snpflags| negA | negB | rowSums(is.na(muA[, 3:5]), na.rm=TRUE) > 0
|
2478 |
2478 |
tmp.objects$snpflags <- snpflags
|
2479 |
2479 |
tmp.objects[["muA"]] <- muA
|
... |
... |
@@ -2497,7 +2497,7 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2497 |
2497 |
AA.A <- GT.A[[1]]
|
2498 |
2498 |
AB.A <- GT.A[[2]]
|
2499 |
2499 |
BB.A <- GT.A[[3]]
|
2500 |
|
-
|
|
2500 |
+
|
2501 |
2501 |
AA.B <- GT.B[[1]]
|
2502 |
2502 |
AB.B <- GT.B[[2]]
|
2503 |
2503 |
BB.B <- GT.B[[3]]
|
... |
... |
@@ -2505,7 +2505,7 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2505 |
2505 |
##batch <- unique(object$batch)
|
2506 |
2506 |
batch <- unique(batch(object))
|
2507 |
2507 |
tau2A <- getParam(object, "tau2A", batch)
|
2508 |
|
- tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
|
2508 |
+ tau2A[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
2509 |
2509 |
DF <- Ns[, "BB"]-1
|
2510 |
2510 |
DF[DF < 1] <- 1
|
2511 |
2511 |
med <- median(tau2A, na.rm=TRUE)
|
... |
... |
@@ -2513,17 +2513,17 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2513 |
2513 |
tau2A[is.na(tau2A) & isSnp(object)] <- med
|
2514 |
2514 |
object <- pr(object, "tau2A", batch, tau2A)
|
2515 |
2515 |
|
2516 |
|
- sig2B <- getParam(object, "sig2B", batch)
|
|
2516 |
+ sig2B <- getParam(object, "sig2B", batch)
|
2517 |
2517 |
x <- BB.B[index.BB, ]
|
2518 |
|
- sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
|
2518 |
+ sig2B[index.BB] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
2519 |
2519 |
med <- median(sig2B, na.rm=TRUE)
|
2520 |
2520 |
sig2B <- (sig2B * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
|
2521 |
2521 |
sig2B[is.na(sig2B) & isSnp(object)] <- med
|
2522 |
|
- object <- pr(object, "sig2B", batch, sig2B)
|
|
2522 |
+ object <- pr(object, "sig2B", batch, sig2B)
|
2523 |
2523 |
|
2524 |
|
- tau2B <- getParam(object, "tau2B", batch)
|
|
2524 |
+ tau2B <- getParam(object, "tau2B", batch)
|
2525 |
2525 |
x <- AA.B[index.AA, ]
|
2526 |
|
- tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
|
2526 |
+ tau2B[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2
|
2527 |
2527 |
DF <- Ns[, "AA"]-1
|
2528 |
2528 |
DF[DF < 1] <- 1
|
2529 |
2529 |
med <- median(tau2B, na.rm=TRUE)
|
... |
... |
@@ -2531,9 +2531,9 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2531 |
2531 |
tau2B[is.na(tau2B) & isSnp(object)] <- med
|
2532 |
2532 |
object <- pr(object, "tau2B", batch, tau2B)
|
2533 |
2533 |
|
2534 |
|
- sig2A <- getParam(object, "sig2A", batch)
|
|
2534 |
+ sig2A <- getParam(object, "sig2A", batch)
|
2535 |
2535 |
x <- AA.A[index.AA, ]
|
2536 |
|
- sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)
|
|
2536 |
+ sig2A[index.AA] <- rowMAD(log2(x), log2(x), na.rm=TRUE)^2##var(log(IA)|AA)
|
2537 |
2537 |
med <- median(sig2A, na.rm=TRUE)
|
2538 |
2538 |
sig2A <- (sig2A * DF + med * DF.PRIOR)/(DF.PRIOR + DF)
|
2539 |
2539 |
sig2A[is.na(sig2A) & isSnp(object)] <- med
|
... |
... |
@@ -2550,7 +2550,7 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2550 |
2550 |
med <- median(corr, na.rm=TRUE)
|
2551 |
2551 |
corr <- (corr*DF + med * DF.PRIOR)/(DF.PRIOR + DF)
|
2552 |
2552 |
corr[is.na(corr) & isSnp(object)] <- med
|
2553 |
|
- object <- pr(object, "corr", batch, corr)
|
|
2553 |
+ object <- pr(object, "corr", batch, corr)
|
2554 |
2554 |
}
|
2555 |
2555 |
corrB.AA <- getParam(object, "corrB.AA", batch)
|
2556 |
2556 |
backgroundB <- AA.B[index.AA, ]
|
... |
... |
@@ -2563,7 +2563,7 @@ locationAndScale <- function(object, cnOptions, tmp.objects){
|
2563 |
2563 |
corrB.AA[is.na(corrB.AA) & isSnp(object)] <- med
|
2564 |
2564 |
object <- pr(object, "corrB.AA", batch, corrB.AA)
|
2565 |
2565 |
|
2566 |
|
- corrA.BB <- getParam(object, "corrA.BB", batch)
|
|
2566 |
+ corrA.BB <- getParam(object, "corrA.BB", batch)
|
2567 |
2567 |
backgroundA <- BB.A[index.BB, ]
|
2568 |
2568 |
signalB <- BB.B[index.BB, ]
|
2569 |
2569 |
corrA.BB[index.BB] <- rowCors(backgroundA, signalB, na.rm=TRUE)
|
... |
... |
@@ -2598,7 +2598,7 @@ coefs <- function(object, cnOptions, tmp.objects){
|
2598 |
2598 |
IB <- muB[, -4]
|
2599 |
2599 |
vA <- vA[, -4]
|
2600 |
2600 |
vB <- vB[, -4]
|
2601 |
|
- Np <- Ns[, -4]
|
|
2601 |
+ Np <- Ns[, -4]
|
2602 |
2602 |
} else{
|
2603 |
2603 |
IA <- muA
|
2604 |
2604 |
IB <- muB
|
... |
... |
@@ -2606,7 +2606,7 @@ coefs <- function(object, cnOptions, tmp.objects){
|
2606 |
2606 |
vB <- vB
|
2607 |
2607 |
Np <- Ns
|
2608 |
2608 |
}
|
2609 |
|
-
|
|
2609 |
+
|
2610 |
2610 |
}
|
2611 |
2611 |
Np[Np < 1] <- 1
|
2612 |
2612 |
vA2 <- vA^2/Np
|
... |
... |
@@ -2617,7 +2617,7 @@ coefs <- function(object, cnOptions, tmp.objects){
|
2617 |
2617 |
YB <- IB*wB
|
2618 |
2618 |
##update lm.coefficients stored in object
|
2619 |
2619 |
object <- nuphiAllele(object, allele="A", Ystar=YA, W=wA, tmp.objects=tmp.objects, cnOptions=cnOptions)
|
2620 |
|
- object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)
|
|
2620 |
+ object <- nuphiAllele(object, allele="B", Ystar=YB, W=wB, tmp.objects=tmp.objects, cnOptions=cnOptions)
|
2621 |
2621 |
##---------------------------------------------------------------------------
|
2622 |
2622 |
##Estimate crosshyb using X chromosome and sequence information
|
2623 |
2623 |
##---------------------------------------------------------------------------
|
... |
... |
@@ -2644,7 +2644,7 @@ polymorphic <- function(object, cnOptions, tmp.objects){
|
2644 |
2644 |
vA <- tmp.objects[["vA"]]
|
2645 |
2645 |
vB <- tmp.objects[["vB"]]
|
2646 |
2646 |
Ns <- tmp.objects[["Ns"]]
|
2647 |
|
-
|
|
2647 |
+
|
2648 |
2648 |
nuA <- getParam(object, "nuA", batch)
|
2649 |
2649 |
nuB <- getParam(object, "nuB", batch)
|
2650 |
2650 |
nuA.se <- getParam(object, "nuA.se", batch)
|
... |
... |
@@ -2663,12 +2663,12 @@ polymorphic <- function(object, cnOptions, tmp.objects){
|
2663 |
2663 |
##---------------------------------------------------------------------------
|
2664 |
2664 |
if(CHR == 23){
|
2665 |
2665 |
phiAX <- getParam(object, "phiAX", batch) ##nonspecific hybridization coef
|
2666 |
|
- phiBX <- getParam(object, "phiBX", batch) ##nonspecific hybridization coef
|
2667 |
|
- phistar <- phiBX/phiA
|
|
2666 |
+ phiBX <- getParam(object, "phiBX", batch) ##nonspecific hybridization coef
|
|
2667 |
+ phistar <- phiBX/phiA
|
2668 |
2668 |
tmp <- (B-nuB - phistar*A + phistar*nuA)/phiB
|
2669 |
2669 |
copyB <- tmp/(1-phistar*phiAX/phiB)
|
2670 |
2670 |
copyA <- (A-nuA-phiAX*copyB)/phiA
|
2671 |
|
- CB(object) <- copyB ## multiplies by 100 and converts to integer
|
|
2671 |
+ CB(object) <- copyB ## multiplies by 100 and converts to integer
|
2672 |
2672 |
CA(object) <- copyA
|
2673 |
2673 |
} else{
|
2674 |
2674 |
CA(object) <- matrix((1/phiA*(A-nuA)), nrow(A), ncol(A))
|
... |
... |
@@ -2693,21 +2693,21 @@ posteriorProbability.snps <- function(object, cnOptions, tmp.objects=list()){
|
2693 |
2693 |
sig2A <- getParam(object, "sig2A", batch)[I]
|
2694 |
2694 |
sig2B <- getParam(object, "sig2B", batch)[I]
|
2695 |
2695 |
tau2A <- getParam(object, "tau2A", batch)[I]
|
2696 |
|
- tau2B <- getParam(object, "tau2B", batch)[I]
|
|
2696 |
+ tau2B <- getParam(object, "tau2B", batch)[I]
|
2697 |
2697 |
corrA.BB <- getParam(object, "corrA.BB", batch)[I]
|
2698 |
2698 |
corrB.AA <- getParam(object, "corrB.AA", batch)[I]
|
2699 |
|
- corr <- getParam(object, "corr", batch)[I]
|
|
2699 |
+ corr <- getParam(object, "corr", batch)[I]
|
2700 |
2700 |
nuA <- getParam(object, "nuA", batch)[I]
|
2701 |
|
- nuB <- getParam(object, "nuB", batch)[I]
|
|
2701 |
+ nuB <- getParam(object, "nuB", batch)[I]
|
2702 |
2702 |
phiA <- getParam(object, "phiA", batch)[I]
|
2703 |
2703 |
phiB <- getParam(object, "phiB", batch)[I]
|
2704 |
2704 |
normal <- tmp.objects[["normal"]][I, ]
|
2705 |
2705 |
prior.prob <- cnOptions$prior.prob
|
2706 |
|
- emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'
|
|
2706 |
+ emit <- array(NA, dim=c(nrow(A), ncol(A), 10))##SNPs x sample x 'truth'
|
2707 |
2707 |
lA <- log2(A)
|
2708 |
|
- lB <- log2(B)
|
|
2708 |
+ lB <- log2(B)
|
2709 |
2709 |
X <- cbind(lA, lB)
|
2710 |
|
- counter <- 1##state counter
|
|
2710 |
+ counter <- 1##state counter
|
2711 |
2711 |
for(CT in 0:3){
|
2712 |
2712 |
for(CA in 0:CT){
|
2713 |
2713 |
cat(".")
|
... |
... |
@@ -2721,7 +2721,7 @@ posteriorProbability.snps <- function(object, cnOptions, tmp.objects=list()){
|
2721 |
2721 |
if(CHR == 23){
|
2722 |
2722 |
##means <- cbind(suppressWarnings(log2(nuA[, p]+CA*phiA[, p] + CB*phiAx[, p])), suppressWarnings(log2(nuB[, p]+CB*phiB[, p] + CA*phiBx[, p])))
|
2723 |
2723 |
meanA <- suppressWarnings(log2(nuA+CA*phiA + CB*phiAX))
|
2724 |
|
- meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))
|
|
2724 |
+ meanB <- suppressWarnings(log2(nuB+CB*phiB + CA*phiBX))
|
2725 |
2725 |
} else{
|
2726 |
2726 |
##means <- cbind(suppressWarnings(log2(nuA+CA*phiA)), suppressWarnings(log2(nuB+CB*phiB)))
|
2727 |
2727 |
meanA <- suppressWarnings(log2(nuA+CA*phiA))
|
... |
... |
@@ -2803,7 +2803,7 @@ bias1 <- function(idxBatch,
|
2803 |
2803 |
prior.prob,
|
2804 |
2804 |
MIN.SAMPLES,
|
2805 |
2805 |
verbose){
|
2806 |
|
-
|
|
2806 |
+
|
2807 |
2807 |
}
|
2808 |
2808 |
|
2809 |
2809 |
bias2 <- function(idxBatch,
|
... |
... |
@@ -2854,7 +2854,7 @@ bias2 <- function(idxBatch,
|
2854 |
2854 |
|
2855 |
2855 |
index <- which(prUp > 0.05 & prUp > prDn)
|
2856 |
2856 |
##if proportion up greater than 5%, trim the high cn est.
|
2857 |
|
- norm[index, k] <- argmax.cn[index, ] > 3
|
|
2857 |
+ norm[index, k] <- argmax.cn[index, ] > 3
|
2858 |
2858 |
|
2859 |
2859 |
index <- which(prDn > 0.05 & prDn > prUp)
|
2860 |
2860 |
norm[index, k] <- argmax.cn[index, ] < 3
|
... |
... |
@@ -2897,7 +2897,7 @@ biasAdjust <- function(object, prior.prob=rep(1/4, 4), MIN.SAMPLES=10, verbose=T
|
2897 |
2897 |
prior.prob=prior.prob,
|
2898 |
2898 |
emit=emit,
|
2899 |
2899 |
MIN.SAMPLES=MIN.SAMPLES,
|
2900 |
|
- verbose=verbose)
|
|
2900 |
+ verbose=verbose)
|
2901 |
2901 |
}
|
2902 |
2902 |
|
2903 |
2903 |
|
... |
... |
@@ -2915,7 +2915,7 @@ biasAdjNP <- function(object, cnOptions, tmp.objects){
|
2915 |
2915 |
nuA <- getParam(object, "nuA", batch)
|
2916 |
2916 |
phiA <- getParam(object, "phiA", batch)
|
2917 |
2917 |
prior.prob <- cnOptions$prior.prob
|
2918 |
|
- emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'
|
|
2918 |
+ emit <- array(NA, dim=c(nrow(A), ncol(A), 4))##SNPs x sample x 'truth'
|
2919 |
2919 |
lT <- log2(A)
|
2920 |
2920 |
I <- isSnp(object)
|
2921 |
2921 |
counter <- 1 ##state counter
|
... |
... |
@@ -2936,7 +2936,7 @@ biasAdjNP <- function(object, cnOptions, tmp.objects){
|
2936 |
2936 |
counter <- counter+1
|
2937 |
2937 |
}
|
2938 |
2938 |
mostLikelyState <- apply(emit, c(1, 2), function(x) order(x, decreasing=TRUE)[1])
|
2939 |
|
-
|
|
2939 |
+
|
2940 |
2940 |
if(CHR == 23){
|
2941 |
2941 |
## the state index for male on chromosome 23 is 2
|
2942 |
2942 |
## add 1 so that the state index is 3 for 'normal' state
|
... |
... |
@@ -2966,11 +2966,11 @@ biasAdjNP <- function(object, cnOptions, tmp.objects){
|
2966 |
2966 |
getParams <- function(object, batch){
|
2967 |
2967 |
##batch <- unique(object$batch)
|
2968 |
2968 |
batch <- unique(batch(object))
|
2969 |
|
- if(length(batch) > 1) stop("batch variable not unique")
|
|
2969 |
+ if(length(batch) > 1) stop("batch variable not unique")
|
2970 |
2970 |
nuA <- as.numeric(fData(object)[, match(paste("nuA", batch, sep="_"), fvarLabels(object))])
|
2971 |
|
- nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])
|
|
2971 |
+ nuB <- as.numeric(fData(object)[, match(paste("nuB", batch, sep="_"), fvarLabels(object))])
|
2972 |
2972 |
phiA <- as.numeric(fData(object)[, match(paste("phiA", batch, sep="_"), fvarLabels(object))])
|
2973 |
|
- phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])
|
|
2973 |
+ phiB <- as.numeric(fData(object)[, match(paste("phiB", batch, sep="_"), fvarLabels(object))])
|
2974 |
2974 |
tau2A <- as.numeric(fData(object)[, match(paste("tau2A", batch, sep="_"), fvarLabels(object))])
|
2975 |
2975 |
tau2B <- as.numeric(fData(object)[, match(paste("tau2B", batch, sep="_"), fvarLabels(object))])
|
2976 |
2976 |
sig2A <- as.numeric(fData(object)[, match(paste("sig2A", batch, sep="_"), fvarLabels(object))])
|
... |
... |
@@ -3012,17 +3012,17 @@ thresholdModelParams <- function(object, cnOptions){
|
3012 |
3012 |
phiB <- getParam(object, "phiB", batch)
|
3013 |
3013 |
if(!all(is.na(phiB))){
|
3014 |
3014 |
phiB[phiB < MIN.PHI] <- MIN.PHI
|
3015 |
|
- object <- pr(object, "phiB", batch, phiB)
|
|
3015 |
+ object <- pr(object, "phiB", batch, phiB)
|
3016 |
3016 |
}
|
3017 |
|
- phiAX <- as.numeric(getParam(object, "phiAX", batch))
|
|
3017 |
+ phiAX <- as.numeric(getParam(object, "phiAX", batch))
|
3018 |
3018 |
if(!all(is.na(phiAX))){
|
3019 |
3019 |
phiAX[phiAX < MIN.PHI] <- MIN.PHI
|
3020 |
|
- object <- pr(object, "phiAX", batch, phiAX)
|
|
3020 |
+ object <- pr(object, "phiAX", batch, phiAX)
|
3021 |
3021 |
}
|
3022 |
3022 |
phiBX <- as.numeric(getParam(object, "phiBX", batch))
|
3023 |
3023 |
if(!all(is.na(phiBX))){
|
3024 |
3024 |
phiBX[phiBX < MIN.PHI] <- MIN.PHI
|
3025 |
|
- object <- pr(object, "phiBX", batch, phiBX)
|
|
3025 |
+ object <- pr(object, "phiBX", batch, phiBX)
|
3026 |
3026 |
}
|
3027 |
3027 |
return(object)
|
3028 |
3028 |
}
|
... |
... |
@@ -3079,11 +3079,11 @@ computeCopynumber.CNSet <- function(object, cnOptions){
|
3079 |
3079 |
verbose <- cnOptions$verbose
|
3080 |
3080 |
##if(verbose) message("Thresholding nu and phi")
|
3081 |
3081 |
object <- thresholdModelParams(object, cnOptions)
|
3082 |
|
- }
|
3083 |
|
- ##if(verbose) message("\nAllele specific copy number")
|
|
3082 |
+ }
|
|
3083 |
+ ##if(verbose) message("\nAllele specific copy number")
|
3084 |
3084 |
object <- polymorphic(object, cnOptions, tmp.objects)
|
3085 |
3085 |
if(any(!isSnp(object))){ ## there are nonpolymorphic probes
|
3086 |
|
- ##if(verbose) message("\nCopy number for nonpolymorphic probes...")
|
|
3086 |
+ ##if(verbose) message("\nCopy number for nonpolymorphic probes...")
|
3087 |
3087 |
object <- nonpolymorphic(object, cnOptions, tmp.objects)
|
3088 |
3088 |
}
|
3089 |
3089 |
##---------------------------------------------------------------------------
|
... |
... |
@@ -3101,9 +3101,9 @@ computeCopynumber.CNSet <- function(object, cnOptions){
|
3101 |
3101 |
object <- pr(object, "phiB", PLATE, getParam(object, "phiB", PLATE))
|
3102 |
3102 |
object <- pr(object, "phiB.se", PLATE, getParam(object, "phiB.se", PLATE))
|
3103 |
3103 |
object <- pr(object, "tau2A", PLATE, getParam(object, "tau2A", PLATE))
|
3104 |
|
- object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))
|
|
3104 |
+ object <- pr(object, "tau2B", PLATE, getParam(object, "tau2B", PLATE))
|
3105 |
3105 |
object <- pr(object, "sig2A", PLATE, getParam(object, "sig2A", PLATE))
|
3106 |
|
- object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))
|
|
3106 |
+ object <- pr(object, "sig2B", PLATE, getParam(object, "sig2B", PLATE))
|
3107 |
3107 |
object <- pr(object, "phiAX", PLATE, as.numeric(getParam(object, "phiAX", PLATE)))
|
3108 |
3108 |
object <- pr(object, "phiBX", PLATE, as.numeric(getParam(object, "phiBX", PLATE)))
|
3109 |
3109 |
object <- pr(object, "corr", PLATE, getParam(object, "corr", PLATE))
|
... |
... |
@@ -3122,3 +3122,50 @@ computeCopynumber.CNSet <- function(object, cnOptions){
|
3122 |
3122 |
|
3123 |
3123 |
|
3124 |
3124 |
|
|
3125 |
+constructFeatureData <- function(gns, cdfName){
|
|
3126 |
+ pkgname <- paste(cdfName, "Crlmm", sep="")
|
|
3127 |
+ path <- system.file("extdata", package=pkgname)
|
|
3128 |
+ load(file.path(path, "cnProbes.rda"))
|
|
3129 |
+ load(file.path(path, "snpProbes.rda"))
|
|
3130 |
+ cnProbes$chr <- chromosome2integer(cnProbes$chr)
|
|
3131 |
+ cnProbes <- as.matrix(cnProbes)
|
|
3132 |
+ snpProbes$chr <- chromosome2integer(snpProbes$chr)
|
|
3133 |
+ snpProbes <- as.matrix(snpProbes)
|
|
3134 |
+ mapping <- rbind(snpProbes, cnProbes, deparse.level=0)
|
|
3135 |
+ mapping <- mapping[match(gns, rownames(mapping)), ]
|
|
3136 |
+ isSnp <- 1L-as.integer(gns %in% rownames(cnProbes))
|
|
3137 |
+ mapping <- cbind(mapping, isSnp, deparse.level=0)
|
|
3138 |
+ stopifnot(identical(rownames(mapping), gns))
|
|
3139 |
+ colnames(mapping) <- c("chromosome", "position", "isSnp")
|
|
3140 |
+ new("AnnotatedDataFrame",
|
|
3141 |
+ data=data.frame(mapping),
|
|
3142 |
+ varMetadata=data.frame(labelDescription=colnames(mapping)))
|
|
3143 |
+}
|
|
3144 |
+constructAssayData <- function(np, snp, object, storage.mode="environment", order.index){
|
|
3145 |
+ stopifnot(identical(snp$gns, featureNames(object)))
|
|
3146 |
+ gns <- c(featureNames(object), np$gns)
|
|
3147 |
+ sns <- np$sns
|
|
3148 |
+ np <- np[1:2]
|
|
3149 |
+ snp <- snp[1:2]
|
|
3150 |
+ stripnames <- function(x) {
|
|
3151 |
+ dimnames(x) <- NULL
|
|
3152 |
+ x
|
|
3153 |
+ }
|
|
3154 |
+ np <- lapply(np, stripnames)
|
|
3155 |
+ snp <- lapply(snp, stripnames)
|
|
3156 |
+ A <- rbind(snp[[1]], np[[1]], deparse.level=0)[order.index, ]
|
|
3157 |
+ B <- rbind(snp[[2]], np[[2]], deparse.level=0)[order.index, ]
|
|
3158 |
+ gt <- stripnames(calls(object))
|
|
3159 |
+ emptyMatrix <- matrix(integer(), nrow(np[[1]]), ncol(A))
|
|
3160 |
+ gt <- rbind(gt, emptyMatrix, deparse.level=0)[order.index,]
|
|
3161 |
+ pr <- stripnames(snpCallProbability(object))
|
|
3162 |
+ pr <- rbind(pr, emptyMatrix, deparse.level=0)[order.index, ]
|
|
3163 |
+ emptyMatrix <- matrix(integer(), nrow(A), ncol(A))
|
|
3164 |
+ aD <- assayDataNew(storage.mode,
|
|
3165 |
+ alleleA=A,
|
|
3166 |
+ alleleB=B,
|
|
3167 |
+ call=gt,
|
|
3168 |
+ callProbability=pr,
|
|
3169 |
+ CA=emptyMatrix,
|
|
3170 |
+ CB=emptyMatrix)
|
|
3171 |
+}
|