Browse code

Commit made by the Bioconductor Git-SVN bridge. Consists of 1 commit.

Commit information:

Commit id: 736588e3d47fd58358f61ce2241bb2c85b4d9a0c

update to vignette to process previously unevaluated codechunks

Committed by: Rob Scharpf
Author Name: Rob Scharpf
Commit date: 2014-09-19 10:14:04 -0400
Author date: 2014-09-19 10:14:04 -0400


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

Rob Scharp authored on 19/09/2014 14:22:04
Showing1 changed files
... ...
@@ -681,7 +681,7 @@ setMethod(OligoSetList, "CNSet", function(object,...){
681 681
 	constructOligoSetListFrom(object, ...)
682 682
 })
683 683
 setMethod(BafLrrSetList, "CNSet", function(object,...){
684
-	constructBafLrrSetListFrom(object, ...)
684
+  constructBafLrrSetListFrom(object, ...)
685 685
 })
686 686
 
687 687
 
... ...
@@ -689,41 +689,41 @@ setMethod(BafLrrSetList, "CNSet", function(object,...){
689 689
 
690 690
 
691 691
 constructOligoSetListFrom <- function(object, ...){
692
-	##row.index <- seq_len(nrow(object))
693
-	##col.index <- seq_len(ncol(object))
694
-	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
695
-	if(is.lds) stopifnot(isPackageLoaded("ff"))
696
-	b.r <- calculateRBaf(object, ...)
697
-	b <- b.r[["baf"]]
698
-	r <- b.r[["lrr"]]
699
-	j <- match(colnames(r[[1]]), sampleNames(object))
700
-	rns <- lapply(r, rownames)
701
-	fDList <- foreach(featureid=rns) %do%{
702
-		featureData(object)[match(featureid, featureNames(object)), ]
703
-	}
704
-	names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
705
-	gtPlist <- gtlist <- vector("list", length(r))
706
-	for(i in seq_along(r)){
707
-		gtlist[[i]] <- initializeBigMatrix("call", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
708
-		gtPlist[[i]] <- initializeBigMatrix("callPr", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
709
-		featureid <- rownames(r[[i]])
710
-		ix <- match(featureid, featureNames(object))
711
-		rownames(gtPlist[[i]]) <- rownames(gtlist[[i]]) <- featureid
712
-		colnames(gtPlist[[i]]) <- colnames(gtlist[[i]]) <- colnames(r[[i]])
713
-		for(k in seq_along(j)){
714
-			gtlist[[i]][, k] <- calls(object)[ix, j[k]]
715
-			gtPlist[[i]][, k] <- snpCallProbability(object)[ix, j[k]]
716
-		}
717
-	}
718
-	ad <- AssayDataList(baf=b, copyNumber=r, call=gtlist, callProbability=gtPlist)
719
-	object <- new("oligoSetList",
720
-		      assayDataList=ad,
721
-		      featureDataList=fDList,
722
-		      chromosome=names(fDList),
723
-		      phenoData=phenoData(object)[j, ],
724
-		      annotation=annotation(object),
725
-		      genome=genomeBuild(object))
726
-	return(object)
692
+  ##row.index <- seq_len(nrow(object))
693
+  ##col.index <- seq_len(ncol(object))
694
+  is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
695
+  if(is.lds) stopifnot(isPackageLoaded("ff"))
696
+  b.r <- calculateRBaf(object, ...)
697
+  b <- b.r[["baf"]]
698
+  r <- b.r[["lrr"]]
699
+  j <- match(colnames(r[[1]]), sampleNames(object))
700
+  rns <- lapply(r, rownames)
701
+  fDList <- foreach(featureid=rns) %do%{
702
+    featureData(object)[match(featureid, featureNames(object)), ]
703
+  }
704
+  names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
705
+  gtPlist <- gtlist <- vector("list", length(r))
706
+  for(i in seq_along(r)){
707
+    gtlist[[i]] <- initializeBigMatrix("call", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
708
+    gtPlist[[i]] <- initializeBigMatrix("callPr", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
709
+    featureid <- rownames(r[[i]])
710
+    ix <- match(featureid, featureNames(object))
711
+    rownames(gtPlist[[i]]) <- rownames(gtlist[[i]]) <- featureid
712
+    colnames(gtPlist[[i]]) <- colnames(gtlist[[i]]) <- colnames(r[[i]])
713
+    for(k in seq_along(j)){
714
+      gtlist[[i]][, k] <- calls(object)[ix, j[k]]
715
+      gtPlist[[i]][, k] <- snpCallProbability(object)[ix, j[k]]
716
+    }
717
+  }
718
+  ad <- AssayDataList(baf=b, copyNumber=r, call=gtlist, callProbability=gtPlist)
719
+  object <- new("oligoSetList",
720
+                assayDataList=ad,
721
+                featureDataList=fDList,
722
+                chromosome=names(fDList),
723
+                phenoData=phenoData(object)[j, ],
724
+                annotation=annotation(object),
725
+                genome=genomeBuild(object))
726
+  return(object)
727 727
 }
728 728
 
729 729
 
Browse code

Removed crlmm::: from .R files and set saveDate=FALSE in genotype.Illumina() man page to clean up Notes and Warnings

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

unknown authored on 28/03/2014 12:43:01
Showing1 changed files
... ...
@@ -571,7 +571,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){
571 571
 		for(j in seq_along(sampleindex)){
572 572
 			bname <- batch.name[j]
573 573
 			J <- sampleindex[[j]]
574
-			res <- crlmm:::calculateRTheta(object=object,
574
+			res <- calculateRTheta(object=object, # crlmm:::
575 575
 						       batch.name=bname,
576 576
 						       feature.index=i)
577 577
 			k <- match(sampleNames(object)[J], sns)
Browse code

merge local changes with dan's changes on bioc

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

Rob Scharp authored on 20/01/2013 18:31:47
Showing1 changed files
... ...
@@ -1,5 +1,6 @@
1 1
 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))
2 2
 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))
3
+##setReplaceMethod("mixtureParams", signature(object="CNSet", value="ANY"), function(object, value) object@mixtureParams <- mixtureParams)
3 4
 
4 5
 linearParamElementReplace <- function(obj, elt, value) {
5 6
     storage.mode <- storageMode(batchStatistics(obj))
... ...
@@ -668,3 +669,87 @@ calculateRBafCNSet <- function(object, batch.name, chrom){
668 669
 ##			      lrr=lrr)
669 670
 ##		  return(res)
670 671
 ##	  })
672
+
673
+
674
+setAs("CNSet", "oligoSetList", function(from, to){
675
+	constructOligoSetListFrom(from)
676
+})
677
+
678
+
679
+
680
+setMethod(OligoSetList, "CNSet", function(object,...){
681
+	constructOligoSetListFrom(object, ...)
682
+})
683
+setMethod(BafLrrSetList, "CNSet", function(object,...){
684
+	constructBafLrrSetListFrom(object, ...)
685
+})
686
+
687
+
688
+
689
+
690
+
691
+constructOligoSetListFrom <- function(object, ...){
692
+	##row.index <- seq_len(nrow(object))
693
+	##col.index <- seq_len(ncol(object))
694
+	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
695
+	if(is.lds) stopifnot(isPackageLoaded("ff"))
696
+	b.r <- calculateRBaf(object, ...)
697
+	b <- b.r[["baf"]]
698
+	r <- b.r[["lrr"]]
699
+	j <- match(colnames(r[[1]]), sampleNames(object))
700
+	rns <- lapply(r, rownames)
701
+	fDList <- foreach(featureid=rns) %do%{
702
+		featureData(object)[match(featureid, featureNames(object)), ]
703
+	}
704
+	names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
705
+	gtPlist <- gtlist <- vector("list", length(r))
706
+	for(i in seq_along(r)){
707
+		gtlist[[i]] <- initializeBigMatrix("call", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
708
+		gtPlist[[i]] <- initializeBigMatrix("callPr", nr=nrow(r[[i]]), nc=length(j), vmode="integer")
709
+		featureid <- rownames(r[[i]])
710
+		ix <- match(featureid, featureNames(object))
711
+		rownames(gtPlist[[i]]) <- rownames(gtlist[[i]]) <- featureid
712
+		colnames(gtPlist[[i]]) <- colnames(gtlist[[i]]) <- colnames(r[[i]])
713
+		for(k in seq_along(j)){
714
+			gtlist[[i]][, k] <- calls(object)[ix, j[k]]
715
+			gtPlist[[i]][, k] <- snpCallProbability(object)[ix, j[k]]
716
+		}
717
+	}
718
+	ad <- AssayDataList(baf=b, copyNumber=r, call=gtlist, callProbability=gtPlist)
719
+	object <- new("oligoSetList",
720
+		      assayDataList=ad,
721
+		      featureDataList=fDList,
722
+		      chromosome=names(fDList),
723
+		      phenoData=phenoData(object)[j, ],
724
+		      annotation=annotation(object),
725
+		      genome=genomeBuild(object))
726
+	return(object)
727
+}
728
+
729
+
730
+
731
+constructBafLrrSetListFrom <- function(object, ...){
732
+	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
733
+	if(is.lds) stopifnot(isPackageLoaded("ff"))
734
+	b.r <- calculateRBaf(object, ...)
735
+	b <- b.r[["baf"]]
736
+	r <- b.r[["lrr"]]
737
+	j <- match(colnames(r[[1]]), sampleNames(object))
738
+	rns <- lapply(r, rownames)
739
+	featureid <- NULL
740
+	fDList <- foreach(featureid=rns) %do%{
741
+		featureData(object)[match(featureid, featureNames(object)), ]
742
+	}
743
+	names(fDList) <- sapply(fDList, function(x) chromosome(x)[1])
744
+	ad <- AssayDataList(baf=b, lrr=r)
745
+	object <- new("BafLrrSetList",
746
+		      assayDataList=ad,
747
+		      featureDataList=fDList,
748
+		      chromosome=names(fDList),
749
+		      phenoData=phenoData(object)[j, ],
750
+		      annotation=annotation(object),
751
+		      genome=genomeBuild(object))
752
+	return(object)
753
+}
754
+
755
+computeBR <- constructBafLrrSetListFrom
Browse code

merge collab branch

* collab:
update vignettes/Makefile
update .gitignore
update table in CopyNumberViewiew
bump dependency on oligoClasses version
Update alias for "[" method for PredictionRegion objects
Replaced CopyNumberOverview and Infrastructure vignettes with those in the release branch. The versions in release have removed AffymetrixPreprocessCN, which is no longer available

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

Rob Scharp authored on 19/11/2012 04:49:34
Showing1 changed files
... ...
@@ -583,6 +583,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){
583 583
 	if(isPackageLoaded("ff")){
584 584
 		pkgs <- c("oligoClasses", "ff", "Biobase", "crlmm")
585 585
 	} else pkgs <- c("oligoClasses", "Biobase", "crlmm")
586
+	i <- NULL
586 587
 	res <- foreach(i=indexlist, chr=names(indexlist), .packages=pkgs) %dopar% {
587 588
 		processByChromosome(i, chr, path)
588 589
 	}
Browse code

merge with collab branch containing fix for dqrls and bug-fix for computeRBaf that can misalign sample index with batch index (when batch is not in alphabetical order)

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

Rob Scharp authored on 17/11/2012 15:24:46
Showing1 changed files
... ...
@@ -548,7 +548,7 @@ calculateRBafCNSet <- function(object, batch.name, chrom){
548 548
 	indexlist <- split(index, chr[index])
549 549
 	J <- which(batch(object) %in% batch.name)
550 550
 	sns <- sampleNames(object)[J]
551
-	sampleindex <- split(J, batch(object)[J])
551
+	sampleindex <- split(J, factor(batch(object)[J], levels=unique(batch(object)[J])))
552 552
 	if(!all(valid.chrs)) warning("Only computing log R ratios and BAFs for autosomes and chr X")
553 553
 	## if ff package is loaded, these will be ff objects
554 554
 	chr <- names(indexlist)
Browse code

Merge branch 'collab'

* collab:
update Rds
revert imputeGender to original api
bug fix for calculateRBafCNSet
update calculateRBafCNSet
Fix documentation for crlmmCopynumber -- added argument fitLinearModel
update crlmmGT2
Put rm(DD, ...) further down in crlmmGT2 function
remove message about cloning A and B
Open and close callsPr and callsGt in crlmmGT2 (when args not missing)
revert removed indices loaded in crlmmGT2
snprmaAffy no longer writes normalized intensities to calls and callProbability slots. crlmmGT2 takes arguments callsGt and callsPr. When present, crlmmGT2 will not overwrite A and B.
explicit coercion to matrix in imputeGender
Assign imputed gender to cnSet$gender within crlmmGT2 function.
Change sum(SNR > SNRmin) to sum(SNR[] > SNRmin) in imputeGender

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

Rob Scharp authored on 14/09/2012 20:07:37
Showing1 changed files
... ...
@@ -535,7 +535,11 @@ setMethod("calculateRBaf", signature(object="CNSet"),
535 535
 	  })
536 536
 
537 537
 calculateRBafCNSet <- function(object, batch.name, chrom){
538
-	if(missing(batch.name)) batch.name <- batchNames(object)
538
+	if(missing(batch.name)) {
539
+		batch.name <- batchNames(object)
540
+		if("grandMean" %in% batch.name)
541
+			batch.name <- batch.name[-length(batch.name)]
542
+	}
539 543
 	if(missing(chrom)) chrom <- unique(chromosome(object))
540 544
 	if(!(all(batch.name %in% batchNames(object)))) stop("batch.name must be belong to batchNames(object)")
541 545
 	chr <- chromosome(object)
... ...
@@ -549,35 +553,42 @@ calculateRBafCNSet <- function(object, batch.name, chrom){
549 553
 	## if ff package is loaded, these will be ff objects
550 554
 	chr <- names(indexlist)
551 555
 	rlist <- blist <- vector("list", length(indexlist))
552
-	for(i in seq_along(indexlist)){
553
-		I <- indexlist[[i]]
554
-		nr <- length(I)
555
-		CHR <- names(indexlist)[i]
556
-		bafname <- paste("baf_chr", CHR, sep="")
557
-		rname <- paste("lrr_chr", CHR, sep="")
556
+	path <- ldPath()
557
+	processByChromosome <- function(i, chr, path){
558
+		ldPath(path)
559
+		nr <- length(i)
560
+		##CHR <- names(i)
561
+		bafname <- paste("baf_chr", chr, sep="")
562
+		rname <- paste("lrr_chr", chr, sep="")
558 563
 		bmatrix <- initializeBigMatrix(bafname, nr=nr, nc=length(sns), vmode="integer")
559 564
 		rmatrix <- initializeBigMatrix(rname, nr=nr, nc=length(sns), vmode="integer")
560 565
 		colnames(rmatrix) <- colnames(bmatrix) <- sns
561 566
 		## put rownames in order of physical position
562
-		ix <- order(position(object)[I])
563
-		I <- I[ix]
564
-		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[I]
567
+		ix <- order(position(object)[i])
568
+		i <- i[ix]
569
+		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[i]
565 570
 		for(j in seq_along(sampleindex)){
566 571
 			bname <- batch.name[j]
567 572
 			J <- sampleindex[[j]]
568
-			res <- calculateRTheta(object=object,
569
-					       batch.name=bname,
570
-					       feature.index=I)
573
+			res <- crlmm:::calculateRTheta(object=object,
574
+						       batch.name=bname,
575
+						       feature.index=i)
571 576
 			k <- match(sampleNames(object)[J], sns)
572 577
 			bmatrix[, k] <- res[["baf"]]
573 578
 			rmatrix[, k] <- res[["lrr"]]
574
-			##colnames(bmatrix)[J] <- colnames(rmatrix)[J] <- sampleNames(object)[J]
575 579
 		}
576
-		blist[[i]] <- bmatrix
577
-		rlist[[i]] <- rmatrix
580
+		list(bmatrix, rmatrix)
581
+	}
582
+	## calcualte R BAF by chromosome
583
+	if(isPackageLoaded("ff")){
584
+		pkgs <- c("oligoClasses", "ff", "Biobase", "crlmm")
585
+	} else pkgs <- c("oligoClasses", "Biobase", "crlmm")
586
+	res <- foreach(i=indexlist, chr=names(indexlist), .packages=pkgs) %dopar% {
587
+		processByChromosome(i, chr, path)
578 588
 	}
579
-	res <- list(baf=blist,
580
-		    lrr=rlist)
589
+	blist <- lapply(res, "[[", 1)
590
+	rlist <- lapply(res, "[[", 2)
591
+	res <- list(baf=blist, lrr=rlist)
581 592
 	return(res)
582 593
 }
583 594
 
Browse code

Merge branch 'collab'

* collab: (34 commits)
revert change to IlluminaPreprocessCN
fix bug in isValidCdfName
print warning when all features in a batch of probes are flagged, but allow processing to continue
add utility cleancdfnames
Add validCdfNames.Rd
export validCdfNames
imputeGender fix when chromosome Y not available
Use splitIndicesByLength(index, ocSamples/getDoParWorkers())
Can not allocate vector of size XG with genotype.Illumina. Use splitIndicesByNode() only if the length of the list is greater than the split from splitIndicesByLength(). Otherwise, split by length using ocSamples()
update .gitignore
Add make.unique for sampleSheet$Sample_ID in readIdatFiles
bug in description
ensure sample ids stored in samplesheet are unique when constructing cnSet object
update oligoClasses dependency
update unit test for genotype.Illumina
revert change in constructInf call from genotype.Illumina
Update genotype.Rd
edit ACN function
1.15.6 use make.unique(basename(arrayNames)) to allow processing of Illumina samples with duplicated barcodes
check that sample identifies are unique in crlmm function
...

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

Rob Scharp authored on 08/07/2012 19:00:03
Showing1 changed files
... ...
@@ -1,63 +1,6 @@
1 1
 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))
2 2
 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))
3 3
 
4
-setAs("CNSet", "oligoSnpSet", function(from, to){
5
-	cnSet2oligoSnpSet(from)
6
-})
7
-
8
-cnSet2oligoSnpSet <- function(object){
9
-	row.index <- seq_len(nrow(object))
10
-	col.index <- seq_len(ncol(object))
11
-	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
12
-	if(is.lds) stopifnot(isPackageLoaded("ff"))
13
-	b.r <- calculateRBaf(object)
14
-	b <- integerMatrix(b.r[[1]], 1000)
15
-	r <- integerMatrix(b.r[[2]], 100)
16
-##	if(is.lds){
17
-##		## initialize a big matrix for raw copy number
18
-##		message("creating an ff object for storing total copy number")
19
-##		tcn <- initializeBigMatrix(name="total_cn", nrow(object), ncol(object), vmode="double")
20
-##		for(j in 1:ncol(object)){
21
-##			tcn[, j] <- totalCopynumber(object, i=row.index, j=j)
22
-##		}
23
-##	} else {
24
-##		if(ncol(object) > 5){
25
-##			##this can be memory intensive, so we try to be careful
26
-##			col.index <- splitIndicesByLength(seq(length=ncol(object)), 5)
27
-##			tcn <- matrix(NA, nrow(object), ncol(object))
28
-##			dimnames(tcn) <- list(featureNames(object), sampleNames(object))
29
-##			rows <- 1:nrow(object)
30
-##			for(i in seq_along(col.index)){
31
-##				cat(".")
32
-##				j <- col.index[[i]]
33
-##				cnSet <- object[, j]
34
-##				tcn[, j] <- totalCopynumber(cnSet, i=row.index, j=1:ncol(cnSet))
35
-##				rm(cnSet); gc()
36
-##			}
37
-##			cat("\n")
38
-##		} else {
39
-##			tcn <- totalCopynumber(object, i=row.index, j=col.index)
40
-##		}
41
-##	}
42
-##	message("Transforming copy number to log2 scale")
43
-##	tcn[tcn < 0.1] <- 0.1
44
-##	tcn[tcn > 8] <- 8
45
-##	log.tcn <- log2(tcn)
46
-	tmp <- new("oligoSnpSet",
47
-		   copyNumber=integerMatrix(b.r[[2]], scale=100),
48
-		   call=calls(object),
49
-		   callProbability=snpCallProbability(object),
50
-		   annotation=annotation(object),
51
-		   featureData=featureData(object),
52
-		   phenoData=phenoData(object),
53
-		   experimentData=experimentData(object),
54
-		   protocolData=protocolData(object),
55
-		   is.integer=FALSE)
56
-	tmp <- assayDataElementReplace(tmp, "baf", integerMatrix(b.r[[1]], scale=1000))
57
-	return(tmp)
58
-}
59
-
60
-
61 4
 linearParamElementReplace <- function(obj, elt, value) {
62 5
     storage.mode <- storageMode(batchStatistics(obj))
63 6
     switch(storage.mode,
... ...
@@ -185,7 +128,7 @@ ACN <- function(object, allele, i , j){
185 128
 	}
186 129
 	## Define batch.index and sample.index
187 130
 	if(!missing.j) {
188
-		batches <- unique(as.character(batch(object))[j])
131
+		batches <- unique(as.character(batch(object)[j]))
189 132
 		##batches <- as.character(batch(object)[j])
190 133
 		batch.index <- match(batches, batchNames(object))
191 134
 	} else {
... ...
@@ -587,78 +530,129 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),
587 530
 })
588 531
 
589 532
 setMethod("calculateRBaf", signature(object="CNSet"),
590
-	  function(object, batch.name){
591
-		  all.autosomes <- all(chromosome(object) < 23)
592
-		  if(!all.autosomes){
593
-			  stop("method currently only defined for chromosomes 1-22")
594
-		  }
595
-		  if(missing(batch.name)) batch.name <- batchNames(object)[1]
596
-		  stopifnot(batch.name %in% batchNames(object))
597
-		  if(length(batch.name) > 1){
598
-			  warning("only the first batch in batch.name processed")
599
-			  batch.name <- batch.name[1]
600
-		  }
601
-		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)
602
-		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)
603
-		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)
604
-
605
-		  J <- which(batch(object) == batch.name)
606
-
607
-		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
608
-		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
609
-		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
610
-
611
-		  a <- A(object)[, J, drop=FALSE]
612
-		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
613
-		  is.np <- !isSnp(object)
614
-		  ##b[is.np, ] <- a[is.np, ]
615
-		  b[is.np, ] <- 0L
616
-		  dns <- dimnames(a)
617
-		  dimnames(a) <- dimnames(b) <- NULL
618
-		  obs.theta <- atan2(b, a)*2/pi
619
-
620
-		  lessAA <- obs.theta < theta.aa
621
-		  lessAB <- obs.theta < theta.ab
622
-		  lessBB <- obs.theta < theta.bb
623
-		  grAA <- !lessAA
624
-		  grAB <- !lessAB
625
-		  grBB <- !lessBB
626
-		  not.na <- !is.na(theta.aa)
627
-		  I1 <- grAA & lessAB & not.na
628
-		  I2 <- grAB & lessBB & not.na
629
-
630
-		  bf <- matrix(NA, nrow(object), ncol(a))
631
-		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
632
-		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
633
-		  bf[lessAA] <- 0
634
-		  bf[grBB] <- 1
635
-
636
-		  r.expected <- matrix(NA, nrow(object), ncol(a))
637
-		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), length(J), byrow=FALSE)
638
-		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), length(J), byrow=FALSE)
639
-		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), length(J), byrow=FALSE)
640
-		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
641
-		  obs.r <- a+b
642
-
643
-		  lessAA <- lessAA & not.na
644
-		  grBB <- grBB & not.na
645
-		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
646
-		  r.expected[I1] <- tmp
647
-		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
648
-		  r.expected[I2] <- tmp
649
-		  r.expected[lessAA] <- r.aa[lessAA]
650
-		  r.expected[grBB] <- r.bb[grBB]
533
+	  function(object, batch.name, chrom){
534
+		  calculateRBafCNSet(object, batch.name, chrom)
535
+	  })
651 536
 
652
-		  index.np <- which(is.np)
653
-		  if(length(index.np) > 0){
654
-			  a.np <- A(object)[index.np, J]
655
-			  meds <- rowMedians(a.np, na.rm=TRUE)
656
-			  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
657
-		  }
658
-		  lrr <- log2(obs.r/r.expected)
537
+calculateRBafCNSet <- function(object, batch.name, chrom){
538
+	if(missing(batch.name)) batch.name <- batchNames(object)
539
+	if(missing(chrom)) chrom <- unique(chromosome(object))
540
+	if(!(all(batch.name %in% batchNames(object)))) stop("batch.name must be belong to batchNames(object)")
541
+	chr <- chromosome(object)
542
+	valid.chrs <- chr <= 23 & chr %in% chrom
543
+	index <- which(valid.chrs)
544
+	indexlist <- split(index, chr[index])
545
+	J <- which(batch(object) %in% batch.name)
546
+	sns <- sampleNames(object)[J]
547
+	sampleindex <- split(J, batch(object)[J])
548
+	if(!all(valid.chrs)) warning("Only computing log R ratios and BAFs for autosomes and chr X")
549
+	## if ff package is loaded, these will be ff objects
550
+	chr <- names(indexlist)
551
+	rlist <- blist <- vector("list", length(indexlist))
552
+	for(i in seq_along(indexlist)){
553
+		I <- indexlist[[i]]
554
+		nr <- length(I)
555
+		CHR <- names(indexlist)[i]
556
+		bafname <- paste("baf_chr", CHR, sep="")
557
+		rname <- paste("lrr_chr", CHR, sep="")
558
+		bmatrix <- initializeBigMatrix(bafname, nr=nr, nc=length(sns), vmode="integer")
559
+		rmatrix <- initializeBigMatrix(rname, nr=nr, nc=length(sns), vmode="integer")
560
+		colnames(rmatrix) <- colnames(bmatrix) <- sns
561
+		## put rownames in order of physical position
562
+		ix <- order(position(object)[I])
563
+		I <- I[ix]
564
+		rownames(rmatrix) <- rownames(bmatrix) <- featureNames(object)[I]
565
+		for(j in seq_along(sampleindex)){
566
+			bname <- batch.name[j]
567
+			J <- sampleindex[[j]]
568
+			res <- calculateRTheta(object=object,
569
+					       batch.name=bname,
570
+					       feature.index=I)
571
+			k <- match(sampleNames(object)[J], sns)
572
+			bmatrix[, k] <- res[["baf"]]
573
+			rmatrix[, k] <- res[["lrr"]]
574
+			##colnames(bmatrix)[J] <- colnames(rmatrix)[J] <- sampleNames(object)[J]
575
+		}
576
+		blist[[i]] <- bmatrix
577
+		rlist[[i]] <- rmatrix
578
+	}
579
+	res <- list(baf=blist,
580
+		    lrr=rlist)
581
+	return(res)
582
+}
659 583
 
660
-		  dimnames(bf) <- dimnames(lrr) <- dns
661
-		  res <- list(baf=bf,
662
-			      lrr=lrr)
663
-		  return(res)
664
-	  })
584
+##setMethod("calculateRBaf", signature(object="CNSet"),
585
+##	  function(object, batch.name){
586
+##		  all.autosomes <- all(chromosome(object) < 23)
587
+##		  if(!all.autosomes){
588
+##			  stop("method currently only defined for chromosomes 1-22")
589
+##		  }
590
+##		  if(missing(batch.name)) batch.name <- batchNames(object)[1]
591
+##		  stopifnot(batch.name %in% batchNames(object))
592
+##		  if(length(batch.name) > 1){
593
+##			  warning("only the first batch in batch.name processed")
594
+##			  batch.name <- batch.name[1]
595
+##		  }
596
+##		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)
597
+##		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)
598
+##		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)
599
+##
600
+##		  J <- which(batch(object) == batch.name)
601
+##
602
+##		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
603
+##		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
604
+##		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
605
+##
606
+##		  a <- A(object)[, J, drop=FALSE]
607
+##		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
608
+##		  is.np <- !isSnp(object)
609
+##		  ##b[is.np, ] <- a[is.np, ]
610
+##		  b[is.np, ] <- 0L
611
+##		  dns <- dimnames(a)
612
+##		  dimnames(a) <- dimnames(b) <- NULL
613
+##		  obs.theta <- atan2(b, a)*2/pi
614
+##
615
+##		  lessAA <- obs.theta < theta.aa
616
+##		  lessAB <- obs.theta < theta.ab
617
+##		  lessBB <- obs.theta < theta.bb
618
+##		  grAA <- !lessAA
619
+##		  grAB <- !lessAB
620
+##		  grBB <- !lessBB
621
+##		  not.na <- !is.na(theta.aa)
622
+##		  I1 <- grAA & lessAB & not.na
623
+##		  I2 <- grAB & lessBB & not.na
624
+##
625
+##		  bf <- matrix(NA, nrow(object), ncol(a))
626
+##		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
627
+##		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
628
+##		  bf[lessAA] <- 0
629
+##		  bf[grBB] <- 1
630
+##
631
+##		  r.expected <- matrix(NA, nrow(object), ncol(a))
632
+##		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), length(J), byrow=FALSE)
633
+##		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), length(J), byrow=FALSE)
634
+##		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), length(J), byrow=FALSE)
635
+##		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
636
+##		  obs.r <- a+b
637
+##
638
+##		  lessAA <- lessAA & not.na
639
+##		  grBB <- grBB & not.na
640
+##		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
641
+##		  r.expected[I1] <- tmp
642
+##		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
643
+##		  r.expected[I2] <- tmp
644
+##		  r.expected[lessAA] <- r.aa[lessAA]
645
+##		  r.expected[grBB] <- r.bb[grBB]
646
+##		  index.np <- which(is.np)
647
+##		  if(length(index.np) > 0){
648
+##			  a.np <- A(object)[index.np, J, drop=FALSE]
649
+##			  meds <- rowMedians(a.np, na.rm=TRUE)
650
+##			  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
651
+##		  }
652
+##		  lrr <- log2(obs.r/r.expected)
653
+##
654
+##		  dimnames(bf) <- dimnames(lrr) <- dns
655
+##		  res <- list(baf=bf,
656
+##			      lrr=lrr)
657
+##		  return(res)
658
+##	  })
Browse code

Merge branch 'collab'

* collab:
bump version
made cnSetExample smaller. Fix notes
Trying to revert bad commit
remove cn-functions. update description
comment most of cn-functions.r
Resaved rdas
update data/cnSetExample.rda and data/cnSetExample2.rda
bump version
coercion method from CNSet to oligoSnpSet makes integer matrices of BAFs and lrr's
import ff_or_matrix from oligoClasses. bump dependency on oligoClasses version. Use library(oligoClasses) in some of the crlmm examples.
Cleaning pkg loading process: work still required
move Biobase and methods to imports

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

Rob Scharp authored on 23/03/2012 03:34:50
Showing1 changed files
... ...
@@ -11,6 +11,8 @@ cnSet2oligoSnpSet <- function(object){
11 11
 	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
12 12
 	if(is.lds) stopifnot(isPackageLoaded("ff"))
13 13
 	b.r <- calculateRBaf(object)
14
+	b <- integerMatrix(b.r[[1]], 1000)
15
+	r <- integerMatrix(b.r[[2]], 100)
14 16
 ##	if(is.lds){
15 17
 ##		## initialize a big matrix for raw copy number
16 18
 ##		message("creating an ff object for storing total copy number")
... ...
@@ -438,36 +440,36 @@ setMethod("calculatePosteriorMean", signature(object="CNSet"),
438 440
 		  return(pm)
439 441
 	  })
440 442
 
441
-		  .bivariateCenter <- function(nu, phi){
442
-			  ##  lexical scope for mus, CA, CB
443
-			  if(CA <= 2 & CB <= 2 & (CA+CB) < 4){
444
-				  mus[,1, ] <- log2(nu[, 1, ] + CA *
445
-						    phi[, 1, ])
446
-				  mus[,2, ] <- log2(nu[, 2, ] + CB *
447
-						    phi[, 2, ])
448
-			  } else { ## CA > 2
449
-				  if(CA > 2){
450
-					  theta <- pi/4*Sigma[,2,]
451
-					  shiftA <- CA/4*phi[, 1, ] * cos(theta)
452
-					  shiftB <- CA/4*phi[, 1, ] * sin(theta)
453
-					  mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA)
454
-					  mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB)
455
-				  }
456
-				  if(CB > 2){
457
-					  ## CB > 2
458
-					  theta <- pi/2-pi/4*Sigma[,2,]
459
-					  shiftA <- CB/4*phi[, 2, ] * cos(theta)
460
-					  shiftB <- CB/4*phi[, 2, ] * sin(theta)
461
-					  mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA)
462
-					  mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB)
463
-				  }
464
-				  if(CA == 2 & CB == 2){
465
-					  mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,])
466
-					  mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,])
467
-				  }
468
-			  }
469
-			  mus
470
-		  }
443
+##.bivariateCenter <- function(nu, phi){
444
+##			  ##  lexical scope for mus, CA, CB
445
+##			  if(CA <= 2 & CB <= 2 & (CA+CB) < 4){
446
+##				  mus[,1, ] <- log2(nu[, 1, ] + CA *
447
+##						    phi[, 1, ])
448
+##				  mus[,2, ] <- log2(nu[, 2, ] + CB *
449
+##						    phi[, 2, ])
450
+##			  } else { ## CA > 2
451
+##				  if(CA > 2){
452
+##					  theta <- pi/4*Sigma[,2,]
453
+##					  shiftA <- CA/4*phi[, 1, ] * cos(theta)
454
+##					  shiftB <- CA/4*phi[, 1, ] * sin(theta)
455
+##					  mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA)
456
+##					  mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB)
457
+##				  }
458
+##				  if(CB > 2){
459
+##					  ## CB > 2
460
+##					  theta <- pi/2-pi/4*Sigma[,2,]
461
+##					  shiftA <- CB/4*phi[, 2, ] * cos(theta)
462
+##					  shiftB <- CB/4*phi[, 2, ] * sin(theta)
463
+##					  mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA)
464
+##					  mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB)
465
+##				  }
466
+##				  if(CA == 2 & CB == 2){
467
+##					  mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,])
468
+##					  mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,])
469
+##				  }
470
+##			  }
471
+##			  mus
472
+##		  }
471 473
 
472 474
 ## for a given copy number, return a named list of bivariate normal prediction regions
473 475
 ##   - elements of list are named by genotype
Browse code

Merge branch 'collab'

* collab:
update cnSet objects in data/
bump version
update version dependencies
coerce copy number and bafs to integer matrix in setAs(CNSet, oligoSnpSet)

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

Rob Scharp authored on 10/03/2012 00:05:37
Showing1 changed files
... ...
@@ -42,17 +42,16 @@ cnSet2oligoSnpSet <- function(object){
42 42
 ##	tcn[tcn > 8] <- 8
43 43
 ##	log.tcn <- log2(tcn)
44 44
 	tmp <- new("oligoSnpSet",
45
-		   ##copyNumber=log.tcn,
46
-		   copyNumber=b.r[[2]],
47
-		   ##baf=b.r[[1]],
45
+		   copyNumber=integerMatrix(b.r[[2]], scale=100),
48 46
 		   call=calls(object),
49 47
 		   callProbability=snpCallProbability(object),
50 48
 		   annotation=annotation(object),
51 49
 		   featureData=featureData(object),
52 50
 		   phenoData=phenoData(object),
53 51
 		   experimentData=experimentData(object),
54
-		   protocolData=protocolData(object))
55
-	tmp <- assayDataElementReplace(tmp, "baf", b.r[[1]])
52
+		   protocolData=protocolData(object),
53
+		   is.integer=FALSE)
54
+	tmp <- assayDataElementReplace(tmp, "baf", integerMatrix(b.r[[1]], scale=1000))
56 55
 	return(tmp)
57 56
 }
58 57
 
Browse code

Merge branch 'collab'

* collab:
coercion of CNSet to oligoSnpSet uses lrr/baf instead of ca+cb/genotypes

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

Rob Scharp authored on 05/02/2012 15:20:45
Showing1 changed files
... ...
@@ -1,6 +1,62 @@
1 1
 setMethod("posteriorMean", signature(object="CNSet"), function(object) assayDataElement(object, "posteriorMean"))
2 2
 setReplaceMethod("posteriorMean", signature(object="CNSet", value="matrix"), function(object, value) assayDataElementReplace(object, "posteriorMean", value))
3 3
 
4
+setAs("CNSet", "oligoSnpSet", function(from, to){
5
+	cnSet2oligoSnpSet(from)
6
+})
7
+
8
+cnSet2oligoSnpSet <- function(object){
9
+	row.index <- seq_len(nrow(object))
10
+	col.index <- seq_len(ncol(object))
11
+	is.lds <- ifelse(is(calls(object), "ff_matrix") | is(calls(object), "ffdf"), TRUE, FALSE)
12
+	if(is.lds) stopifnot(isPackageLoaded("ff"))
13
+	b.r <- calculateRBaf(object)
14
+##	if(is.lds){
15
+##		## initialize a big matrix for raw copy number
16
+##		message("creating an ff object for storing total copy number")
17
+##		tcn <- initializeBigMatrix(name="total_cn", nrow(object), ncol(object), vmode="double")
18
+##		for(j in 1:ncol(object)){
19
+##			tcn[, j] <- totalCopynumber(object, i=row.index, j=j)
20
+##		}
21
+##	} else {
22
+##		if(ncol(object) > 5){
23
+##			##this can be memory intensive, so we try to be careful
24
+##			col.index <- splitIndicesByLength(seq(length=ncol(object)), 5)
25
+##			tcn <- matrix(NA, nrow(object), ncol(object))
26
+##			dimnames(tcn) <- list(featureNames(object), sampleNames(object))
27
+##			rows <- 1:nrow(object)
28
+##			for(i in seq_along(col.index)){
29
+##				cat(".")
30
+##				j <- col.index[[i]]
31
+##				cnSet <- object[, j]
32
+##				tcn[, j] <- totalCopynumber(cnSet, i=row.index, j=1:ncol(cnSet))
33
+##				rm(cnSet); gc()
34
+##			}
35
+##			cat("\n")
36
+##		} else {
37
+##			tcn <- totalCopynumber(object, i=row.index, j=col.index)
38
+##		}
39
+##	}
40
+##	message("Transforming copy number to log2 scale")
41
+##	tcn[tcn < 0.1] <- 0.1
42
+##	tcn[tcn > 8] <- 8
43
+##	log.tcn <- log2(tcn)
44
+	tmp <- new("oligoSnpSet",
45
+		   ##copyNumber=log.tcn,
46
+		   copyNumber=b.r[[2]],
47
+		   ##baf=b.r[[1]],
48
+		   call=calls(object),
49
+		   callProbability=snpCallProbability(object),
50
+		   annotation=annotation(object),
51
+		   featureData=featureData(object),
52
+		   phenoData=phenoData(object),
53
+		   experimentData=experimentData(object),
54
+		   protocolData=protocolData(object))
55
+	tmp <- assayDataElementReplace(tmp, "baf", b.r[[1]])
56
+	return(tmp)
57
+}
58
+
59
+
4 60
 linearParamElementReplace <- function(obj, elt, value) {
5 61
     storage.mode <- storageMode(batchStatistics(obj))
6 62
     switch(storage.mode,
Browse code

Compute atan2(0, A)*2/pi for nonpolymorphic markers

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

Rob Scharp authored on 17/01/2012 19:11:45
Showing1 changed files
... ...
@@ -554,7 +554,8 @@ setMethod("calculateRBaf", signature(object="CNSet"),
554 554
 		  a <- A(object)[, J, drop=FALSE]
555 555
 		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
556 556
 		  is.np <- !isSnp(object)
557
-		  b[is.np, ] <- a[is.np, ]
557
+		  ##b[is.np, ] <- a[is.np, ]
558
+		  b[is.np, ] <- 0L
558 559
 		  dns <- dimnames(a)
559 560
 		  dimnames(a) <- dimnames(b) <- NULL
560 561
 		  obs.theta <- atan2(b, a)*2/pi
Browse code

fix dimension of a.np in calculateRBaf

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

Rob Scharp authored on 29/11/2011 20:01:20
Showing1 changed files
... ...
@@ -592,9 +592,11 @@ setMethod("calculateRBaf", signature(object="CNSet"),
592 592
 		  r.expected[grBB] <- r.bb[grBB]
593 593
 
594 594
 		  index.np <- which(is.np)
595
-		  a.np <- A(object)[index.np, ]
596
-		  meds <- rowMedians(a.np, na.rm=TRUE)
597
-		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
595
+		  if(length(index.np) > 0){
596
+			  a.np <- A(object)[index.np, J]
597
+			  meds <- rowMedians(a.np, na.rm=TRUE)
598
+			  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
599
+		  }
598 600
 		  lrr <- log2(obs.r/r.expected)
599 601
 
600 602
 		  dimnames(bf) <- dimnames(lrr) <- dns
Browse code

Fix 'nonconformable arrays' bug in calculateRBaf

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

Rob Scharp authored on 29/11/2011 20:01:16
Showing1 changed files
... ...
@@ -576,9 +576,9 @@ setMethod("calculateRBaf", signature(object="CNSet"),
576 576
 		  bf[grBB] <- 1
577 577
 
578 578
 		  r.expected <- matrix(NA, nrow(object), ncol(a))
579
-		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
580
-		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
581
-		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
579
+		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), length(J), byrow=FALSE)
580
+		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), length(J), byrow=FALSE)
581
+		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), length(J), byrow=FALSE)
582 582
 		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
583 583
 		  obs.r <- a+b
584 584
 
Browse code

update calculateRTheta for NP markers

use atan2(x,x)*2/pi=1

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

Rob Scharp authored on 24/11/2011 12:50:18
Showing1 changed files
... ...
@@ -519,6 +519,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
519 519
 		  ##predictRegion is an argument of ABpanel
520 520
 		  xyplot(x, df, predictRegion=predictRegion, ...)
521 521
 	  })
522
+
522 523
 setMethod("xyplot", signature(x="formula", data="CNSet"),
523 524
 	  function(x, data, ...){
524 525
 		  if("predictRegion" %in% names(list(...))){
... ...
@@ -552,6 +553,8 @@ setMethod("calculateRBaf", signature(object="CNSet"),
552 553
 
553 554
 		  a <- A(object)[, J, drop=FALSE]
554 555
 		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
556
+		  is.np <- !isSnp(object)
557
+		  b[is.np, ] <- a[is.np, ]
555 558
 		  dns <- dimnames(a)
556 559
 		  dimnames(a) <- dimnames(b) <- NULL
557 560
 		  obs.theta <- atan2(b, a)*2/pi
... ...
@@ -588,7 +591,7 @@ setMethod("calculateRBaf", signature(object="CNSet"),
588 591
 		  r.expected[lessAA] <- r.aa[lessAA]
589 592
 		  r.expected[grBB] <- r.bb[grBB]
590 593
 
591
-		  index.np <- which(!isSnp(object))
594
+		  index.np <- which(is.np)
592 595
 		  a.np <- A(object)[index.np, ]
593 596
 		  meds <- rowMedians(a.np, na.rm=TRUE)
594 597
 		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
Browse code

Check that only autosomes are in object passed to calculateRBaf

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

Rob Scharp authored on 01/11/2011 16:32:14
Showing1 changed files
... ...
@@ -530,68 +530,72 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),
530 530
 
531 531
 setMethod("calculateRBaf", signature(object="CNSet"),
532 532
 	  function(object, batch.name){
533
-	if(missing(batch.name)) batch.name <- batchNames(object)[1]
534
-	stopifnot(batch.name %in% batchNames(object))
535
-	if(length(batch.name) > 1){
536
-		warning("only the first batch in batch.name processed")
537
-		batch.name <- batch.name[1]
538
-	}
539
-	RTheta.aa <- calculateRTheta(object, "AA", batch.name)
540
-	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
541
-	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
533
+		  all.autosomes <- all(chromosome(object) < 23)
534
+		  if(!all.autosomes){
535
+			  stop("method currently only defined for chromosomes 1-22")
536
+		  }
537
+		  if(missing(batch.name)) batch.name <- batchNames(object)[1]
538
+		  stopifnot(batch.name %in% batchNames(object))
539
+		  if(length(batch.name) > 1){
540
+			  warning("only the first batch in batch.name processed")
541
+			  batch.name <- batch.name[1]
542
+		  }
543
+		  RTheta.aa <- calculateRTheta(object, "AA", batch.name)
544
+		  RTheta.ab <- calculateRTheta(object, "AB", batch.name)
545
+		  RTheta.bb <- calculateRTheta(object, "BB", batch.name)
542 546
 
543
-	J <- which(batch(object) == batch.name)
547
+		  J <- which(batch(object) == batch.name)
544 548
 
545
-	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
546
-	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
547
-	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
549
+		  theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
550
+		  theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
551
+		  theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
548 552
 
549
-	a <- A(object)[, J, drop=FALSE]
550
-	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
551
-	dns <- dimnames(a)
552
-	dimnames(a) <- dimnames(b) <- NULL
553
-	obs.theta <- atan2(b, a)*2/pi
553
+		  a <- A(object)[, J, drop=FALSE]
554
+		  b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
555
+		  dns <- dimnames(a)
556
+		  dimnames(a) <- dimnames(b) <- NULL
557
+		  obs.theta <- atan2(b, a)*2/pi
554 558
 
555
-	lessAA <- obs.theta < theta.aa
556
-	lessAB <- obs.theta < theta.ab
557
-	lessBB <- obs.theta < theta.bb
558
-	grAA <- !lessAA
559
-	grAB <- !lessAB
560
-	grBB <- !lessBB
561
-	not.na <- !is.na(theta.aa)
562
-	I1 <- grAA & lessAB & not.na
563
-	I2 <- grAB & lessBB & not.na
559
+		  lessAA <- obs.theta < theta.aa
560
+		  lessAB <- obs.theta < theta.ab
561
+		  lessBB <- obs.theta < theta.bb
562
+		  grAA <- !lessAA
563
+		  grAB <- !lessAB
564
+		  grBB <- !lessBB
565
+		  not.na <- !is.na(theta.aa)
566
+		  I1 <- grAA & lessAB & not.na
567
+		  I2 <- grAB & lessBB & not.na
564 568
 
565
-	bf <- matrix(NA, nrow(object), ncol(a))
566
-	bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
567
-	bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
568
-	bf[lessAA] <- 0
569
-	bf[grBB] <- 1
569
+		  bf <- matrix(NA, nrow(object), ncol(a))
570
+		  bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
571
+		  bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
572
+		  bf[lessAA] <- 0
573
+		  bf[grBB] <- 1
570 574
 
571
-	r.expected <- matrix(NA, nrow(object), ncol(a))
572
-	r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
573
-	r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
574
-	r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
575
-	rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
576
-	obs.r <- a+b
575
+		  r.expected <- matrix(NA, nrow(object), ncol(a))
576
+		  r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
577
+		  r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
578
+		  r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
579
+		  rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
580
+		  obs.r <- a+b
577 581
 
578
-	lessAA <- lessAA & not.na
579
-	grBB <- grBB & not.na
580
-	tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
581
-	r.expected[I1] <- tmp
582
-	tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
583
-	r.expected[I2] <- tmp
584
-	r.expected[lessAA] <- r.aa[lessAA]
585
-	r.expected[grBB] <- r.bb[grBB]
582
+		  lessAA <- lessAA & not.na
583
+		  grBB <- grBB & not.na
584
+		  tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
585
+		  r.expected[I1] <- tmp
586
+		  tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
587
+		  r.expected[I2] <- tmp
588
+		  r.expected[lessAA] <- r.aa[lessAA]
589
+		  r.expected[grBB] <- r.bb[grBB]
586 590
 
587
-	index.np <- which(!isSnp(object))
588
-	a.np <- A(object)[index.np, ]
589
-	meds <- rowMedians(a.np, na.rm=TRUE)
590
-	r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
591
-	lrr <- log2(obs.r/r.expected)
591
+		  index.np <- which(!isSnp(object))
592
+		  a.np <- A(object)[index.np, ]
593
+		  meds <- rowMedians(a.np, na.rm=TRUE)
594
+		  r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
595
+		  lrr <- log2(obs.r/r.expected)
592 596
 
593
-	dimnames(bf) <- dimnames(lrr) <- dns
594
-	res <- list(baf=bf,
595
-		    lrr=lrr)
596
-	return(res)
597
-})
597
+		  dimnames(bf) <- dimnames(lrr) <- dns
598
+		  res <- list(baf=bf,
599
+			      lrr=lrr)
600
+		  return(res)
601
+	  })
Browse code

bug fig in calculateRTheta

theta.aa, theta.ab, theta.bb should have columns equal to the number of samples in the batch

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

Rob Scharp authored on 01/11/2011 16:31:57
Showing1 changed files
... ...
@@ -540,11 +540,12 @@ setMethod("calculateRBaf", signature(object="CNSet"),
540 540
 	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
541 541
 	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
542 542
 
543
-	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), ncol(object), byrow=FALSE)
544
-	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), ncol(object), byrow=FALSE)
545
-	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), ncol(object), byrow=FALSE)
546
-
547 543
 	J <- which(batch(object) == batch.name)
544
+
545
+	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), length(J), byrow=FALSE)
546
+	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), length(J), byrow=FALSE)
547
+	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), length(J), byrow=FALSE)
548
+
548 549
 	a <- A(object)[, J, drop=FALSE]
549 550
 	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
550 551
 	dns <- dimnames(a)
Browse code

calculateRBaf adds dimnames to the emission probability array just prior to its return.

previously we assigned the dimnames of the object 'a', but this was NULL.

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

Rob Scharp authored on 04/10/2011 16:12:43
Showing1 changed files
... ...
@@ -547,6 +547,7 @@ setMethod("calculateRBaf", signature(object="CNSet"),
547 547
 	J <- which(batch(object) == batch.name)
548 548
 	a <- A(object)[, J, drop=FALSE]
549 549
 	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
550
+	dns <- dimnames(a)
550 551
 	dimnames(a) <- dimnames(b) <- NULL
551 552
 	obs.theta <- atan2(b, a)*2/pi
552 553
 
... ...
@@ -588,7 +589,7 @@ setMethod("calculateRBaf", signature(object="CNSet"),
588 589
 	r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
589 590
 	lrr <- log2(obs.r/r.expected)
590 591
 
591
-	dimnames(bf) <- dimnames(lrr) <- dimnames(a)
592
+	dimnames(bf) <- dimnames(lrr) <- dns
592 593
 	res <- list(baf=bf,
593 594
 		    lrr=lrr)
594 595
 	return(res)
Browse code

crlmm on github

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

Rob Scharp authored on 04/10/2011 16:11:02
Showing1 changed files
... ...
@@ -20,7 +20,6 @@ linearParamElementReplace <- function(obj, elt, value) {
20 20
 }
21 21
 
22 22
 ## parameters
23
-
24 23
 ## allele A
25 24
 ##   autosome SNPs
26 25
 ##   autosome NPs
Browse code

Add calculateRBaf method

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

Rob Scharp authored on 01/10/2011 04:50:37
Showing1 changed files
... ...
@@ -528,3 +528,69 @@ setMethod("xyplot", signature(x="formula", data="CNSet"),
528 528
 			  callNextMethod()
529 529
 		  }
530 530
 })
531
+
532
+setMethod("calculateRBaf", signature(object="CNSet"),
533
+	  function(object, batch.name){
534
+	if(missing(batch.name)) batch.name <- batchNames(object)[1]
535
+	stopifnot(batch.name %in% batchNames(object))
536
+	if(length(batch.name) > 1){
537
+		warning("only the first batch in batch.name processed")
538
+		batch.name <- batch.name[1]
539
+	}
540
+	RTheta.aa <- calculateRTheta(object, "AA", batch.name)
541
+	RTheta.ab <- calculateRTheta(object, "AB", batch.name)
542
+	RTheta.bb <- calculateRTheta(object, "BB", batch.name)
543
+
544
+	theta.aa <- matrix(RTheta.aa[, "theta"], nrow(object), ncol(object), byrow=FALSE)
545
+	theta.ab <- matrix(RTheta.ab[, "theta"], nrow(object), ncol(object), byrow=FALSE)
546
+	theta.bb <- matrix(RTheta.bb[, "theta"], nrow(object), ncol(object), byrow=FALSE)
547
+
548
+	J <- which(batch(object) == batch.name)
549
+	a <- A(object)[, J, drop=FALSE]
550
+	b <- B(object)[, J, drop=FALSE] ## NA's for b where nonpolymorphic
551
+	dimnames(a) <- dimnames(b) <- NULL
552
+	obs.theta <- atan2(b, a)*2/pi
553
+
554
+	lessAA <- obs.theta < theta.aa
555
+	lessAB <- obs.theta < theta.ab
556
+	lessBB <- obs.theta < theta.bb
557
+	grAA <- !lessAA
558
+	grAB <- !lessAB
559
+	grBB <- !lessBB
560
+	not.na <- !is.na(theta.aa)
561
+	I1 <- grAA & lessAB & not.na
562
+	I2 <- grAB & lessBB & not.na
563
+
564
+	bf <- matrix(NA, nrow(object), ncol(a))
565
+	bf[I1] <- 0.5 * ((obs.theta-theta.aa)/(theta.ab-theta.aa))[I1]
566
+	bf[I2] <- (.5 * (obs.theta - theta.ab) / (theta.bb - theta.ab))[I2] + 0.5
567
+	bf[lessAA] <- 0
568
+	bf[grBB] <- 1
569
+
570
+	r.expected <- matrix(NA, nrow(object), ncol(a))
571
+	r.aa <- matrix(RTheta.aa[, "R"], nrow(object), ncol(object), byrow=FALSE)
572
+	r.ab <- matrix(RTheta.ab[, "R"], nrow(object), ncol(object), byrow=FALSE)
573
+	r.bb <- matrix(RTheta.bb[, "R"], nrow(object), ncol(object), byrow=FALSE)
574
+	rm(RTheta.aa, RTheta.ab, RTheta.bb); gc()
575
+	obs.r <- a+b
576
+
577
+	lessAA <- lessAA & not.na
578
+	grBB <- grBB & not.na
579
+	tmp <- ((obs.theta - theta.aa) * (r.ab-r.aa)/(theta.ab-theta.aa))[I1] + r.aa[I1]
580
+	r.expected[I1] <- tmp
581
+	tmp <- ((obs.theta - theta.ab) * (r.bb - r.ab)/(theta.bb-theta.ab))[I2] + r.ab[I2]
582
+	r.expected[I2] <- tmp
583
+	r.expected[lessAA] <- r.aa[lessAA]
584
+	r.expected[grBB] <- r.bb[grBB]
585
+
586
+	index.np <- which(!isSnp(object))
587
+	a.np <- A(object)[index.np, ]
588
+	meds <- rowMedians(a.np, na.rm=TRUE)
589
+	r.expected[index.np, ] <- matrix(meds, length(index.np), ncol(a.np))
590
+	lrr <- log2(obs.r/r.expected)
591
+
592
+	dimnames(bf) <- dimnames(lrr) <- dimnames(a)
593
+	res <- list(baf=bf,
594
+		    lrr=lrr)
595
+	return(res)
596
+})
Browse code

Use original version of bivariateCenter.

Defining the bivariate center to acknowledge background through correlation did not help the signal to noise ratio in the trisomy data. This code is now in .bivariateCenter.

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

Rob Scharp authored on 01/10/2011 04:49:31
Showing1 changed files
... ...
@@ -384,6 +384,36 @@ setMethod("calculatePosteriorMean", signature(object="CNSet"),
384 384
 		  return(pm)
385 385
 	  })
386 386
 
387
+		  .bivariateCenter <- function(nu, phi){
388
+			  ##  lexical scope for mus, CA, CB
389
+			  if(CA <= 2 & CB <= 2 & (CA+CB) < 4){
390
+				  mus[,1, ] <- log2(nu[, 1, ] + CA *
391
+						    phi[, 1, ])
392
+				  mus[,2, ] <- log2(nu[, 2, ] + CB *
393
+						    phi[, 2, ])
394
+			  } else { ## CA > 2
395
+				  if(CA > 2){
396
+					  theta <- pi/4*Sigma[,2,]
397
+					  shiftA <- CA/4*phi[, 1, ] * cos(theta)
398
+					  shiftB <- CA/4*phi[, 1, ] * sin(theta)
399
+					  mus[, 1, ] <- log2(nu[, 1, ] + 2 * phi[,1,]+shiftA)
400
+					  mus[, 2, ] <- log2(nu[, 2, ] + CB *phi[,2,]+shiftB)
401
+				  }
402
+				  if(CB > 2){
403
+					  ## CB > 2
404
+					  theta <- pi/2-pi/4*Sigma[,2,]
405
+					  shiftA <- CB/4*phi[, 2, ] * cos(theta)
406
+					  shiftB <- CB/4*phi[, 2, ] * sin(theta)
407
+					  mus[, 1, ] <- log2(nu[, 1, ] + CA*phi[,1,]+shiftA)
408
+					  mus[, 2, ] <- log2(nu[, 2, ]+ 2*phi[,2,]+shiftB)
409
+				  }
410
+				  if(CA == 2 & CB == 2){
411
+					  mus[, 1, ] <- log2(nu[, 1, ] + 1/2*CA*phi[,1,])
412
+					  mus[, 2, ] <- log2(nu[, 2, ]+ 1/2*CB*phi[,2,])
413
+				  }
414
+			  }
415
+			  mus
416
+		  }
387 417
 
388 418
 ## for a given copy number, return a named list of bivariate normal prediction regions
389 419
 ##   - elements of list are named by genotype
... ...
@@ -428,9 +458,11 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
428 458
 				 bnms))
429 459
 		  bivariateCenter <- function(nu, phi){
430 460
 			  ##  lexical scope for mus, CA, CB
431
-			  mus[,1, ] <- log2(nu[, 1, ] + CA * phi[, 1, ])
432
-			  mus[,2, ] <- log2(nu[, 2, ] + CB * phi[, 2, ])
433
-			  mus
461
+			  mus[,1, ] <- log2(nu[, 1, ] + CA *
462
+					    phi[, 1, ])
463
+			  mus[,2, ] <- log2(nu[, 2, ] + CB *
464
+					    phi[, 2, ])
465
+			  return(mus)
434 466
 		  }
435 467
 		  np.index <- which(!isSnp(object))
436 468
 		  for(i in seq_along(copyNumber)){
... ...
@@ -452,14 +484,20 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
452 484
 					  Sigma[, 1, ] <- getVar(object, "tau2A.AA")
453 485
 					  Sigma[np.index, 2, ] <- NA
454 486
 				  }
455
-				  res[[gt]]$mu <- bivariateCenter(nu, phi)
487
+				  res[[gt]]$mu <- bivariateCenter(nu,phi)
488
+				  ## adjust correlation
489
+##				  if(CA == 0 | CB == 0){
490
+##					  Sigma[,2,] <- 0
491
+##				  }
456 492
 				  res[[gt]]$cov <- Sigma
457 493
 			  }
458 494
 		  }
495
+		  res <- as(res, "PredictionRegion")
459 496
 		  return(res)
460 497
 	  })
461 498
 
462 499
 
500
+
463 501
 setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="list"),
464 502
 	  function(x, data, predictRegion, ...){
465 503
 		  fns <- featureNames(data)
... ...
@@ -479,7 +517,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
479 517
 			  x$cov <- x$cov[, , batch.index, drop=FALSE]
480 518
 			  return(x)
481 519
 		  }, bns=bns)
482
-		  ##df <- as.data.frame(data)
520
+		  ##predictRegion is an argument of ABpanel
483 521
 		  xyplot(x, df, predictRegion=predictRegion, ...)
484 522
 	  })
485 523
 setMethod("xyplot", signature(x="formula", data="CNSet"),
Browse code

Update posteriorProbability for handling nonpolymorphic markers

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

Rob Scharp authored on 01/10/2011 04:49:02
Showing1 changed files
... ...
@@ -298,36 +298,69 @@ setMethod("posteriorProbability", signature(object="CNSet"),
298 298
 		  prob <- array(NA, dim=c(nrow(object), ncol(object), length(copyNumber)),
299 299
 				dimnames=list(NULL,NULL, paste("copynumber",copyNumber, sep="")))
300 300
 		  bns <- batchNames(object)
301
-		  for(i in seq_along(copyNumber)){
302
-			  G <- genotypes(copyNumber[i])
303
-			  P <- array(NA, dim=c(nrow(object), ncol(object), length(G)),
304
-				     dimnames=list(NULL,NULL,G))
305
-			  for(g in seq_along(G)){
306
-				  gt <- G[g]
307
-				  for(j in seq_along(bns)){
308
-					  this.batch <- bns[j]
309
-					  sample.index <- which(batch(object) == this.batch)
310
-					  mu <- predictRegion[[gt]]$mu[, , this.batch, drop=FALSE]
311
-					  dim(mu) <- dim(mu)[1:2]
312
-					  if(all(is.na(mu))){
313
-						  P[, sample.index, gt] <- NA
314
-					  } else {
315
-						  Sigma <- predictRegion[[gt]]$cov[, , this.batch, drop=FALSE]
316
-						  dim(Sigma) <- dim(Sigma)[1:2]
317
-						  P[, sample.index, gt] <- dbvn(x=logI[, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
301
+		  snp.index <- which(isSnp(object))
302
+		  if(length(snp.index) > 0){
303
+			  for(i in seq_along(copyNumber)){
304
+				  G <- genotypes(copyNumber[i])
305
+				  P <- array(NA, dim=c(length(snp.index), ncol(object), length(G)),
306
+					     dimnames=list(NULL,NULL,G))
307
+				  for(g in seq_along(G)){
308
+					  gt <- G[g]
309
+					  for(j in seq_along(bns)){
310
+						  this.batch <- bns[j]
311
+						  sample.index <- which(batch(object) == this.batch)
312
+						  mu <- predictRegion[[gt]]$mu[snp.index, , this.batch, drop=FALSE]
313
+						  dim(mu) <- dim(mu)[1:2]
314
+						  if(all(is.na(mu))){
315
+							  P[, sample.index, gt] <- NA
316
+						  } else {
317
+							  Sigma <- predictRegion[[gt]]$cov[snp.index, , this.batch, drop=FALSE]
318
+							  dim(Sigma) <- dim(Sigma)[1:2]
319
+							  P[, sample.index, gt] <- dbvn(x=logI[snp.index, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
320
+						  }
318 321
 					  }
319 322
 				  }
323
+				  ## the marginal probability for the total copy number
324
+				  ##  -- integrate out the probability for the different genotypes
325
+				  for(j in 1:ncol(P)){
326
+					  PP <- P[, j, , drop=FALSE]
327
+					  dim(PP) <- dim(PP)[c(1, 3)]
328
+					  prob[snp.index, j, i] <- rowSums(PP, na.rm=TRUE)
329
+				  }
320 330
 			  }
321
-			  ## the marginal probability for the total copy number
322
-			  ##  -- integrate out the probability for the different genotypes
323
-			  for(j in 1:ncol(P)){
324
-				  PP <- P[, j, , drop=FALSE]
325
-				  dim(PP) <- dim(PP)[c(1, 3)]
326
-				  prob[, j, i] <- rowSums(PP, na.rm=TRUE)
327
-			  }
328
-		  }
329
-
330
-		  ## scale by prior weights and divide by normalizing
331
+		  } ## length(snp.index) > 0
332
+		  ##---------------------------------------------------------------------------
333
+		  ##
334
+		  ##  nonpolymorphic markers
335
+		  ##
336
+		  ##---------------------------------------------------------------------------
337
+		  np.index <- which(!isSnp(object))
338
+		  if(length(np.index) > 0){
339
+			  for(i in seq_along(copyNumber)){
340
+				  G <- genotypes(copyNumber[i], is.snp=FALSE)
341
+				  P <- array(NA, dim=c(length(np.index), ncol(object), length(G)),
342
+					     dimnames=list(NULL,NULL,G))
343
+				  for(g in seq_along(G)){
344
+					  gt <- G[g]
345
+					  for(j in seq_along(bns)){
346
+						  this.batch <- bns[j]
347
+						  sample.index <- which(batch(object) == this.batch)
348
+						  mu <- predictRegion[[gt]]$mu[np.index, 1, this.batch, drop=FALSE]
349
+						  dim(mu) <- dim(mu)[1]
350
+						  if(all(is.na(mu))){
351
+							  P[, sample.index, gt] <- NA
352
+						  } else {
353
+							  Sigma <- predictRegion[[gt]]$cov[np.index, 1, this.batch, drop=FALSE]
354
+							  dim(Sigma) <- dim(Sigma)[1]
355
+							  P[, sample.index, gt] <- dnorm(logI[np.index, sample.index, 1], mean=mu, sd=Sigma)
356
+						  }
357
+					  } ## j in seq_along(bns)
358
+				  } ## g in seq_along(G)
359
+				  ## the marginal probability for the total copy number
360
+				  ##  -- integrate out the probability for the different genotypes
361
+				  prob[np.index, , i] <- P
362
+			  } ## seq_along(copyNumber)
363
+		  } ## length(np.index) > 0
331 364
 		  ##constant.  Probs across copy number states must
332 365
 		  ##sum to one.  nc <- apply(prob, c(1,3), sum,
333 366
 		  ##na.rm=TRUE)
... ...
@@ -399,6 +432,7 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
399 432
 			  mus[,2, ] <- log2(nu[, 2, ] + CB * phi[, 2, ])
400 433
 			  mus
401 434
 		  }
435
+		  np.index <- which(!isSnp(object))
402 436
 		  for(i in seq_along(copyNumber)){
403 437
 			  G <- genotypes(copyNumber[i])
404 438
 			  tmp <- vector("list", length(G))
... ...
@@ -414,6 +448,10 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
414 448
 				  Sigma[, 1, ] <- getVar(object, nma)
415 449
 				  Sigma[, 3, ] <- getVar(object, nmb)
416 450
 				  Sigma[, 2, ] <- getCor(object, gt.corr)
451
+				  if(length(np.index) > 0){
452
+					  Sigma[, 1, ] <- getVar(object, "tau2A.AA")
453
+					  Sigma[np.index, 2, ] <- NA
454
+				  }
417 455
 				  res[[gt]]$mu <- bivariateCenter(nu, phi)
418 456
 				  res[[gt]]$cov <- Sigma
419 457
 			  }
Browse code

Update Rd files. xyplot method checks for existance of 'predictionRange'. If not found, we use callNextMethod() [Should help avoid conflicts with VanillaICE method for xyplot]

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

Rob Scharp authored on 01/10/2011 04:48:57
Showing1 changed files
... ...
@@ -446,5 +446,9 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
446 446
 	  })
447 447
 setMethod("xyplot", signature(x="formula", data="CNSet"),
448 448
 	  function(x, data, ...){
449
-		  xyplotcrlmm(x, data, ...)
449
+		  if("predictRegion" %in% names(list(...))){
450
+			  xyplotcrlmm(x, data, ...)
451
+		  } else{
452
+			  callNextMethod()
453
+		  }
450 454
 })
Browse code

Use factor for snpid in xyplotcrlmm (otherwise plots are out of order

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

Rob Scharp authored on 01/10/2011 04:47:47
Showing1 changed files
... ...
@@ -433,6 +433,7 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
433 433
 			     gt.conf=as.numeric(confs(data)),
434 434
 			     snpid=fns)#, snp=snpId)
435 435
 		  df <- as.data.frame(df)
436
+		  df$snpid <- factor(df$snpid, levels=unique(df$snpid), ordered=TRUE)
436 437
 		  bns <- batchNames(data)
437 438
 		  predictRegion <- lapply(predictRegion, function(x, bns){
438 439
 			  batch.index <- match(bns, dimnames(x$mu)[[3]])
... ...
@@ -440,7 +441,6 @@ setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="lis
440 441
 			  x$cov <- x$cov[, , batch.index, drop=FALSE]
441 442
 			  return(x)
442 443
 		  }, bns=bns)
443
-		  df$snpid <- as.character(df$snpid)
444 444
 		  ##df <- as.data.frame(data)
445 445
 		  xyplot(x, df, predictRegion=predictRegion, ...)
446 446
 	  })
Browse code

Update generics for posteriorProbability and calculatePosteriorMean (exclude weight from generic of calculatePosteriorMean)

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

Rob Scharp authored on 01/10/2011 04:47:35
Showing1 changed files
... ...
@@ -284,7 +284,9 @@ setMethod("totalCopynumber", signature=signature(object="CNSet"),
284 284
 rawCopynumber <- totalCopynumber
285 285
 
286 286
 setMethod("posteriorProbability", signature(object="CNSet"),
287
-	  function(object, predictRegion, copyNumber=0:4){
287
+	  function(object, predictRegion, copyNumber=0:4, w){
288
+		  if(missing(w)) w <- rep(1/length(copyNumber),length(copyNumber))
289
+		  stopifnot(sum(w)==1)
288 290
 		  logI <- array(NA, dim=c(nrow(object), ncol(object), 2), dimnames=list(NULL, NULL, LETTERS[1:2]))
289 291
 		  getIntensity <- function(object){
290 292
 			  logI[, , 1] <- log2(A(object))
... ...
@@ -295,38 +297,56 @@ setMethod("posteriorProbability", signature(object="CNSet"),
295 297
 		  ##gts <- lapply(as.list(0:4), genotypes)
296 298
 		  prob <- array(NA, dim=c(nrow(object), ncol(object), length(copyNumber)),
297 299
 				dimnames=list(NULL,NULL, paste("copynumber",copyNumber, sep="")))
300
+		  bns <- batchNames(object)
298 301
 		  for(i in seq_along(copyNumber)){
299 302
 			  G <- genotypes(copyNumber[i])
300 303
 			  P <- array(NA, dim=c(nrow(object), ncol(object), length(G)),
301 304
 				     dimnames=list(NULL,NULL,G))
302 305
 			  for(g in seq_along(G)){
303 306
 				  gt <- G[g]
304
-				  P[, , gt] <- dbvn(x=logI, mu=predictRegion[[gt]]$mu, Sigma=predictRegion[[gt]]$cov)
307
+				  for(j in seq_along(bns)){
308
+					  this.batch <- bns[j]
309
+					  sample.index <- which(batch(object) == this.batch)
310
+					  mu <- predictRegion[[gt]]$mu[, , this.batch, drop=FALSE]
311
+					  dim(mu) <- dim(mu)[1:2]
312
+					  if(all(is.na(mu))){
313
+						  P[, sample.index, gt] <- NA
314
+					  } else {
315
+						  Sigma <- predictRegion[[gt]]$cov[, , this.batch, drop=FALSE]
316
+						  dim(Sigma) <- dim(Sigma)[1:2]
317
+						  P[, sample.index, gt] <- dbvn(x=logI[, sample.index, , drop=FALSE], mu=mu, Sigma=Sigma)
318
+					  }
319
+				  }
305 320
 			  }
306 321
 			  ## the marginal probability for the total copy number
307 322
 			  ##  -- integrate out the probability for the different genotypes
308 323
 			  for(j in 1:ncol(P)){
309
-				  prob[, j, i] <- rowSums(as.matrix(P[, j, ]), na.rm=TRUE)
324
+				  PP <- P[, j, , drop=FALSE]
325
+				  dim(PP) <- dim(PP)[c(1, 3)]
326
+				  prob[, j, i] <- rowSums(PP, na.rm=TRUE)
310 327
 			  }
311 328
 		  }
312
-		  ## divide by normalizing constant.  Probs across copy number states must sum to one.
313
-		  ##nc <- apply(prob, c(1,3), sum, na.rm=TRUE)
329
+
330
+		  ## scale by prior weights and divide by normalizing
331
+		  ##constant.  Probs across copy number states must
332
+		  ##sum to one.  nc <- apply(prob, c(1,3), sum,
333
+		  ##na.rm=TRUE)
334
+		  wm <- matrix(w, nrow(object), length(copyNumber), byrow=TRUE)
314 335
 		  nc <- matrix(NA, nrow(object), ncol(object))
315 336
 		  for(j in seq_len(ncol(object))){
316
-			  nc <- rowSums(as.matrix(prob[,j, ]), na.rm=TRUE)
317
-			  prob[, j, ] <- prob[, j, ]/nc
337
+			  pp <- prob[, j, ] * wm
338
+			  nc <- rowSums(pp, na.rm=TRUE)
339
+			  prob[, j, ] <- pp/nc
318 340
 		  }
319 341
 		  return(prob)
320 342
 	  })
321 343
 
322 344
 setMethod("calculatePosteriorMean", signature(object="CNSet"),
323
-	  function(object, posteriorProb, copyNumber=0:4, w, ...){
324
-		  if(missing(w)) w <- rep(1/length(copyNumber),length(copyNumber))
325
-		  stopifnot(sum(w)==1)
345
+	  function(object, posteriorProb, copyNumber=0:4, ...){
326 346
 		  stopifnot(dim(posteriorProb)[[3]] == length(copyNumber))
327 347
 		  pm <- matrix(0, nrow(object), ncol(object))
328 348
 		  for(i in seq_along(copyNumber)){
329
-			  pm <- pm + posteriorProb[, , i] * copyNumber[i] * w[i]
349
+			  pm <- pm + posteriorProb[, , i] * copyNumber[i]
330 350
 		  }
331 351
 		  return(pm)
332 352
 	  })
... ...
@@ -404,7 +424,23 @@ setMethod("predictionRegion", signature(object="CNSet", copyNumber="integer"),
404 424
 
405 425
 setMethod("xyplotcrlmm", signature(x="formula", data="CNSet", predictRegion="list"),
406 426
 	  function(x, data, predictRegion, ...){
407
-		  df <- data.frame(A=log2(A(data)), B=log2(B(data)), gt=calls(data), gt.conf=confs(data))#, snp=snpId)
427
+		  fns <- featureNames(data)
428
+		  fns <- matrix(fns, nrow(data), ncol(data), byrow=FALSE)
429
+		  fns <- as.character(fns)
430
+		  df <- list(A=as.numeric(log2(A(data))),
431
+			     B=as.numeric(log2(B(data))),
432
+			     gt=as.integer(calls(data)),
433
+			     gt.conf=as.numeric(confs(data)),
434
+			     snpid=fns)#, snp=snpId)
435
+		  df <- as.data.frame(df)
436
+		  bns <- batchNames(data)
437
+		  predictRegion <- lapply(predictRegion, function(x, bns){
438
+			  batch.index <- match(bns, dimnames(x$mu)[[3]])
439
+			  x$mu <- x$mu[, , batch.index, drop=FALSE]
440
+			  x$cov <- x$cov[, , batch.index, drop=FALSE]
441
+			  return(x)
442
+		  }, bns=bns)
443
+		  df$snpid <- as.character(df$snpid)
408 444
 		  ##df <- as.data.frame(data)
409 445
 		  xyplot(x, df, predictRegion=predictRegion, ...)
410 446
 	  })
Browse code

Update predictionRegion to return arrays in the mu and cov list elements (prediction regions are batch-specific)

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

Rob Scharp authored on 01/10/2011 04:47:24
Showing1 changed files