Browse code

changes to shrinkGenotypeSummarries and C3 function for estimation of copy number at SNPs on chromosome X

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@52815 bc3139a8-67e5-0310-9ffc-ced21a209358

Rob Scharp authored on 16/02/2011 15:57:05
Showing3 changed files

... ...
@@ -351,7 +351,6 @@ shrinkGenotypeSummaries <- function(strata, index.list, object, MIN.OBS, MIN.SAM
351 351
 	for(k in seq(along=batches)){
352 352
 		B <- batches[[k]]
353 353
 		this.batch <- unique(as.character(batch(object)[B]))
354
-
355 354
 		medianA[[k]] <- cbind(medianA.AA[, k], medianA.AB[, k], medianA.BB[, k])
356 355
 		medianB[[k]] <- cbind(medianB.AA[, k], medianB.AB[, k], medianB.BB[, k])
357 356
 		madA <- cbind(madA.AA[, k], madA.AB[, k], madA.BB[, k])
... ...
@@ -393,8 +392,11 @@ shrinkGenotypeSummaries <- function(strata, index.list, object, MIN.OBS, MIN.SAM
393 392
 		unobservedAB <- NN[, 2] < MIN.OBS
394 393
 		unobservedBB <- NN[, 3] < MIN.OBS
395 394
 		unobserved.index <- vector("list", 3)
395
+		## AB, BB observed
396 396
 		unobserved.index[[1]] <- which(unobservedAA & (NN[, 2] >= MIN.OBS & NN[, 3] >= MIN.OBS))
397
+		## AA and BB observed (strange)
397 398
 		unobserved.index[[2]] <- which(unobservedAB & (NN[, 1] >= MIN.OBS & NN[, 3] >= MIN.OBS))
399
+		## AB and AA observed
398 400
 		unobserved.index[[3]] <- which(unobservedBB & (NN[, 2] >= MIN.OBS & NN[, 1] >= MIN.OBS))
399 401
 		res <- imputeCenter(medianA[[k]], medianB[[k]], index.complete, unobserved.index)
400 402
 		medianA[[k]] <- res[[1]]
... ...
@@ -619,7 +621,8 @@ summarizeMaleXGenotypes <- function(marker.index,
619 621
 				    DF.PRIOR,...){
620 622
 	nr <- length(marker.index)
621 623
 	nc <- length(batchNames(object))
622
-	NN.Mlist <- imputed.medianA <- imputed.medianB <- shrink.madA <- shrink.madB <- vector("list", nc)
624
+##	NN.Mlist <- imputed.medianA <- imputed.medianB <- shrink.madA <- shrink.madB <- vector("list", nc)
625
+	NN.Mlist <- medianA <- medianB <- shrink.madA <- shrink.madB <- vector("list", nc)
623 626
 	gender <- object$gender
624 627
 	GG <- as.matrix(calls(object)[marker.index, gender==1])
625 628
 	CP <- as.matrix(snpCallProbability(object)[marker.index, gender==1])
... ...
@@ -658,43 +661,46 @@ summarizeMaleXGenotypes <- function(marker.index,
658 661
 			tmp
659 662
 		}
660 663
 		statsA.AA <- summaryStats(G.AA, A, FUNS=c("rowMedians", "rowMAD"))
661
-		statsA.AB <- summaryStats(G.AB, A, FUNS=c("rowMedians", "rowMAD"))
664
+		##statsA.AB <- summaryStats(G.AB, A, FUNS=c("rowMedians", "rowMAD"))
662 665
 		statsA.BB <- summaryStats(G.BB, A, FUNS=c("rowMedians", "rowMAD"))
663 666
 		statsB.AA <- summaryStats(G.AA, B, FUNS=c("rowMedians", "rowMAD"))
664
-		statsB.AB <- summaryStats(G.AB, B, FUNS=c("rowMedians", "rowMAD"))
667
+		##statsB.AB <- summaryStats(G.AB, B, FUNS=c("rowMedians", "rowMAD"))
665 668
 		statsB.BB <- summaryStats(G.BB, B, FUNS=c("rowMedians", "rowMAD"))
666
-		medianA <- cbind(statsA.AA[, 1], statsA.AB[, 1], statsA.BB[, 1])
667
-		medianB <- cbind(statsB.AA[, 1], statsB.AB[, 1], statsB.BB[, 1])
668
-		madA <- cbind(statsA.AA[, 1], statsA.AB[, 1], statsA.BB[, 1])
669
-		madB <- cbind(statsB.AA[, 1], statsB.AB[, 1], statsB.BB[, 1])
670
-		rm(statsA.AA, statsA.AB, statsA.BB, statsB.AA, statsB.AB, statsB.BB)
671
-
672
-		NN.M <- cbind(N.AA.M, N.AB.M, N.BB.M)
669
+		medianA[[k]] <- cbind(statsA.AA[, 1], ##statsA.AB[, 1],
670
+				      statsA.BB[, 1])
671
+		medianB[[k]] <- cbind(statsB.AA[, 1], ##statsB.AB[, 1],
672
+				      statsB.BB[, 1])
673
+		madA <- cbind(statsA.AA[, 2],  ##statsA.AB[, 2],
674
+			      statsA.BB[, 2])
675
+		madB <- cbind(statsB.AA[, 2], ##statsB.AB[, 2],
676
+			      statsB.BB[, 2])
677
+		rm(statsA.AA, ##statsA.AB,
678
+		   statsA.BB, statsB.AA,
679
+		   ##statsB.AB,
680
+		   statsB.BB)
681
+		NN.M <- cbind(N.AA.M, ##N.AB.M,
682
+			      N.BB.M)
673 683
 		NN.Mlist[[k]] <- NN.M
674
-
675 684
 		shrink.madA[[k]] <- shrink(madA, NN.M, DF.PRIOR)
676 685
 		shrink.madB[[k]] <- shrink(madB, NN.M, DF.PRIOR)
677
-
678 686
 		##---------------------------------------------------------------------------
679 687
 		## SNPs that we'll use for imputing location/scale of unobserved genotypes
680 688
 		##---------------------------------------------------------------------------
681
-		index.complete <- indexComplete(NN.M[, -2], medianA, medianB, MIN.OBS)
682
-		if(length(index.complete) == 1){
683
-			if(index.complete == FALSE) return()
684
-		}
689
+		##index.complete <- indexComplete(NN.M[, -2], medianA[[k]], medianB[[k]], MIN.OBS)
690
+		index.complete <- indexComplete(NN.M, medianA[[k]], medianB[[k]], MIN.OBS)
685 691
 
686 692
 		##---------------------------------------------------------------------------
687 693
 		## Impute sufficient statistics for unobserved genotypes (plate-specific)
688 694
 		##---------------------------------------------------------------------------
689
-		res <- imputeCenterX(medianA, medianB, NN.M, index.complete, MIN.OBS)
690
-		imputed.medianA[[k]] <- res[[1]]
691
-		imputed.medianB[[k]] <- res[[2]]
695
+		res <- imputeCenterX(medianA[[k]], medianB[[k]], NN.M, index.complete, MIN.OBS)
696
+		medianA[[k]] <- res[[1]]
697
+		medianB[[k]] <- res[[2]]
692 698
 	}
693 699
 	return(list(madA=shrink.madA,
694 700
 		    madB=shrink.madB,
695 701
 		    NN.M=NN.Mlist,
696
-		    medianA=imputed.medianA,
697
-		    medianB=imputed.medianB))
702
+		    medianA=medianA,
703
+		    medianB=medianB))
698 704
 }
699 705
 
700 706
 ## X chromosome, SNPs
... ...
@@ -729,11 +735,14 @@ fit.lm3 <- function(strata,
729 735
 	phiA2 <- as.matrix(phiPrimeA(object)[marker.index, ])
730 736
 	phiB2 <- as.matrix(phiPrimeB(object)[marker.index, ])
731 737
 	if(enough.males){
732
-		res <- summarizeMaleXGenotypes(marker.index=marker.index, batches=batches,
733
-					       object=object, GT.CONF.THR=GT.CONF.THR,
738
+		res <- summarizeMaleXGenotypes(marker.index=marker.index,
739
+					       batches=batches,
740
+					       object=object,
741
+					       GT.CONF.THR=GT.CONF.THR,
734 742
 					       MIN.SAMPLES=MIN.SAMPLES,
735 743
 					       MIN.OBS=MIN.OBS,
736
-					       verbose=verbose, is.lds=is.lds,
744
+					       verbose=verbose,
745
+					       is.lds=is.lds,
737 746
 					       DF.PRIOR=DF.PRIOR/2)
738 747
 		madA.Mlist <- res[["madA"]]
739 748
 		madB.Mlist <- res[["madB"]]
... ...
@@ -803,18 +812,18 @@ fit.lm3 <- function(strata,
803 812
 			phiB2[, k] <- betas[3, ]
804 813
 		}
805 814
 		if(enough.men & !enough.women){
806
-			betas <- fit.wls(NN.M[, c(1,3)],
807
-					 sigma=madA.M[, c(1,3)],
815
+			betas <- fit.wls(NN.M,##[, c(1,3)],
816
+					 sigma=madA.M,##[, c(1,3)],
808 817
 					 allele="A",
809
-					 Y=medianA.M[, c(1,3)],
818
+					 Y=medianA.M,##[, c(1,3)],
810 819
 					 autosome=FALSE,
811
-					 X=cbind(1, c(0, 1)))
820
+					 X=cbind(1, c(1, 0)))
812 821
 			nuA[, k] <- betas[1, ]
813 822
 			phiA[, k] <- betas[2, ]
814
-			betas <- fit.wls(NN.M[, c(1,3)],
815
-					 sigma=madB.M[, c(1,3)],
823
+			betas <- fit.wls(NN.M,##[, c(1,3)],
824
+					 sigma=madB.M,#[, c(1,3)],
816 825
 					 allele="B",
817
-					 Y=medianB.M[, c(1,3)],
826
+					 Y=medianB.M,#[, c(1,3)],
818 827
 					 autosome=FALSE,
819 828
 					 X=cbind(1, c(0, 1)))
820 829
 			nuB[, k] <- betas[1, ]
... ...
@@ -1050,17 +1059,20 @@ imputeCenter <- function(muA, muB, index.complete, index.missing){
1050 1059
 	index <- index.missing
1051 1060
 	mnA <- muA
1052 1061
 	mnB <- muB
1053
-	for(j in 1:3){
1054
-		if(length(index[[j]]) == 0) next()
1055
-		X <- cbind(1, mnA[index.complete,  -j, drop=FALSE], mnB[index.complete,  -j, drop=FALSE])
1056
-		Y <- cbind(mnA[index.complete, j], mnB[index.complete,  j])
1057
-		betahat <- solve(crossprod(X), crossprod(X,Y))
1058
-		X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
1059
-		mus <- X %*% betahat
1060
-		muA[index[[j]], j] <- mus[, 1]
1061
-		muB[index[[j]], j] <- mus[, 2]
1062
+	if(length(index.complete) >= 100){
1063
+		for(j in 1:3){
1064
+			if(length(index[[j]]) == 0) next()
1065
+			X <- cbind(1, mnA[index.complete,  -j, drop=FALSE], mnB[index.complete,  -j, drop=FALSE])
1066
+			Y <- cbind(mnA[index.complete, j], mnB[index.complete,  j])
1067
+			betahat <- solve(crossprod(X), crossprod(X,Y))
1068
+			X <- cbind(1, mnA[index[[j]],  -j, drop=FALSE],  mnB[index[[j]],  -j, drop=FALSE])
1069
+			mus <- X %*% betahat
1070
+			muA[index[[j]], j] <- mus[, 1]
1071
+			muB[index[[j]], j] <- mus[, 2]
1072
+		}
1062 1073
 	}
1063
-	list(muA, muB)
1074
+	results <- list(muA, muB)
1075
+	return(results)
1064 1076
 }
1065 1077
 
1066 1078
 indexComplete <- function(NN, medianA, medianB, MIN.OBS){
... ...
@@ -1069,15 +1081,19 @@ indexComplete <- function(NN, medianA, medianB, MIN.OBS){
1069 1081
 	index.complete <- intersect(Nindex, correct.order)
1070 1082
 	size <- min(5000, length(index.complete))
1071 1083
 	if(size == 5000) index.complete <- sample(index.complete, 5000, replace=TRUE)
1072
-	if(length(index.complete) < 100){
1073
-		warning("fewer than 100 snps pass criteria for imputing unobserve dgenotype location/scale")
1074
-		return(FALSE)
1075
-	}
1084
+##	if(length(index.complete) < 100){
1085
+##		stop("fewer than 100 snps pass criteria for imputing unobserved genotype location/scale")
1086
+##	}
1076 1087
 	return(index.complete)
1077 1088
 }
1078 1089
 
1079 1090
 imputeCentersForMonomorphicSnps <- function(medianA, medianB, index.complete, unobserved.index){
1080 1091
 	cols <- c(3, 1, 2)
1092
+	if(length(index.complete) < 100){
1093
+		##message("index.complete less than 100.  No imputation")
1094
+		results <- list(medianA=medianA, medianB=medianB)
1095
+		return(results)
1096
+	}
1081 1097
 	for(j in 1:3){
1082 1098
 		if(length(unobserved.index[[j]]) == 0) next()
1083 1099
 		kk <- cols[j]
... ...
@@ -1090,36 +1106,40 @@ imputeCentersForMonomorphicSnps <- function(medianA, medianB, index.complete, un
1090 1106
 		medianA[unobserved.index[[j]], -kk] <- mus[, 1:2]
1091 1107
 		medianB[unobserved.index[[j]], -kk] <- mus[, 3:4]
1092 1108
 	}
1093
-	list(medianA=medianA, medianB=medianB)
1109
+	results <- list(medianA=medianA, medianB=medianB)
1110
+	return(results)
1094 1111
 }
1095 1112
 
1096 1113
 
1097 1114
 imputeCenterX <- function(muA, muB, Ns, index.complete, MIN.OBS){
1098
-	index1 <- which(Ns[, 1] == 0 & Ns[, 3] > MIN.OBS)
1115
+	##index1 <- which(Ns[, 1] == 0 & Ns[, 3] > MIN.OBS)
1116
+	index1 <- which(Ns[, 1] == 0 & Ns[, 2] > MIN.OBS)
1099 1117
 	if(length(index1) > 0){
1100
-		X <- cbind(1, muA[index.complete, 3], muB[index.complete, 3])
1118
+		##X <- cbind(1, muA[index.complete, 3], muB[index.complete, 3])
1119
+		X <- cbind(1, muA[index.complete, 2], muB[index.complete, 2])
1101 1120
 		Y <- cbind(1, muA[index.complete, 1], muB[index.complete, 1])
1102 1121
 ##		X <- cbind(1, muA[index.complete[[1]], 3], muB[index.complete[[1]], 3])
1103 1122
 ##		Y <- cbind(1, muA[index.complete[[1]], 1], muB[index.complete[[1]], 1])
1104 1123
 		betahat <- solve(crossprod(X), crossprod(X,Y))
1105 1124
 		##now with the incomplete SNPs
1106
-		X <- cbind(1, muA[index1, 3], muB[index1, 3])
1125
+		##X <- cbind(1, muA[index1, 3], muB[index1, 3])
1126
+		X <- cbind(1, muA[index1, 2], muB[index1, 2])
1107 1127
 		mus <- X %*% betahat
1108 1128
 		muA[index1, 1] <- mus[, 2]
1109 1129
 		muB[index1, 1] <- mus[, 3]
1110 1130
 	}
1111
-	index1 <- which(Ns[, 3] == 0)
1131
+	index1 <- which(Ns[, 2] == 0)
1112 1132
 	if(length(index1) > 0){
1113 1133
 		X <- cbind(1, muA[index.complete, 1], muB[index.complete, 1])
1114
-		Y <- cbind(1, muA[index.complete, 3], muB[index.complete, 3])
1134
+		Y <- cbind(1, muA[index.complete, 2], muB[index.complete, 2])
1115 1135
 ##		X <- cbind(1, muA[index.complete[[2]], 1], muB[index.complete[[2]], 1])
1116 1136
 ##		Y <- cbind(1, muA[index.complete[[2]], 3], muB[index.complete[[2]], 3])
1117 1137
 		betahat <- solve(crossprod(X), crossprod(X,Y))
1118 1138
 		##now with the incomplete SNPs
1119 1139
 		X <- cbind(1, muA[index1, 1], muB[index1, 1])
1120 1140
 		mus <- X %*% betahat
1121
-		muA[index1, 3] <- mus[, 2]
1122
-		muB[index1, 3] <- mus[, 3]
1141
+		muA[index1, 2] <- mus[, 2]
1142
+		muB[index1, 2] <- mus[, 3]
1123 1143
 	}
1124 1144
 	list(muA, muB)
1125 1145
 }
... ...
@@ -1311,11 +1331,12 @@ genotypeSummary <- function(object,
1311 1331
 			    marker.index,
1312 1332
 			    is.lds){
1313 1333
 	if(type == "X.SNP" | type=="X.NP"){
1314
-		gender <- object$gender
1315
-		if(sum(gender == 2) < 3) {
1316
-			message("too few females to estimate within genotype summary statistics on CHR X")
1317
-			return(object)
1318
-		}
1334
+##		gender <- object$gender
1335
+##		## the number of women in each batch could be less than 3...
1336
+##		if(sum(gender == 2) < 3) {
1337
+##			message("too few females to estimate within genotype summary statistics on CHR X")
1338
+##			return(object)
1339
+##		}
1319 1340
 		CHR.X <- TRUE
1320 1341
 	} else CHR.X <- FALSE
1321 1342
 	if(missing(marker.index)){
... ...
@@ -1436,6 +1457,9 @@ summarizeSnps <- function(strata,
1436 1457
 ##	} else {
1437 1458
 ##		batches <- split(seq_along(batch(object)), as.character(batch(object)))
1438 1459
 ##	}
1460
+	if(CHR.X){
1461
+		if(verbose) message("        biallelic cluster medians are estimated using only the women for SNPs on chr. X")
1462
+	}
1439 1463
 	batches <- split(seq_along(batch(object)), as.character(batch(object)))
1440 1464
 	batchnames <- batchNames(object)
1441 1465
 	nr <- length(index)
... ...
@@ -1451,23 +1475,17 @@ summarizeSnps <- function(strata,
1451 1475
 	FL <- as.matrix(flags(object)[index, ])
1452 1476
 	if(verbose) message("        Computing summaries...")
1453 1477
 	for(k in seq_along(batches)){
1454
-		B <- batches[[k]]
1455
-		this.batch <- unique(as.character(batch(object)[B]))
1478
+		##note that the genotype frequency for AA would include 'A' on chromosome X for men
1479
+		sample.index <- batches[[k]]
1480
+		this.batch <- unique(as.character(batch(object)[sample.index]))
1456 1481
 		j <- match(this.batch, batchnames)
1457
-		G <- GG[, B]
1482
+		G <- GG[, sample.index]
1458 1483
 		##NORM <- normal.index[, k]
1459
-		xx <- CP[, B]
1484
+		xx <- CP[, sample.index]
1460 1485
 		highConf <- (1-exp(-xx/1000)) > GT.CONF.THR
1461 1486
 		G <- G*highConf
1462
-		## Some markers may have genotype confidences scores that are ALL below the threshold
1463
-		## For these markers, provide statistical summaries based on all the samples and flag
1464
-		## Provide summaries for these markers and flag to indicate the problem
1465
-		ii <- which(rowSums(G) == 0)
1466
-		G[ii, ] <- GG[ii, B]
1467
-		## table(rowSums(G==0))
1468
-		##G <- G*highConf*NORM
1469
-		A <- AA[, B]
1470
-		B <- BB[, B]
1487
+		A <- AA[, sample.index]
1488
+		B <- BB[, sample.index]
1471 1489
 		## this can be time consuming...do only once
1472 1490
 		G.AA <- G==1
1473 1491
 		G.AA[G.AA==FALSE] <- NA
... ...
@@ -1478,6 +1496,11 @@ summarizeSnps <- function(strata,
1478 1496
 		Ns.AA[, k] <- rowSums(G.AA, na.rm=TRUE)
1479 1497
 		Ns.AB[, k] <- rowSums(G.AB, na.rm=TRUE)
1480 1498
 		Ns.BB[, k] <- rowSums(G.BB, na.rm=TRUE)
1499
+		if(CHR.X){
1500
+			gender <- object$gender[sample.index]
1501
+			sample.index <- sample.index[gender == 2]
1502
+			if(length(sample.index) == 0) next()
1503
+		}
1481 1504
 		summaryStats <- function(X, INT, FUNS){
1482 1505
 			tmp <- matrix(NA, nrow(X), length(FUNS))
1483 1506
 			for(j in seq_along(FUNS)){
... ...
@@ -1486,12 +1509,29 @@ summarizeSnps <- function(strata,
1486 1509
 			}
1487 1510
 			tmp
1488 1511
 		}
1489
-		statsA.AA[[k]] <- summaryStats(G.AA, A, FUNS=c("rowMedians", "rowMAD"))
1490
-		statsA.AB[[k]] <- summaryStats(G.AB, A, FUNS=c("rowMedians", "rowMAD"))
1491
-		statsA.BB[[k]] <- summaryStats(G.BB, A, FUNS=c("rowMedians", "rowMAD"))
1492
-		statsB.AA[[k]] <- summaryStats(G.AA, B, FUNS=c("rowMedians", "rowMAD"))
1493
-		statsB.AB[[k]] <- summaryStats(G.AB, B, FUNS=c("rowMedians", "rowMAD"))
1512
+		stats <- summaryStats(G.AA, A, FUNS=c("rowMedians", "rowMAD"))
1513
+		medianA.AA(object)[index, k] <- stats[, 1]
1514
+		madA.AA(object)[index, k] <- stats[, 2]
1515
+		stats <- summaryStats(G.AB, A, FUNS=c("rowMedians", "rowMAD"))
1516
+		medianA.AB(object)[index, k] <- stats[, 1]
1517
+		madA.AB(object)[index, k] <- stats[, 2]
1518
+
1519
+		stats <- summaryStats(G.BB, A, FUNS=c("rowMedians", "rowMAD"))
1520
+		medianA.BB(object)[index, k] <- stats[, 1]
1521
+		madA.BB(object)[index, k] <- stats[, 2]
1522
+
1523
+		stats <- summaryStats(G.AA, B, FUNS=c("rowMedians", "rowMAD"))
1524
+		medianB.AA(object)[index, k] <- stats[, 1]
1525
+		madB.AA(object)[index, k] <- stats[, 2]
1526
+
1527
+		stats <- summaryStats(G.AB, B, FUNS=c("rowMedians", "rowMAD"))
1528
+		medianB.AB(object)[index, k] <- stats[, 1]
1529
+		madB.AB(object)[index, k] <- stats[, 2]
1530
+
1494 1531
 		statsB.BB[[k]] <- summaryStats(G.BB, B, FUNS=c("rowMedians", "rowMAD"))
1532
+		medianB.BB(object)[index, k] <- stats[, 1]
1533
+		madB.BB(object)[index, k] <- stats[, 2]
1534
+
1495 1535
 		## log2 Transform Intensities
1496 1536
 		A <- log2(A); B <- log2(B)
1497 1537
 		tau2A.AA[, k] <- summaryStats(G.AA, A, FUNS="rowMAD")^2
... ...
@@ -1510,19 +1550,19 @@ summarizeSnps <- function(strata,
1510 1550
 	corrAA(object)[index,] <- corrAA
1511 1551
 	corrAB(object)[index,] <- corrAB
1512 1552
 	corrBB(object)[index,] <- corrBB
1513
-	medianA.AA(object)[index,] <- do.call(cbind, lapply(statsA.AA, function(x) x[, 1]))
1514
-	medianA.AB(object)[index,] <- do.call(cbind, lapply(statsA.AB, function(x) x[, 1]))
1515
-	medianA.BB(object)[index,] <- do.call(cbind, lapply(statsA.BB, function(x) x[, 1]))
1516
-	medianB.AA(object)[index,] <- do.call(cbind, lapply(statsB.AA, function(x) x[, 1]))
1517
-	medianB.AB(object)[index,] <- do.call(cbind, lapply(statsB.AB, function(x) x[, 1]))
1518
-	medianB.BB(object)[index,] <- do.call(cbind, lapply(statsB.BB, function(x) x[, 1]))
1519
-
1520
-	madA.AA(object)[index,] <- do.call(cbind, lapply(statsA.AA, function(x) x[, 2]))
1521
-	madA.AB(object)[index,] <- do.call(cbind, lapply(statsA.AB, function(x) x[, 2]))
1522
-	madA.BB(object)[index,] <- do.call(cbind, lapply(statsA.BB, function(x) x[, 2]))
1523
-	madB.AA(object)[index,] <- do.call(cbind, lapply(statsB.AA, function(x) x[, 2]))
1524
-	madB.AB(object)[index,] <- do.call(cbind, lapply(statsB.AB, function(x) x[, 2]))
1525
-	madB.BB(object)[index,] <- do.call(cbind, lapply(statsB.BB, function(x) x[, 2]))
1553
+##	medianA.AA(object)[index,] <- do.call(cbind, lapply(statsA.AA, function(x) x[, 1]))
1554
+##	medianA.AB(object)[index,] <- do.call(cbind, lapply(statsA.AB, function(x) x[, 1]))
1555
+##	medianA.BB(object)[index,] <- do.call(cbind, lapply(statsA.BB, function(x) x[, 1]))
1556
+##	medianB.AA(object)[index,] <- do.call(cbind, lapply(statsB.AA, function(x) x[, 1]))
1557
+##	medianB.AB(object)[index,] <- do.call(cbind, lapply(statsB.AB, function(x) x[, 1]))
1558
+##	medianB.BB(object)[index,] <- do.call(cbind, lapply(statsB.BB, function(x) x[, 1]))
1559
+##
1560
+##	madA.AA(object)[index,] <- do.call(cbind, lapply(statsA.AA, function(x) x[, 2]))
1561
+##	madA.AB(object)[index,] <- do.call(cbind, lapply(statsA.AB, function(x) x[, 2]))
1562
+##	madA.BB(object)[index,] <- do.call(cbind, lapply(statsA.BB, function(x) x[, 2]))
1563
+##	madB.AA(object)[index,] <- do.call(cbind, lapply(statsB.AA, function(x) x[, 2]))
1564
+##	madB.AB(object)[index,] <- do.call(cbind, lapply(statsB.AB, function(x) x[, 2]))
1565
+##	madB.BB(object)[index,] <- do.call(cbind, lapply(statsB.BB, function(x) x[, 2]))
1526 1566
 	tau2A.AA(object)[index, ] <- tau2A.AA
1527 1567
 ##	tau2A.AB(object)[index, ] <- tau2A.AB
1528 1568
 	tau2A.BB(object)[index, ] <- tau2A.BB
... ...
@@ -356,6 +356,7 @@ C3 <- function(object, allele, marker.index, batch.index, sample.index){
356 356
 		l <- batch.index[k]
357 357
 		##j <- which(as.character(batch(object))[sample.index] == batchNames(object)[l])
358 358
 		jj <- sample.index[as.character(batch(object))[sample.index] == batchNames(object)[l]]
359
+		##phiA2 and phiB2 are not always estimable  -- need both men and women
359 360
 		phiA2 <- phiPrimeA(object)[marker.index, l]
360 361
 		phiB2 <- phiPrimeB(object)[marker.index, l]
361 362
 		phiA <- phiA(object)[marker.index, l]
... ...
@@ -364,28 +365,28 @@ C3 <- function(object, allele, marker.index, batch.index, sample.index){
364 365
 		nuB <- nuB(object)[marker.index, l]
365 366
 		IA <- A(object)[marker.index, jj]
366 367
 		IB <- B(object)[marker.index, jj]
368
+		## I = nu + acn * phi
369
+		## acn = 1/phi* (I-nu)
367 370
 		phistar <- phiB2/phiA
368 371
 		tmp <- (IB - nuB - phistar*IA + phistar*nuA)/phiB
369 372
 		CB <- tmp/(1-phistar*phiA2/phiB)
370
-		##CB <- 1/(1-phiA2*phiB2/(phiA*phiB)) * 1/phiB * (IA-nuB-phiB2/phiA*(IA-nuA))
371
-		CA <- (IA-nuA-phiA2*CB)/phiA
373
+		index <- which(is.na(phiB2) | is.na(phiA2))
374
+		if(length(index) > 0){
375
+			cb <- 1/phiB[index] * (IB[index, ] - nuB[index])
376
+			CB[index, ] <- cb
377
+		}
372 378
 		if(allele == "B"){
373 379
 			acn[, match(jj, sample.index)] <- CB
374 380
 			##acn[[k]] <- CB
375 381
 		}
376 382
 		if(allele == "A"){
377
-			acn[, match(jj, sample.index)] <- (IA-nuA-phiA2*CB)/phiA
378
-		}
379
-		if(allele == "AandB"){
380
-			CA <- tmp/(1-phistar*phiA2/phiB)
381
-			CB <- (IA-nuA-phiA2*CB)/phiA
382
-			acn[, match(jj, sample.index)] <- (IA-nuA-phiA2*CB)/phiA
383
+			CA <- (IA-nuA-phiA2*CB)/phiA
384
+			if(length(index) > 0){
385
+				ca <- 1/phiA[index] * (IA[index, ] - nuA[index])
386
+ 				CA[index, ] <- ca
387
+			}
388
+			acn[, match(jj, sample.index)] <- CA
383 389
 		}
384
-##		if(allele=="AandB")
385
-##			CA <- tmp/(1-phistar*phiA2/phiB)
386
-##			CB <- (IA-nuA-phiA2*CB)/phiA
387
-##			acn[[k]] <- CA+CB
388
-##		}
389 390
 	}
390 391
 ##	if(length(acn) > 1){
391 392
 ##		acn <- do.call("cbind", acn)
... ...
@@ -1,3 +1,4 @@
1
-rsync -avuzb --exclude '~*' --exclude '.git*' -e ssh ~/madman/Rpacks/git/crlmm/ enigma2.jhsph.edu:~/madman/Rpacks/release/crlmm
2
-rsync -avuzb --exclude '.git*' -e ssh enigma2.jhsph.edu:~/madman/Rpacks/release/crlmm/inst/scripts/copynumber.pdf ~/madman/Rpacks/git/crlmm/inst/scripts/copynumber.pdf
1
+rsync -avuzb --exclude '*~' --exclude '.git*' -e ssh ~/madman/Rpacks/git/crlmm/ enigma2.jhsph.edu:~/madman/Rpacks/git/crlmm
2
+rsync -avuzb --exclude '*~' --exclude '.git*' -e ssh enigma2.jhsph.edu:~/madman/Rpacks/release/crlmm/inst/scripts/*opynumber.pdf ~/madman/Rpacks/git/crlmm/inst/scripts/
3
+
3 4