git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/crlmm@43010 bc3139a8-67e5-0310-9ffc-ced21a209358
Rob Scharp authored on 15/11/2009 10:46:23
| ... | ... |
@@ -312,3 +312,11 @@ is decoded and scanned |
| 312 | 312 |
2009-11-14 B Carvalho - committed version 1.5.2 |
| 313 | 313 |
|
| 314 | 314 |
* code cleanup - removed unneeded C code |
| 315 |
+ |
|
| 316 |
+2009-11-15 R. Scharpf - committed version 1.5.3 |
|
| 317 |
+ |
|
| 318 |
+ * numerous changes to the copy number estimation |
|
| 319 |
+ - ABset and CrlmmSetList classes are gone. Using SnpCallSetPlus and CrlmmSet |
|
| 320 |
+ - Added class SegmentSet that extends CrlmmSet directly. |
|
| 321 |
+ |
|
| 322 |
+ |
| ... | ... |
@@ -1,15 +1,15 @@ |
| 1 | 1 |
Package: crlmm |
| 2 | 2 |
Type: Package |
| 3 | 3 |
Title: Genotype Calling (CRLMM) and Copy Number Analysis tool for Affymetrix SNP 5.0 and 6.0 and Illumina arrays. |
| 4 |
-Version: 1.5.1 |
|
| 5 |
-Date: 2009-10-25 |
|
| 4 |
+Version: 1.5.3 |
|
| 5 |
+Date: 2009-11-15 |
|
| 6 | 6 |
Author: Rafael A Irizarry, Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU> |
| 7 | 7 |
Maintainer: Benilton S Carvalho <bcarvalh@jhsph.edu>, Robert Scharpf <rscharpf@jhsph.edu>, Matt Ritchie <mritchie@wehi.EDU.AU> |
| 8 |
-Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 6.0. |
|
| 8 |
+Description: Faster implementation of CRLMM specific to SNP 5.0 and 6.0 arrays, as well as a copy number tool specific to 5.0, 6.0, and Illumina platforms |
|
| 9 | 9 |
License: Artistic-2.0 |
| 10 | 10 |
Depends: methods, Biobase (>= 2.5.5), R (>= 2.10.0) |
| 11 |
-Imports: affyio, preprocessCore, utils, stats, genefilter, splines, mvtnorm, oligoClasses, ellipse, methods, SNPchip |
|
| 12 |
-Suggests: hapmapsnp5, hapmapsnp6, genomewidesnp5Crlmm (>= 1.0.2),genomewidesnp6Crlmm (>= 1.0.2), snpMatrix, VanillaICE (>= 1.7.8), GGdata |
|
| 11 |
+Imports: affyio, preprocessCore, utils, stats, genefilter, splines, mvtnorm, oligoClasses, ellipse, methods, SNPchip, oligoClasses, VanillaICE (>= 1.9.1) |
|
| 12 |
+Suggests: hapmapsnp5, hapmapsnp6, genomewidesnp5Crlmm (>= 1.0.2),genomewidesnp6Crlmm (>= 1.0.2), snpMatrix, GGdata |
|
| 13 | 13 |
Collate: AllClasses.R |
| 14 | 14 |
AllGenerics.R |
| 15 | 15 |
methods-ABset.R |
| ... | ... |
@@ -1,11 +1,7 @@ |
| 1 |
-## Class definition |
|
| 2 |
-setClass("ABset", contains="eSet",
|
|
| 3 |
- prototype = prototype( |
|
| 4 |
- new("VersionedBiobase",
|
|
| 5 |
- versions=c(classVersion("eSet"), SnpSet="1.0.0"))))
|
|
| 6 |
-setClass("crlmmSet", contains="eSet")
|
|
| 7 |
-setClass("CrlmmSetList", contains="list")
|
|
| 8 |
-setClass("CopyNumberSet", contains="eSet",
|
|
| 9 |
- prototype = prototype( |
|
| 10 |
- new("VersionedBiobase",
|
|
| 11 |
- versions=c(classVersion("eSet"), SnpSet="1.0.0"))))
|
|
| 1 |
+setClass("CopyNumberSet", contains="SnpLevelSet")
|
|
| 2 |
+setClass("CrlmmSet", contains=c("SnpCallSetPlus", "CopyNumberSet"))
|
|
| 3 |
+setClass("SnpCallSetPlusFF", contains="SnpCallSetPlus")
|
|
| 4 |
+setClass("CrlmmSetFF", contains="CrlmmSet")
|
|
| 5 |
+setClass("SegmentSet", contains="CrlmmSet",
|
|
| 6 |
+ representation(emissionPr="array", |
|
| 7 |
+ segmentData="data.frame")) |
| ... | ... |
@@ -3,20 +3,27 @@ setGeneric("B", function(object) standardGeneric("B"))
|
| 3 | 3 |
setGeneric("A<-", function(object, value) standardGeneric("A<-"))
|
| 4 | 4 |
setGeneric("addFeatureAnnotation", function(object, ...) standardGeneric("addFeatureAnnotation"))
|
| 5 | 5 |
setGeneric("B<-", function(object, value) standardGeneric("B<-"))
|
| 6 |
-setGeneric("batch", function(object) standardGeneric("batch"))
|
|
| 7 |
-##setGeneric("calls", function(x) standardGeneric("calls"))
|
|
| 8 | 6 |
setGeneric("confs", function(object) standardGeneric("confs"))
|
| 9 |
-setGeneric("CA", function(object, ...) standardGeneric("CA"))
|
|
| 10 |
-setGeneric("CB", function(object, ...) standardGeneric("CB"))
|
|
| 7 |
+setGeneric("CA", function(object) standardGeneric("CA"))
|
|
| 8 |
+setGeneric("CB", function(object) standardGeneric("CB"))
|
|
| 11 | 9 |
setGeneric("CA<-", function(object, value) standardGeneric("CA<-"))
|
| 12 | 10 |
setGeneric("CB<-", function(object, value) standardGeneric("CB<-"))
|
| 13 |
-##setGeneric("chromosome", function(object) standardGeneric("chromosome"))
|
|
| 14 |
-setGeneric("cnIndex", function(object, ...) standardGeneric("cnIndex"))
|
|
| 15 |
-setGeneric("cnNames", function(object, ...) standardGeneric("cnNames"))
|
|
| 16 |
-##setGeneric("copyNumber", function(object) standardGeneric("copyNumber"))
|
|
| 17 |
-##setGeneric("position", function(object) standardGeneric("position"))
|
|
| 11 |
+setGeneric("emissionPr", function(object) standardGeneric("emissionPr"))
|
|
| 12 |
+setGeneric("emissionPr<-", function(object, value) standardGeneric("emissionPr<-"))
|
|
| 13 |
+setGeneric("getParam", function(object, name, batch) standardGeneric("getParam"))
|
|
| 14 |
+setGeneric("GT<-", function(object, value) standardGeneric("GT<-"))
|
|
| 15 |
+setGeneric("cnIndex", function(object) standardGeneric("cnIndex"))
|
|
| 16 |
+setGeneric("cnNames", function(object) standardGeneric("cnNames"))
|
|
| 17 |
+setGeneric("computeCopynumber", function(object, cnOptions) standardGeneric("computeCopynumber"))
|
|
| 18 |
+setGeneric("computeEmission", function(object, hmmOptions) standardGeneric("computeEmission"))
|
|
| 19 |
+setGeneric("computeHmm", function(object, hmmOptions) standardGeneric("computeHmm"))
|
|
| 20 |
+setGeneric("confs<-", function(object, value) standardGeneric("confs<-"))
|
|
| 21 |
+setGeneric("GT", function(object, ...) standardGeneric("GT"))
|
|
| 18 | 22 |
setGeneric(".harmonizeDimnames", function(object) standardGeneric(".harmonizeDimnames"))
|
| 19 |
-setGeneric("snpIndex", function(object, ...) standardGeneric("snpIndex"))
|
|
| 20 |
-setGeneric("snpNames", function(object, ...) standardGeneric("snpNames"))
|
|
| 23 |
+setGeneric("isSnp", function(object) standardGeneric("isSnp"))
|
|
| 24 |
+setGeneric("pr", function(object, name, batch, value) standardGeneric("pr"))
|
|
| 25 |
+setGeneric("segmentData", function(object) standardGeneric("segmentData"))
|
|
| 26 |
+setGeneric("segmentData<-", function(object, value) standardGeneric("segmentData<-"))
|
|
| 27 |
+setGeneric("snpIndex", function(object) standardGeneric("snpIndex"))
|
|
| 28 |
+setGeneric("snpNames", function(object) standardGeneric("snpNames"))
|
|
| 21 | 29 |
setGeneric("splitByChromosome", function(object, ...) standardGeneric("splitByChromosome"))
|
| 22 |
- |
| 23 | 30 |
deleted file mode 100644 |
| ... | ... |
@@ -1,18 +0,0 @@ |
| 1 |
-setValidity("ABset", function(object) {
|
|
| 2 |
- ##msg <- validMsg(NULL, Biobase:::isValidVersion(object, "CopyNumberSet")) |
|
| 3 |
- msg <- validMsg(NULL, assayDataValidMembers(assayData(object), c("A", "B")))
|
|
| 4 |
- if (is.null(msg)) TRUE else msg |
|
| 5 |
-}) |
|
| 6 |
-setMethod("A", "ABset", function(object) assayData(object)[["A"]])
|
|
| 7 |
-setMethod("B", "ABset", function(object) assayData(object)[["B"]])
|
|
| 8 |
-setReplaceMethod("A", signature(object="ABset", value="matrix"),
|
|
| 9 |
- function(object, value){
|
|
| 10 |
- assayDataElementReplace(object, "A", value) |
|
| 11 |
- }) |
|
| 12 |
-setReplaceMethod("B", signature(object="ABset", value="matrix"),
|
|
| 13 |
- function(object, value){
|
|
| 14 |
- assayDataElementReplace(object, "B", value) |
|
| 15 |
- }) |
|
| 16 |
- |
|
| 17 |
- |
|
| 18 |
- |
| ... | ... |
@@ -3,49 +3,26 @@ setValidity("CopyNumberSet", function(object) {
|
| 3 | 3 |
if (is.null(msg)) TRUE else msg |
| 4 | 4 |
}) |
| 5 | 5 |
##may want to allow thresholding here (... arg) |
| 6 |
-setMethod("CA", "CopyNumberSet", function(object, ...) assayData(object)[["CA"]]/100)
|
|
| 7 |
-setMethod("CB", "CopyNumberSet", function(object, ...) assayData(object)[["CB"]]/100)
|
|
| 8 |
- |
|
| 9 |
-setReplaceMethod("CA", signature(object="CopyNumberSet", value="matrix"), function(object, value){
|
|
| 10 |
- dns <- dimnames(value) |
|
| 11 |
- value <- matrix(as.integer(value*100), nrow(value), ncol(value)) |
|
| 12 |
- dimnames(value) <- dns |
|
| 13 |
- assayDataElementReplace(object, "CA", value) |
|
| 14 |
-}) |
|
| 15 |
- |
|
| 16 |
-setReplaceMethod("CB", signature(object="CopyNumberSet", value="matrix"), function(object, value){
|
|
| 17 |
- dns <- dimnames(value) |
|
| 18 |
- value <- matrix(as.integer(value*100), nrow(value), ncol(value)) |
|
| 19 |
- dimnames(value) <- dns |
|
| 20 |
- assayDataElementReplace(object, "CB", value) |
|
| 21 |
-}) |
|
| 22 |
- |
|
| 6 |
+setMethod("CA", "CopyNumberSet", function(object) assayData(object)[["CA"]]/100)
|
|
| 7 |
+setMethod("CB", "CopyNumberSet", function(object) assayData(object)[["CB"]]/100)
|
|
| 23 | 8 |
|
| 9 |
+setReplaceMethod("CA", signature(object="CopyNumberSet", value="matrix"),
|
|
| 10 |
+ function(object, value){
|
|
| 11 |
+ assayDataElementReplace(object, "CA", value) |
|
| 12 |
+ }) |
|
| 24 | 13 |
|
| 25 |
- |
|
| 26 |
-setMethod("batch", "CopyNumberSet", function(object){
|
|
| 27 |
- if("batch" %in% varLabels(object)){
|
|
| 28 |
- result <- object$batch |
|
| 29 |
- } else {
|
|
| 30 |
- stop("batch not in varLabels of CopyNumberSet")
|
|
| 31 |
- } |
|
| 32 |
- return(result) |
|
| 33 |
-}) |
|
| 14 |
+setReplaceMethod("CB", signature(object="CopyNumberSet", value="matrix"),
|
|
| 15 |
+ function(object, value){
|
|
| 16 |
+ assayDataElementReplace(object, "CB", value) |
|
| 17 |
+ }) |
|
| 34 | 18 |
|
| 35 | 19 |
setMethod("copyNumber", "CopyNumberSet", function(object){
|
| 36 |
- require(paste(annotation(object), "Crlmm", sep=""), character.only=TRUE) || stop(paste("Annotation package ", annotation(object), "Crlmm not available", sep=""))
|
|
| 37 |
- ##ensure that 2 + NA = 2 by replacing NA's with zero |
|
| 38 |
- ##the above results in copy number 0, 1, or 2 depending on the genotype....safer just to drop |
|
| 20 |
+ I <- isSnp(object) |
|
| 39 | 21 |
CA <- CA(object) |
| 40 | 22 |
CB <- CB(object) |
| 41 |
- ##nas <- is.na(CA) & is.na(CB) |
|
| 42 |
- ##CA[is.na(CA)] <- 0 |
|
| 43 |
- ##CB[is.na(CB)] <- 0 |
|
| 44 | 23 |
CN <- CA + CB |
| 45 | 24 |
##For nonpolymorphic probes, CA is the total copy number |
| 46 |
- CN[cnIndex(object, annotation(object)), ] <- CA(object)[cnIndex(object, annotation(object)), ] |
|
| 47 |
- ##if both CA and CB are NA, report NA |
|
| 48 |
- ##CN[nas] <- NA |
|
| 25 |
+ CN[!I, ] <- CA(object)[!I, ] |
|
| 49 | 26 |
CN |
| 50 | 27 |
}) |
| 51 | 28 |
|
| ... | ... |
@@ -98,3 +75,5 @@ ellipse.CopyNumberSet <- function(x, copynumber, ...){
|
| 98 | 75 |
|
| 99 | 76 |
|
| 100 | 77 |
|
| 78 |
+ |
|
| 79 |
+ |
| 101 | 80 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,298 @@ |
| 1 |
+setMethod("initialize", "CrlmmSet",
|
|
| 2 |
+ function(.Object, |
|
| 3 |
+ phenoData, |
|
| 4 |
+ featureData, |
|
| 5 |
+ annotation, |
|
| 6 |
+ experimentData, |
|
| 7 |
+ protocolData, |
|
| 8 |
+ calls=new("matrix"),
|
|
| 9 |
+ callsConfidence=new("matrix"),
|
|
| 10 |
+ senseThetaA=new("matrix"),
|
|
| 11 |
+ senseThetaB=new("matrix"),
|
|
| 12 |
+ CA=new("matrix"),
|
|
| 13 |
+ CB=new("matrix"), ... ){
|
|
| 14 |
+ ad <- assayDataNew("lockedEnvironment",
|
|
| 15 |
+ calls=calls, |
|
| 16 |
+ callsConfidence=callsConfidence, |
|
| 17 |
+ senseThetaA=senseThetaA, |
|
| 18 |
+ senseThetaB=senseThetaB, |
|
| 19 |
+ CA=CA, |
|
| 20 |
+ CB=CB) |
|
| 21 |
+ assayData(.Object) <- ad |
|
| 22 |
+ if (missing(phenoData)){
|
|
| 23 |
+ phenoData(.Object) <- annotatedDataFrameFrom(calls, byrow=FALSE) |
|
| 24 |
+ } else phenoData(.Object) <- phenoData |
|
| 25 |
+ if (missing(featureData)){
|
|
| 26 |
+ featureData(.Object) <- annotatedDataFrameFrom(calls, byrow=TRUE) |
|
| 27 |
+ } else featureData(.Object) <- featureData |
|
| 28 |
+ if(!missing(annotation)) annotation(.Object) <- annotation |
|
| 29 |
+ if(!missing(experimentData)) experimentData(.Object) <- experimentData |
|
| 30 |
+ if(!missing(protocolData)) protocolData(.Object) <- protocolData |
|
| 31 |
+ .Object |
|
| 32 |
+ }) |
|
| 33 |
+ |
|
| 34 |
+setAs("SnpCallSetPlus", "CrlmmSet",
|
|
| 35 |
+ function(from, to){
|
|
| 36 |
+ CA <- CB <- matrix(NA, nrow(from), ncol(from)) |
|
| 37 |
+ dimnames(CA) <- dimnames(CB) <- list(featureNames(from), sampleNames(from)) |
|
| 38 |
+ new("CrlmmSet",
|
|
| 39 |
+ calls=calls(from), |
|
| 40 |
+ callsConfidence=callsConfidence(from), |
|
| 41 |
+ senseThetaA=A(from), |
|
| 42 |
+ senseThetaB=B(from), |
|
| 43 |
+ CA=CA, |
|
| 44 |
+ CB=CB, |
|
| 45 |
+ phenoData=phenoData(from), |
|
| 46 |
+ experimentData=experimentData(from), |
|
| 47 |
+ annotation=annotation(from), |
|
| 48 |
+ protocolData=protocolData(from), |
|
| 49 |
+ featureData=featureData(from)) |
|
| 50 |
+ }) |
|
| 51 |
+ |
|
| 52 |
+setValidity("CrlmmSet", function(object) {
|
|
| 53 |
+ assayDataValidMembers(assayData(object), c("CA", "CB", "calls", "callsConfidence", "senseThetaA", "senseThetaB"))
|
|
| 54 |
+}) |
|
| 55 |
+ |
|
| 56 |
+ |
|
| 57 |
+ |
|
| 58 |
+setMethod("computeCopynumber", "CrlmmSet", function(object, cnOptions){
|
|
| 59 |
+ computeCopynumber.CrlmmSet(object, cnOptions) |
|
| 60 |
+}) |
|
| 61 |
+ |
|
| 62 |
+setMethod("computeCopynumber", "character", function(object, cnOptions){
|
|
| 63 |
+ crlmmFile <- object |
|
| 64 |
+ isCrlmmSet <- length(grep("crlmmSet", crlmmFile[1])) > 0
|
|
| 65 |
+ for(i in seq(along=crlmmFile)){
|
|
| 66 |
+ cat("Processing ", crlmmFile[i], "...\n")
|
|
| 67 |
+ load(crlmmFile[i]) |
|
| 68 |
+ if(isCrlmmSet){
|
|
| 69 |
+ object <- get("crlmmSet")
|
|
| 70 |
+ } else {
|
|
| 71 |
+ object <- get("callSetPlus")
|
|
| 72 |
+ } |
|
| 73 |
+ CHR <- unique(chromosome(object)) |
|
| 74 |
+ ##if(length(CHR) > 1) stop("More than one chromosome in the object. This method requires one chromosome at a time.")
|
|
| 75 |
+ if(all(CHR==24)){
|
|
| 76 |
+ message("skipping chromosome 24")
|
|
| 77 |
+ next() |
|
| 78 |
+ } |
|
| 79 |
+ cat("----------------------------------------------------------------------------\n")
|
|
| 80 |
+ cat("- Estimating copy number for chromosome", CHR, "\n")
|
|
| 81 |
+ cat("----------------------------------------------------------------------------\n")
|
|
| 82 |
+ crlmmSet <- computeCopynumber(object, cnOptions) |
|
| 83 |
+ save(crlmmSet, file=file.path(dirname(crlmmFile), paste("crlmmSet_", CHR, ".rda", sep="")))
|
|
| 84 |
+ if(!isCrlmmSet) if(cnOptions[["unlink"]]) unlink(crlmmFile[i]) |
|
| 85 |
+ rm(object, crlmmSet); gc(); |
|
| 86 |
+ } |
|
| 87 |
+}) |
|
| 88 |
+ |
|
| 89 |
+setMethod("pr", signature(object="CrlmmSet",
|
|
| 90 |
+ name="character", |
|
| 91 |
+ batch="ANY", |
|
| 92 |
+ value="numeric"), |
|
| 93 |
+ function(object, name, batch, value){
|
|
| 94 |
+ label <- paste(name, batch, sep="_") |
|
| 95 |
+ colindex <- match(label, fvarLabels(object)) |
|
| 96 |
+ if(length(colindex) == 1){
|
|
| 97 |
+ fData(object)[, colindex] <- value |
|
| 98 |
+ } |
|
| 99 |
+ if(is.na(colindex)){
|
|
| 100 |
+ stop(paste(label, " not found in object")) |
|
| 101 |
+ } |
|
| 102 |
+ if(length(colindex) > 1){
|
|
| 103 |
+ stop(paste(label, " not unique")) |
|
| 104 |
+ } |
|
| 105 |
+ object |
|
| 106 |
+ }) |
|
| 107 |
+ |
|
| 108 |
+setMethod("computeEmission", "CrlmmSet", function(object, hmmOptions){
|
|
| 109 |
+ computeEmission.CrlmmSet(object, hmmOptions) |
|
| 110 |
+}) |
|
| 111 |
+ |
|
| 112 |
+setMethod("computeEmission", "character", function(object, hmmOptions){
|
|
| 113 |
+ filename <- object |
|
| 114 |
+ chrom <- gsub(".rda", "", strsplit(filename, "_")[[1]][[2]])
|
|
| 115 |
+ if(hmmOptions[["verbose"]]) |
|
| 116 |
+ message("Compute emission probabilities for chromosome ", chrom)
|
|
| 117 |
+ if(file.exists(filename)){
|
|
| 118 |
+ load(filename) |
|
| 119 |
+ crlmmSet <- get("crlmmSet")
|
|
| 120 |
+ } else {
|
|
| 121 |
+ stop("File ", filename, " does not exist.")
|
|
| 122 |
+ } |
|
| 123 |
+ emission <- computeEmission(crlmmSet, hmmOptions) |
|
| 124 |
+ message("Saving ", file.path(dirname(filename), paste("emission_", chrom, ".rda", sep="")))
|
|
| 125 |
+ if(hmmOptions[["save.it"]]){
|
|
| 126 |
+ save(emission, |
|
| 127 |
+ file=file.path(dirname(filename), paste("emission_", chrom, ".rda", sep="")))
|
|
| 128 |
+ } |
|
| 129 |
+ return(emission) |
|
| 130 |
+}) |
|
| 131 |
+ |
|
| 132 |
+setMethod("computeHmm", "CrlmmSet", function(object, hmmOptions){
|
|
| 133 |
+ computeHmm.CrlmmSet(object, hmmOptions) |
|
| 134 |
+}) |
|
| 135 |
+ |
|
| 136 |
+##setMethod("computeHmm", "SnpCallSetPlus", function(object, hmmOptions){
|
|
| 137 |
+## crlmmSet <- computeCopynumber(object, hmmOptions) |
|
| 138 |
+## computeHmm(crlmmSet, hmmOptions) |
|
| 139 |
+##}) |
|
| 140 |
+ |
|
| 141 |
+## Genotype everything to get callSetPlus objects |
|
| 142 |
+## Go from callSets to Segments sets, writing only the segment set to file |
|
| 143 |
+## Safe, but very inefficient. Writes the quantile normalized data to file several times... |
|
| 144 |
+setMethod("computeHmm", "character", function(object, hmmOptions){
|
|
| 145 |
+ outdir <- cnOptions[["outdir"]] |
|
| 146 |
+ hmmOptions <- hmmOptions[["hmmOpts"]] |
|
| 147 |
+ filenames <- object |
|
| 148 |
+ for(i in seq(along=filenames)){
|
|
| 149 |
+ chrom <- gsub(".rda", "", strsplit(filenames[i], "_")[[1]][[2]])
|
|
| 150 |
+ if(hmmOptions[["verbose"]]) |
|
| 151 |
+ message("Fitting HMM to chromosome ", chrom)
|
|
| 152 |
+ if(file.exists(filenames[i])){
|
|
| 153 |
+ message("Loading ", filenames[i])
|
|
| 154 |
+ load(filenames[i]) |
|
| 155 |
+ crlmmSet <- get("crlmmSet")
|
|
| 156 |
+ } else {
|
|
| 157 |
+ stop("File ", filenames[i], " does not exist.")
|
|
| 158 |
+ } |
|
| 159 |
+ hmmOptions$emission <- computeEmission(filenames[i], hmmOptions) |
|
| 160 |
+ segmentSet <- computeHmm(crlmmSet, hmmOptions) |
|
| 161 |
+ ##MIN.MARKERS <- hmmOptions[["MIN.MARKERS"]] |
|
| 162 |
+ ##segmentSet <- segments[segments$nprobes >= MIN.MARKERS, ] |
|
| 163 |
+ message("Saving ", file.path(outdir, paste("segmentSet_", chrom, ".rda", sep="")))
|
|
| 164 |
+ save(segmentSet, |
|
| 165 |
+ file=file.path(outdir, paste("segmentSet_", chrom, ".rda", sep="")))
|
|
| 166 |
+ unlink(file.path(outdir, paste("crlmmSet_", chrom, ".rda", sep="")))
|
|
| 167 |
+ } |
|
| 168 |
+ fns <- list.files(outdir, pattern="segmentSet", full.names=TRUE) |
|
| 169 |
+ return(fns) |
|
| 170 |
+}) |
|
| 171 |
+ |
|
| 172 |
+##initializeCrlmmSet <- function(annotation, sns, backingfile="./"){
|
|
| 173 |
+## require(bigmemory) |
|
| 174 |
+## message("Initializing CrlmmSet file--be patient...")
|
|
| 175 |
+## path <- system.file("extdata", package=paste(annotation, "Crlmm", sep=""))
|
|
| 176 |
+## load(file.path(path, "snpProbes.rda")) |
|
| 177 |
+## load(file.path(path, "cnProbes.rda")) |
|
| 178 |
+## nr <- nrow(snpProbes)+nrow(cnProbes) |
|
| 179 |
+## pos <- c(snpProbes[, "position"], cnProbes[, "position"]) |
|
| 180 |
+## chr <- c(snpProbes[, "chrom"], cnProbes[, "chrom"]) |
|
| 181 |
+## index <- order(chr, pos) |
|
| 182 |
+## fns <- c(rownames(snpProbes), rownames(cnProbes))[index] |
|
| 183 |
+#### fns <- fns[1:20] |
|
| 184 |
+## nr <- length(fns) |
|
| 185 |
+## nc <- length(sns) |
|
| 186 |
+#### if(!file.exists(file.path(backingfile, "A.bin"))){
|
|
| 187 |
+## AA <- filebacked.big.matrix(nr, nc, |
|
| 188 |
+## type="integer", |
|
| 189 |
+## init=0, |
|
| 190 |
+## backingpath=backingfile, |
|
| 191 |
+## backingfile="A.bin", |
|
| 192 |
+## descriptorfile="A.desc", |
|
| 193 |
+## dimnames=list(fns, sns)) |
|
| 194 |
+#### } else {
|
|
| 195 |
+#### AA <- attach.big.matrix("A.desc", backingpath=backingfile)
|
|
| 196 |
+#### } |
|
| 197 |
+#### if(!file.exists(file.path(backingfile, "B.bin"))){
|
|
| 198 |
+## BB <- filebacked.big.matrix(nr, nc, |
|
| 199 |
+## type="integer", |
|
| 200 |
+## init=0, |
|
| 201 |
+## backingpath=backingfile, |
|
| 202 |
+## backingfile="B.bin", |
|
| 203 |
+## descriptorfile="B.desc", |
|
| 204 |
+## dimnames=list(fns, sns)) |
|
| 205 |
+#### } else {
|
|
| 206 |
+#### BB <- attach.big.matrix("B.desc", backingpath=backingfile)
|
|
| 207 |
+#### } |
|
| 208 |
+#### if(!file.exists(file.path(backingfile, "GT.bin"))){
|
|
| 209 |
+## gt <- filebacked.big.matrix(nr, nc, |
|
| 210 |
+## type="integer", |
|
| 211 |
+## init=0, |
|
| 212 |
+## backingpath=backingfile, |
|
| 213 |
+## backingfile="GT.bin", |
|
| 214 |
+## descriptorfile="GT.desc", |
|
| 215 |
+## dimnames=list(fns, sns)) |
|
| 216 |
+## confs <- filebacked.big.matrix(nr, nc, |
|
| 217 |
+## type="integer", |
|
| 218 |
+## init=0, |
|
| 219 |
+## backingpath=backingfile, |
|
| 220 |
+## backingfile="confs.bin", |
|
| 221 |
+## descriptorfile="confs.desc", |
|
| 222 |
+## dimnames=list(fns, sns)) |
|
| 223 |
+#### } else {
|
|
| 224 |
+#### gt <- attach.big.matrix("GT.desc", backingpath=backingfile)
|
|
| 225 |
+#### } |
|
| 226 |
+#### if(!file.exists(file.path(backingfile, "CA.bin"))){
|
|
| 227 |
+## ca <- filebacked.big.matrix(nr, nc, |
|
| 228 |
+## type="integer", |
|
| 229 |
+## init=0, |
|
| 230 |
+## backingpath=backingfile, |
|
| 231 |
+## backingfile="CA.bin", |
|
| 232 |
+## descriptorfile="CA.desc", |
|
| 233 |
+## dimnames=list(fns, sns)) |
|
| 234 |
+#### } else {
|
|
| 235 |
+#### ca <- attach.big.matrix("CA.desc", backingpath=backingfile)
|
|
| 236 |
+#### } |
|
| 237 |
+#### if(!file.exists(file.path(backingfile, "CB.bin"))){
|
|
| 238 |
+## cb <- filebacked.big.matrix(nr, nc, |
|
| 239 |
+## type="integer", |
|
| 240 |
+## init=0, |
|
| 241 |
+## backingpath=backingfile, |
|
| 242 |
+## backingfile="CB.bin", |
|
| 243 |
+## descriptorfile="CB.desc", |
|
| 244 |
+## dimnames=list(fns, sns)) |
|
| 245 |
+#### } else {
|
|
| 246 |
+#### cb <- attach.big.matrix("CB.desc", backingpath=backingfile)
|
|
| 247 |
+#### } |
|
| 248 |
+## return(new("CrlmmSet", A=AA, B=BB, CA=ca, CB=cb, GT=gt, confs=confs))
|
|
| 249 |
+##} |
|
| 250 |
+ |
|
| 251 |
+ |
|
| 252 |
+##setMethod("A", signature(object="CrlmmSet"),
|
|
| 253 |
+## function(object) assayDataElement(object,"A")) |
|
| 254 |
+##setMethod("B", signature(object="CrlmmSet"),
|
|
| 255 |
+## function(object) assayDataElement(object,"A")) |
|
| 256 |
+##setMethod("CA", signature(object="CrlmmSet"),
|
|
| 257 |
+## function(object) assayDataElement(object,"A")) |
|
| 258 |
+##setMethod("CB", signature(object="CrlmmSet"),
|
|
| 259 |
+## function(object) assayDataElement(object,"A")) |
|
| 260 |
+##setMethod("GT", signature(object="CrlmmSet"),
|
|
| 261 |
+## function(object) assayDataElement(object,"GT")) |
|
| 262 |
+ |
|
| 263 |
+##setReplaceMethod("CA", signature(object="CrlmmSetList", value="big.matrix"),
|
|
| 264 |
+## function(object, value){
|
|
| 265 |
+## assayDataElementReplace(object, "A", value) |
|
| 266 |
+## }) |
|
| 267 |
+##setReplaceMethod("CB", signature(object="CrlmmSetList", value="big.matrix"),
|
|
| 268 |
+## function(object, value){
|
|
| 269 |
+## CB(object[[3]]) <- value |
|
| 270 |
+## object |
|
| 271 |
+## }) |
|
| 272 |
+## |
|
| 273 |
+##setReplaceMethod("A", signature(object="CrlmmSet", value="big.matrix"),
|
|
| 274 |
+## function(object, value){
|
|
| 275 |
+## assayDataElementReplace(object, "A", value) |
|
| 276 |
+## }) |
|
| 277 |
+ |
|
| 278 |
+##setReplaceMethod("GT", signature(object="CrlmmSet", value="big.matrix"),
|
|
| 279 |
+## function(object, value){
|
|
| 280 |
+## assayDataElementReplace(object, "GT", value) |
|
| 281 |
+## }) |
|
| 282 |
+##setReplaceMethod("confs", signature(object="CrlmmSet", value="big.matrix"),
|
|
| 283 |
+## function(object, value){
|
|
| 284 |
+## assayDataElementReplace(object, "callProbability", value) |
|
| 285 |
+## }) |
|
| 286 |
+##setReplaceMethod("confs", signature(object="CrlmmSet", value="matrix"),
|
|
| 287 |
+## function(object, value){
|
|
| 288 |
+## assayDataElementReplace(object, "callProbability", value) |
|
| 289 |
+## }) |
|
| 290 |
+##setReplaceMethod("GT", signature(object="CrlmmSet", value="matrix"),
|
|
| 291 |
+## function(object, value){
|
|
| 292 |
+## assayDataElementReplace(object, "GT", value) |
|
| 293 |
+## }) |
|
| 294 |
+##setReplaceMethod("B", signature(object="CrlmmSetList", value="big.matrix"),
|
|
| 295 |
+## function(object, value){
|
|
| 296 |
+## B(object[[1]]) <- value |
|
| 297 |
+## object |
|
| 298 |
+## }) |
| 0 | 299 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,39 @@ |
| 1 |
+##setAs("SnpCallSetPlus", "CrlmmSetFF",
|
|
| 2 |
+## function(from, to){
|
|
| 3 |
+## CA <- CB <- matrix(NA, nrow(from), ncol(from)) |
|
| 4 |
+## dimnames(CA) <- dimnames(CB) <- list(featureNames(from), sampleNames(from)) |
|
| 5 |
+## new("CrlmmSetFF",
|
|
| 6 |
+## calls=calls(from), |
|
| 7 |
+## callsConfidence=callsConfidence(from), |
|
| 8 |
+## senseThetaA=A(from), |
|
| 9 |
+## senseThetaB=B(from), |
|
| 10 |
+## CA=CA, |
|
| 11 |
+## CB=CB, |
|
| 12 |
+## phenoData=phenoData(from), |
|
| 13 |
+## experimentData=experimentData(from), |
|
| 14 |
+## annotation=annotation(from), |
|
| 15 |
+## protocolData=protocolData(from), |
|
| 16 |
+## featureData=featureData(from)) |
|
| 17 |
+## }) |
|
| 18 |
+ |
|
| 19 |
+##setMethod("[", "CrlmmSetFF", function(x, i, j, ..., drop = FALSE) {
|
|
| 20 |
+ ## does NOT subset any of the assayDataElements |
|
| 21 |
+ ## does subset the featureData, phenoData, protocolData |
|
| 22 |
+ |
|
| 23 |
+ ## accessors to the assayData work as follows |
|
| 24 |
+ ## object2 <- object[1:10, ] |
|
| 25 |
+ ## |
|
| 26 |
+ ## A(object2)[1:10, ] |
|
| 27 |
+ ## - i. A(object2) gets the file |
|
| 28 |
+ ## - ii. A is subset by the dimNames information as follows |
|
| 29 |
+ ## dimNames[[1]] %in% featureNames(object2) |
|
| 30 |
+ ## dimNames[[2]] %in% sampleNames(object2) |
|
| 31 |
+ ## iii. first 10 elements are retrieved |
|
| 32 |
+##}) |
|
| 33 |
+ |
|
| 34 |
+##setMethod("A", "CrlmmSetFF", function(object) assayData(object)[["senseThetaA"]])
|
|
| 35 |
+##setMethod("B", "CrlmmSetFF", function(object) assayData(object)[["senseThetaB"]])
|
|
| 36 |
+##setMethod("CA", "CrlmmSetFF", function(object) assayData(object)[["CA"]]/100)
|
|
| 37 |
+##setMethod("CB", "CrlmmSetFF", function(object) assayData(object)[["CB"]]/100)
|
|
| 38 |
+ |
|
| 39 |
+ |
| 0 | 40 |
deleted file mode 100644 |
| ... | ... |
@@ -1,317 +0,0 @@ |
| 1 |
-setMethod("[", "CrlmmSetList", function(x, i, j, ..., drop = FALSE){
|
|
| 2 |
- if (missing(drop)) drop <- FALSE |
|
| 3 |
- if (missing(i) && missing(j)) |
|
| 4 |
- {
|
|
| 5 |
- if (length(list(...))!=0) |
|
| 6 |
- stop("specify genes or samples to subset; use '",
|
|
| 7 |
- substitute(x), "$", names(list(...))[[1]], |
|
| 8 |
- "' to access phenoData variables") |
|
| 9 |
- return(x) |
|
| 10 |
- } |
|
| 11 |
- if (!missing(j)){
|
|
| 12 |
- f1 <- function(x, j){
|
|
| 13 |
- x <- x[, j] |
|
| 14 |
- } |
|
| 15 |
- x <- lapply(x, f1, j) |
|
| 16 |
- } |
|
| 17 |
- if(!missing(i)){
|
|
| 18 |
- f2 <- function(x, i){
|
|
| 19 |
- x <- x[i, ] |
|
| 20 |
- } |
|
| 21 |
- x <- lapply(x, f2, i) |
|
| 22 |
- } |
|
| 23 |
- as(x, "CrlmmSetList") |
|
| 24 |
-}) |
|
| 25 |
- |
|
| 26 |
-setMethod("$", "CrlmmSetList", function(x, name) {
|
|
| 27 |
- ##if(!(name %in% .parameterNames()[output(x) != 0])){
|
|
| 28 |
- if(length(x) != 3){
|
|
| 29 |
- stop("'$' operature reserved for accessing parameter names in CopyNumberSet object. CrlmmSetList must be of length 3")
|
|
| 30 |
- } |
|
| 31 |
- j <- grep(name, fvarLabels(x[[3]])) |
|
| 32 |
- if(length(j) < 1) |
|
| 33 |
- stop(name, " not in fvarLabels of CopyNumberSet object") |
|
| 34 |
- if(length(j) > 1){
|
|
| 35 |
- warning("Multiple instances of ", name, " in fvarLabels. Using the first instance")
|
|
| 36 |
- j <- j[1] |
|
| 37 |
- } |
|
| 38 |
- param <- fData(x[[3]])[, j] |
|
| 39 |
- param |
|
| 40 |
-}) |
|
| 41 |
-setMethod(".harmonizeDimnames", "CrlmmSetList", function(object){
|
|
| 42 |
- i <- length(object) |
|
| 43 |
- while(i > 1){
|
|
| 44 |
- object[[i-1]] <- harmonizeDimnamesTo(object[[i-1]], object[[i]]) |
|
| 45 |
- i <- i-1 |
|
| 46 |
- } |
|
| 47 |
- object |
|
| 48 |
-}) |
|
| 49 |
- |
|
| 50 |
-setMethod("A", "CrlmmSetList", function(object) A(object[[1]]))
|
|
| 51 |
- |
|
| 52 |
-setMethod("addFeatureAnnotation", "CrlmmSetList", function(object, ...){
|
|
| 53 |
- ##if(missing(CHR)) stop("Must specificy chromosome")
|
|
| 54 |
- cdfName <- annotation(object) |
|
| 55 |
- pkgname <- paste(cdfName, "Crlmm", sep="") |
|
| 56 |
- path <- system.file("extdata", package=pkgname)
|
|
| 57 |
- loader("cnProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
|
|
| 58 |
- cnProbes <- get("cnProbes", envir=.crlmmPkgEnv)
|
|
| 59 |
- loader("snpProbes.rda", pkgname=pkgname, envir=.crlmmPkgEnv)
|
|
| 60 |
- snpProbes <- get("snpProbes", envir=.crlmmPkgEnv)
|
|
| 61 |
- |
|
| 62 |
- ##Feature Data |
|
| 63 |
- snps <- featureNames(object)[snpIndex(object)] |
|
| 64 |
- nps <- featureNames(object)[cnIndex(object)] |
|
| 65 |
- position.snp <- snpProbes[match(snps, rownames(snpProbes)), "position"] |
|
| 66 |
- names(position.snp) <- snps |
|
| 67 |
- position.np <- cnProbes[match(nps, rownames(cnProbes)), "position"] |
|
| 68 |
- names(position.np) <- nps |
|
| 69 |
- |
|
| 70 |
- J <- grep("chr", colnames(snpProbes))
|
|
| 71 |
- chr.snp <- snpProbes[match(snps, rownames(snpProbes)), J] |
|
| 72 |
- chr.np <- cnProbes[match(nps, rownames(cnProbes)), J] |
|
| 73 |
- |
|
| 74 |
- position <- c(position.snp, position.np) |
|
| 75 |
- chrom <- c(chr.snp, chr.np) |
|
| 76 |
- |
|
| 77 |
- ##We may not have annotation for all of the snps |
|
| 78 |
- if(!all(featureNames(object) %in% names(position))){
|
|
| 79 |
- message("Dropping loci for which physical position is not available.")
|
|
| 80 |
- object <- object[featureNames(object) %in% names(position), ] |
|
| 81 |
- } |
|
| 82 |
- ix <- match(featureNames(object), names(position)) |
|
| 83 |
- position <- position[ix] |
|
| 84 |
- chrom <- chrom[ix] |
|
| 85 |
- ##require(SNPchip) |
|
| 86 |
- chrom <- chromosome2integer(chrom) |
|
| 87 |
- |
|
| 88 |
- stopifnot(identical(names(position), featureNames(object))) |
|
| 89 |
- if(sum(duplicated(names(position))) > 0){
|
|
| 90 |
- warning("Removing rows with NA identifiers...")
|
|
| 91 |
- ##RS: fix this |
|
| 92 |
- I <- which(!is.na(names(position))) |
|
| 93 |
- } else I <- seq(along=names(position)) |
|
| 94 |
- fd <- data.frame(cbind(chrom[I], |
|
| 95 |
- position[I])) |
|
| 96 |
- colnames(fd) <- c("chromosome", "position")
|
|
| 97 |
- rownames(fd) <- featureNames(object) |
|
| 98 |
- fD <- new("AnnotatedDataFrame",
|
|
| 99 |
- data=fd, |
|
| 100 |
- varMetadata=data.frame(labelDescription=colnames(fd))) |
|
| 101 |
- return(fD) |
|
| 102 |
-}) |
|
| 103 |
- |
|
| 104 |
-setMethod("annotation", "CrlmmSetList", function(object) annotation(object[[1]]))
|
|
| 105 |
-setMethod("B", "CrlmmSetList", function(object) B(object[[1]]))
|
|
| 106 |
-setMethod("batch", "CrlmmSetList", function(object) batch(object[[3]]))
|
|
| 107 |
-setMethod("CA", "CrlmmSetList", function(object, ...) CA(object[[3]], ...))
|
|
| 108 |
-setMethod("CB", "CrlmmSetList", function(object, ...) CB(object[[3]], ...))
|
|
| 109 |
-setMethod("calls", "CrlmmSetList", function(object) calls(object[[2]]))
|
|
| 110 |
-setMethod("chromosome", "CrlmmSetList", function(object){
|
|
| 111 |
- chr <- NULL |
|
| 112 |
- for(i in 1:length(object)){
|
|
| 113 |
- if(length(fvarLabels(object[[i]])) > 0){
|
|
| 114 |
- if("chromosome" %in% fvarLabels(object[[i]])){
|
|
| 115 |
- chr <- chromosome(object[[i]]) |
|
| 116 |
- break() |
|
| 117 |
- } |
|
| 118 |
- } |
|
| 119 |
- } |
|
| 120 |
- if(is.null(chr)) warning("fvarLabel 'chromosome' not in any element of the CrlmmSetList object")
|
|
| 121 |
- return(chr) |
|
| 122 |
- ##chromosome(object[[3]]) |
|
| 123 |
- }) |
|
| 124 |
- |
|
| 125 |
-setMethod("position", "CrlmmSetList", function(object){
|
|
| 126 |
- pos <- NULL |
|
| 127 |
- for(i in 1:length(object)){
|
|
| 128 |
- if(length(fvarLabels(object[[i]])) > 0){
|
|
| 129 |
- if("position" %in% fvarLabels(object[[i]])){
|
|
| 130 |
- pos <- position(object[[i]]) |
|
| 131 |
- break() |
|
| 132 |
- } |
|
| 133 |
- } else next() |
|
| 134 |
- } |
|
| 135 |
- if(is.null(pos)) warning("fvarLabel 'position' not in any element of the CrlmmSetList object")
|
|
| 136 |
- return(pos) |
|
| 137 |
- }) |
|
| 138 |
- |
|
| 139 |
-setMethod("cnIndex", "CrlmmSetList", function(object, ...) {
|
|
| 140 |
- match(cnNames(object[[1]], annotation(object)), featureNames(object)) |
|
| 141 |
-}) |
|
| 142 |
-setMethod("combine", signature=signature(x="CrlmmSetList", y="CrlmmSetList"),
|
|
| 143 |
- function(x, y, ...){
|
|
| 144 |
- x.abset <- x[[1]] |
|
| 145 |
- y.abset <- y[[1]] |
|
| 146 |
- x.snpset <- x[[2]] |
|
| 147 |
- y.snpset <- y[[2]] |
|
| 148 |
- abset <- combine(x.abset, y.abset) |
|
| 149 |
- ##we have hijacked the featureData slot to store parameters. Biobase will not allow combining our 'feature' data. |
|
| 150 |
- warning("the featureData is not easily combined... removing the featureData")
|
|
| 151 |
- ##fd1 <- featureData(x.snpset) |
|
| 152 |
- ##fd2 <- featureData(y.snpset) |
|
| 153 |
- featureData(x.snpset) <- annotatedDataFrameFrom(calls(x.snpset), byrow=TRUE) |
|
| 154 |
- featureData(y.snpset) <- annotatedDataFrameFrom(calls(y.snpset), byrow=TRUE) |
|
| 155 |
- snpset <- combine(x.snpset, y.snpset) |
|
| 156 |
- merged <- list(abset, snpset) |
|
| 157 |
- merged <- as(merged, "CrlmmSetList") |
|
| 158 |
- merged |
|
| 159 |
- }) |
|
| 160 |
-setMethod("confs", "CrlmmSetList", function(object) confs(object[[2]]))
|
|
| 161 |
-setMethod("copyNumber", "CrlmmSetList", function(object) copyNumber(object[[3]]))
|
|
| 162 |
-setMethod("dims", "CrlmmSetList", function(object) sapply(object, dim))
|
|
| 163 |
-setMethod("featureNames", "CrlmmSetList", function(object) featureNames(object[[1]]))
|
|
| 164 |
-setMethod("ncol", signature(x="CrlmmSetList"), function(x) ncol(x[[1]]))
|
|
| 165 |
-setMethod("nrow", signature(x="CrlmmSetList"), function(x) nrow(x[[1]]))
|
|
| 166 |
-setMethod("plot", signature(x="CrlmmSetList"),
|
|
| 167 |
- function(x, y, ...){
|
|
| 168 |
- A <- log2(A(x)) |
|
| 169 |
- B <- log2(B(x)) |
|
| 170 |
- plot(A, B, ...) |
|
| 171 |
- }) |
|
| 172 |
-setMethod("points", signature(x="CrlmmSetList"),
|
|
| 173 |
- function(x, y, ...){
|
|
| 174 |
- A <- log2(A(x)) |
|
| 175 |
- B <- log2(B(x)) |
|
| 176 |
- points(A, B, ...) |
|
| 177 |
- }) |
|
| 178 |
-setMethod("sampleNames", "CrlmmSetList", function(object) sampleNames(object[[1]]))
|
|
| 179 |
-setMethod("show", "CrlmmSetList", function(object){
|
|
| 180 |
- cat("\n Elements in CrlmmSetList object: \n")
|
|
| 181 |
- cat("\n")
|
|
| 182 |
- for(i in 1:length(object)){
|
|
| 183 |
- cat("class: ", class(object[[i]]), "\n")
|
|
| 184 |
- cat("assayData elements: ", ls(assayData(object[[i]])), "\n")
|
|
| 185 |
- cat("Dimensions: ", dim(object[[i]]))
|
|
| 186 |
- cat("\n \n")
|
|
| 187 |
- } |
|
| 188 |
-}) |
|
| 189 |
-setMethod("snpIndex", "CrlmmSetList", function(object, ...){
|
|
| 190 |
- match(snpNames(object[[1]], annotation(object)), featureNames(object)) |
|
| 191 |
-}) |
|
| 192 |
-setMethod("splitByChromosome", "CrlmmSetList", function(object, cdfName, outdir){
|
|
| 193 |
- path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))
|
|
| 194 |
- load(file.path(path, "snpProbes.rda")) |
|
| 195 |
- load(file.path(path, "cnProbes.rda")) |
|
| 196 |
- k <- grep("chr", colnames(snpProbes))
|
|
| 197 |
- if(length(k) < 1) stop("chr or chromosome not in colnames(snpProbes)")
|
|
| 198 |
- for(CHR in 1:24){
|
|
| 199 |
- cat("Chromosome ", CHR, "\n")
|
|
| 200 |
- snps <- rownames(snpProbes)[snpProbes[, k] == CHR] |
|
| 201 |
- cnps <- rownames(cnProbes)[cnProbes[, k] == CHR] |
|
| 202 |
- index <- c(match(snps, featureNames(object)), |
|
| 203 |
- match(cnps, featureNames(object))) |
|
| 204 |
- index <- index[!is.na(index)] |
|
| 205 |
- crlmmSetList <- object[index, ] |
|
| 206 |
- save(crlmmSetList, file=file.path(outdir, paste("crlmmSetList_", CHR, ".rda", sep="")))
|
|
| 207 |
- } |
|
| 208 |
-}) |
|
| 209 |
- |
|
| 210 |
-setMethod("update", "CrlmmSetList", function(object, ...){
|
|
| 211 |
- if(length(crlmmSetList) == 3){
|
|
| 212 |
- message("copy number object already present. Nothing to do.")
|
|
| 213 |
- return() |
|
| 214 |
- } |
|
| 215 |
- CHR <- unique(chromosome(crlmmSetList[[1]])) |
|
| 216 |
- if(length(CHR) > 1) stop("More than one chromosome in the object. This method requires one chromosome at a time.")
|
|
| 217 |
- if(CHR == 24){
|
|
| 218 |
- message("A solution for chromosome 24 is not yet available.")
|
|
| 219 |
- return() |
|
| 220 |
- } |
|
| 221 |
- computeCopynumber(object, CHR=CHR, ...) |
|
| 222 |
-}) |
|
| 223 |
- |
|
| 224 |
-setReplaceMethod("CA", signature(object="CrlmmSetList", value="matrix"),
|
|
| 225 |
- function(object, value){
|
|
| 226 |
- CA(object[[3]]) <- value |
|
| 227 |
- object |
|
| 228 |
- }) |
|
| 229 |
-setReplaceMethod("CB", signature(object="CrlmmSetList", value="matrix"),
|
|
| 230 |
- function(object, value){
|
|
| 231 |
- CB(object[[3]]) <- value |
|
| 232 |
- object |
|
| 233 |
- }) |
|
| 234 |
- |
|
| 235 |
-setReplaceMethod("A", signature(object="CrlmmSetList", value="matrix"),
|
|
| 236 |
- function(object, value){
|
|
| 237 |
- A(object[[1]]) <- value |
|
| 238 |
- object |
|
| 239 |
- }) |
|
| 240 |
-setReplaceMethod("B", signature(object="CrlmmSetList", value="matrix"),
|
|
| 241 |
- function(object, value){
|
|
| 242 |
- B(object[[1]]) <- value |
|
| 243 |
- object |
|
| 244 |
- }) |
|
| 245 |
- |
|
| 246 |
- |
|
| 247 |
- |
|
| 248 |
-setMethod("boxplot", "CrlmmSetList", function(x, ...){
|
|
| 249 |
-##boxplot.CrlmmSetList <- function(x, ...){
|
|
| 250 |
- if(length(x) != 3) stop("elements of list should be of class ABset, SnpSet, and CopyNumberSet, respectively.")
|
|
| 251 |
- genotypes <- calls(x)-1 |
|
| 252 |
- A1 <- A(x) |
|
| 253 |
- B1 <- B(x) |
|
| 254 |
- Alist <- split(A1, genotypes) |
|
| 255 |
- Alist <- as(rev(Alist), "data.frame") |
|
| 256 |
- Blist <- split(B1, genotypes) |
|
| 257 |
- ylim <- range(unlist(Alist)) |
|
| 258 |
- boxplot(Alist, xaxt="n", ylab=expression(I[A]), |
|
| 259 |
- cex.axis=0.6, |
|
| 260 |
- xlab="", |
|
| 261 |
- ylim=range(unlist(Alist), na.rm=TRUE), |
|
| 262 |
- border="grey50", xaxs="i", at=0:2, xlim=c(-0.5, 2.5), |
|
| 263 |
- cex.main=0.9, xaxt="n", |
|
| 264 |
- yaxt="n", |
|
| 265 |
- col=cols) |
|
| 266 |
- axis(2, at=pretty(ylim), cex.axis=0.8) |
|
| 267 |
- axis(1, at=0:2, labels=rev(c("2A (AA genotype)", "1A (AB genotype)", "0A (BB genotype)")), cex.axis=0.7)
|
|
| 268 |
- ##extracts nuA for first batch |
|
| 269 |
- suppressWarnings(nuA <- x$nuA) |
|
| 270 |
- segments(-1, nuA, |
|
| 271 |
- 0, nuA, lty=2, col="blue") |
|
| 272 |
- suppressWarnings(phiA <- x$phiA) |
|
| 273 |
- ##phiA <- fData(x)[["phiA_A"]] |
|
| 274 |
- segments(0, nuA, |
|
| 275 |
- 2.5, nuA+2.5*phiA, lwd=2, col="blue") |
|
| 276 |
- axis(2, at=nuA, labels=expression(hat(nu[A])), cex.axis=0.9) |
|
| 277 |
- text(0, ylim[1], labels=paste("n =", length(Alist[[1]])), cex=0.8)
|
|
| 278 |
- text(1, ylim[1], labels=paste("n =", length(Alist[[2]])), cex=0.8)
|
|
| 279 |
- text(2, ylim[1], labels=paste("n =", length(Alist[[3]])), cex=0.8)
|
|
| 280 |
-## segments(0.5, nuA+0.5*phiA, |
|
| 281 |
-## 1, pretty(ylim, n=5)[2], lty=2, col="blue") |
|
| 282 |
-## segments(1, pretty(ylim, n=5)[2], |
|
| 283 |
-## 1.2, pretty(ylim)[2], lty=2, col="blue") |
|
| 284 |
-## text(1.25, pretty(ylim)[2], pos=4, |
|
| 285 |
-## labels=expression(hat(nu[A]) + c[A]*hat(phi[A]))) |
|
| 286 |
- legend("topleft", fill=rev(cols), legend=c("AA", "AB", "BB"))##, title="diallelic genotypes")
|
|
| 287 |
- |
|
| 288 |
- ylim <- range(unlist(Blist)) |
|
| 289 |
- boxplot(Blist, xaxt="n", ylab=expression(I[B]), |
|
| 290 |
- ##xlab="diallelic genotypes", |
|
| 291 |
- cex.axis=0.6, |
|
| 292 |
- ylim=range(unlist(Blist), na.rm=TRUE), |
|
| 293 |
- border="grey50", xaxs="i", |
|
| 294 |
- at=0:2, xlim=c(-0.5, 2.5), |
|
| 295 |
- cex.main=0.9, yaxt="n", |
|
| 296 |
- col=rev(cols)) |
|
| 297 |
- axis(2, at=pretty(ylim), cex.axis=0.8) |
|
| 298 |
- axis(1, at=0:2, labels=c("0B (AA genotype)", "1B (AB genotype)", "2B (BB genotype)"), cex.axis=0.7)
|
|
| 299 |
- ##nuB <- fData(x)[["nuB_A"]] |
|
| 300 |
- suppressWarnings(nuB <- x$nuB) |
|
| 301 |
- segments(-1, nuB, |
|
| 302 |
- 0, nuB, lty=2, col="blue") |
|
| 303 |
- ##phiB <- fData(x[i,])[["phiB_A"]] |
|
| 304 |
- suppressWarnings(phiB <- x$phiB) |
|
| 305 |
- segments(0, nuB, |
|
| 306 |
- 2.5, nuB+2.5*phiB, lwd=2, col="blue") |
|
| 307 |
- axis(2, at=nuB, labels=expression(hat(nu[B])), cex.axis=0.9) |
|
| 308 |
- text(0, ylim[1], labels=paste("n =", length(Blist[[1]])), cex=0.8)
|
|
| 309 |
- text(1, ylim[1], labels=paste("n =", length(Blist[[2]])), cex=0.8)
|
|
| 310 |
- text(2, ylim[1], labels=paste("n =", length(Blist[[3]])), cex=0.8)
|
|
| 311 |
-## segments(0.5, nuB+0.5*phiB, |
|
| 312 |
-## 1, pretty(ylim,n=5)[2], lty=2, col="blue") |
|
| 313 |
-## segments(1, pretty(ylim, n=5)[2], |
|
| 314 |
-## 1.2, pretty(ylim,n=5)[2], lty=2, col="blue") |
|
| 315 |
-## text(1.25, pretty(ylim,n=5)[2], pos=4, labels=expression(hat(nu[B]) + c[B]*hat(phi[B]))) |
|
| 316 |
- mtext(featureNames(x)[i], 3, outer=TRUE, line=1) |
|
| 317 |
-}) |
| 318 | 0 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,109 @@ |
| 1 |
+##How to make the initialization platform-specific? |
|
| 2 |
+setMethod("initialize", "SnpCallSetPlus",
|
|
| 3 |
+ function(.Object, |
|
| 4 |
+ phenoData, featureData, |
|
| 5 |
+ calls=new("matrix"),
|
|
| 6 |
+ callsConfidence=new("matrix"),
|
|
| 7 |
+ senseThetaA=new("matrix"),
|
|
| 8 |
+ senseThetaB=new("matrix"),
|
|
| 9 |
+ annotation, |
|
| 10 |
+ experimentData, |
|
| 11 |
+ protocolData, ... ){
|
|
| 12 |
+ ad <- assayDataNew("lockedEnvironment",
|
|
| 13 |
+ calls=calls, |
|
| 14 |
+ callsConfidence=callsConfidence, |
|
| 15 |
+ senseThetaA=senseThetaA, |
|
| 16 |
+ senseThetaB=senseThetaB) |
|
| 17 |
+ assayData(.Object) <- ad |
|
| 18 |
+ if (missing(phenoData)){
|
|
| 19 |
+ phenoData(.Object) <- annotatedDataFrameFrom(calls, byrow=FALSE) |
|
| 20 |
+ } else{
|
|
| 21 |
+ phenoData(.Object) <- phenoData |
|
| 22 |
+ } |
|
| 23 |
+ if (missing(featureData)){
|
|
| 24 |
+ featureData(.Object) <- annotatedDataFrameFrom(calls, byrow=TRUE) |
|
| 25 |
+ } else{
|
|
| 26 |
+ featureData(.Object) <- featureData |
|
| 27 |
+ } |
|
| 28 |
+ if(!missing(annotation)) annotation(.Object) <- annotation |
|
| 29 |
+ if(!missing(experimentData)) experimentData(.Object) <- experimentData |
|
| 30 |
+ if(!missing(protocolData)) protocolData(.Object) <- protocolData |
|
| 31 |
+ .Object |
|
| 32 |
+ }) |
|
| 33 |
+getParam.SnpCallSetPlus <- function(object, name, batch){
|
|
| 34 |
+ label <- paste(name, batch, sep="_") |
|
| 35 |
+ colindex <- grep(label, fvarLabels(object)) |
|
| 36 |
+ if(length(colindex) == 1){
|
|
| 37 |
+ param <- fData(object)[, colindex] |
|
| 38 |
+ } |
|
| 39 |
+ if(length(colindex) < 1){
|
|
| 40 |
+ param <- NULL |
|
| 41 |
+ } |
|
| 42 |
+ if(is.na(colindex)){
|
|
| 43 |
+ stop(paste(label, " not found in object")) |
|
| 44 |
+ } |
|
| 45 |
+ if(length(colindex) > 1){
|
|
| 46 |
+ stop(paste(label, " not unique")) |
|
| 47 |
+ } |
|
| 48 |
+ return(param) |
|
| 49 |
+ } |
|
| 50 |
+ |
|
| 51 |
+setMethod("getParam", signature(object="SnpCallSetPlus",
|
|
| 52 |
+ name="character", |
|
| 53 |
+ batch="ANY"), |
|
| 54 |
+ function(object, name, batch){
|
|
| 55 |
+ if(length(batch) > 1){
|
|
| 56 |
+ warning("batch argument to getParam should have length 1. Only using the first")
|
|
| 57 |
+ batch <- batch[1] |
|
| 58 |
+ } |
|
| 59 |
+ getParam.SnpCallSetPlus(object, name, batch) |
|
| 60 |
+}) |
|
| 61 |
+ |
|
| 62 |
+setMethod("splitByChromosome", "SnpCallSetPlus", function(object, cnOptions){
|
|
| 63 |
+ tmpdir <- cnOptions[["tmpdir"]] |
|
| 64 |
+ outdir <- cnOptions[["outdir"]] |
|
| 65 |
+ save.it <- cnOptions[["save.it"]] |
|
| 66 |
+ path <- system.file("extdata", package=paste(annotation(object), "Crlmm", sep=""))
|
|
| 67 |
+ load(file.path(path, "snpProbes.rda")) |
|
| 68 |
+ snpProbes <- get("snpProbes")
|
|
| 69 |
+ load(file.path(path, "cnProbes.rda")) |
|
| 70 |
+ cnProbes <- get("cnProbes")
|
|
| 71 |
+ k <- grep("chr", colnames(snpProbes))
|
|
| 72 |
+ if(length(k) < 1) stop("chr or chromosome not in colnames(snpProbes)")
|
|
| 73 |
+ for(CHR in 1:24){
|
|
| 74 |
+ cat("Chromosome ", CHR, "\n")
|
|
| 75 |
+ snps <- rownames(snpProbes)[snpProbes[, k] == CHR] |
|
| 76 |
+ cnps <- rownames(cnProbes)[cnProbes[, k] == CHR] |
|
| 77 |
+ index <- c(match(snps, featureNames(object)), |
|
| 78 |
+ match(cnps, featureNames(object))) |
|
| 79 |
+ index <- index[!is.na(index)] |
|
| 80 |
+ callSetPlus <- object[index, ] |
|
| 81 |
+ if(CHR != 24){
|
|
| 82 |
+ crlmmSet <- computeCopynumber(callSetPlus, cnOptions) |
|
| 83 |
+ |
|
| 84 |
+ } else{
|
|
| 85 |
+ message("Copy number estimates not available for chromosome Y. Saving only the 'callSetPlus' object for this chromosome")
|
|
| 86 |
+ save(callSetPlus, file=file.path(outdir, paste("callSetPlus_", CHR, ".rda", sep="")))
|
|
| 87 |
+ } |
|
| 88 |
+ if(cnOptions[["hiddenMarkovModel"]] & CHR != 24){
|
|
| 89 |
+ segmentSet <- computeHmm(crlmmSet, cnOptions) |
|
| 90 |
+ save(segmentSet, file=file.path(outdir, paste("segmentSet_", CHR, ".rda", sep="")))
|
|
| 91 |
+ saved.objects <- list.files(outdir, pattern="segmentSet", full.names=TRUE) |
|
| 92 |
+ } else{ ## save crlmmSet to outdir
|
|
| 93 |
+ save(crlmmSet, file=file.path(outdir, paste("crlmmSet_", CHR, ".rda", sep="")))
|
|
| 94 |
+ saved.objects <- list.files(outdir, pattern="crlmmSet", full.names=TRUE) |
|
| 95 |
+ } |
|
| 96 |
+ } |
|
| 97 |
+ saved.objects |
|
| 98 |
+}) |
|
| 99 |
+setMethod("addFeatureAnnotation", "SnpCallSetPlus", function(object, ...){
|
|
| 100 |
+ addFeatureAnnotation.SnpCallSetPlus(object, ...) |
|
| 101 |
+}) |
|
| 102 |
+setMethod("computeCopynumber", "SnpCallSetPlus",
|
|
| 103 |
+ function(object, cnOptions){
|
|
| 104 |
+ computeCopynumber.SnpCallSetPlus(object, cnOptions) |
|
| 105 |
+ }) |
|
| 106 |
+setMethod("chromosome", "SnpCallSetPlus", function(object) fData(object)$chromosome)
|
|
| 107 |
+setMethod("position", "SnpCallSetPlus", function(object) fData(object)$position)
|
|
| 108 |
+gtConfidence <- function(object) 1-exp(-callsConfidence(object)/1000) |
|
| 109 |
+ |
| 0 | 4 |
new file mode 100644 |
| ... | ... |
@@ -0,0 +1,44 @@ |
| 1 |
+##setValidity("SnpQSet", function(object) {
|
|
| 2 |
+## ##msg <- validMsg(NULL, Biobase:::isValidVersion(object, "CopyNumberSet")) |
|
| 3 |
+## msg <- validMsg(NULL, assayDataValidMembers(assayData(object), c("A", "B")))
|
|
| 4 |
+## if (is.null(msg)) TRUE else msg |
|
| 5 |
+##}) |
|
| 6 |
+ |
|
| 7 |
+setMethod("initialize", "SnpQSet",
|
|
| 8 |
+ function(.Object, |
|
| 9 |
+ assayData = assayDataNew(senseThetaA=senseThetaA, senseThetaB=senseThetaB), |
|
| 10 |
+ senseThetaA=new("matrix"),
|
|
| 11 |
+ senseThetaB=new("matrix"),
|
|
| 12 |
+ phenoData=annotatedDataFrameFrom(assayData, byrow=FALSE), |
|
| 13 |
+ featureData = annotatedDataFrameFrom(assayData, byrow=TRUE), |
|
| 14 |
+ experimentData=new("MIAME"),
|
|
| 15 |
+ annotation=new("character")){
|
|
| 16 |
+ .Object <- callNextMethod(.Object, |
|
| 17 |
+ assayData = assayDataNew(senseThetaA=senseThetaA, senseThetaB=senseThetaB), |
|
| 18 |
+ phenoData=phenoData, |
|
| 19 |
+ experimentData=experimentData, |
|
| 20 |
+ annotation=annotation) |
|
| 21 |
+ .Object |
|
| 22 |
+ }) |
|
| 23 |
+ |
|
| 24 |
+setValidity("SnpQSet",
|
|
| 25 |
+ function(object) |
|
| 26 |
+ assayDataValidMembers(assayData(object), |
|
| 27 |
+ c("senseThetaA",
|
|
| 28 |
+ "senseThetaB") |
|
| 29 |
+ )) |
|
| 30 |
+setMethod("A", "SnpQSet", function(object) senseThetaA(object))
|
|
| 31 |
+setMethod("B", "SnpQSet", function(object) senseThetaB(object))
|
|
| 32 |
+setReplaceMethod("A", signature(object="SnpQSet", value="matrix"),
|
|
| 33 |
+ function(object, value){
|
|
| 34 |
+ assayDataElementReplace(object, "senseThetaA", value) |
|
| 35 |
+ }) |
|
| 36 |
+setReplaceMethod("B", signature(object="SnpQSet", value="matrix"),
|
|
| 37 |
+ function(object, value){
|
|
| 38 |
+ assayDataElementReplace(object, "senseThetaB", value) |
|
| 39 |
+ }) |
|
| 40 |
+setMethod("isSnp", "SnpQSet", function(object) {
|
|
| 41 |
+ isSnp.SnpQSet(object) |
|
| 42 |
+}) |
|
| 43 |
+ |
|
| 44 |
+ |
| ... | ... |
@@ -1,5 +1,5 @@ |
| 1 | 1 |
## Methods for crlmm |
| 2 |
- |
|
| 3 |
- |
|
| 4 | 2 |
setMethod("calls", "SnpSet", function(object) assayData(object)$call)
|
| 3 |
+setReplaceMethod("calls", signature(object="SnpSet", value="matrix"), function(object, value) assayDataElementReplace(object, "call", value))
|
|
| 5 | 4 |
setMethod("confs", "SnpSet", function(object) assayData(object)$callProbability)
|
| 5 |
+setReplaceMethod("confs", signature(object="SnpSet", value="matrix"), function(object, value) assayDataElementReplace(object, "callProbability", value))
|
| ... | ... |
@@ -1,22 +1,39 @@ |
| 1 |
-setMethod("cnIndex", "eSet", function(object, cdfName) match(cnNames(object, cdfName), featureNames(object)))
|
|
| 2 |
-setMethod("cnNames", "eSet", function(object, cdfName) {
|
|
| 3 |
- path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))
|
|
| 1 |
+setMethod("cnIndex", "eSet", function(object){
|
|
| 2 |
+ index <- match(cnNames(object), featureNames(object), nomatch=0) |
|
| 3 |
+ index[index != 0] |
|
| 4 |
+ }) |
|
| 5 |
+ |
|
| 6 |
+setMethod("cnNames", "eSet", function(object) {
|
|
| 7 |
+ path <- system.file("extdata", package=paste(annotation(object), "Crlmm", sep=""))
|
|
| 4 | 8 |
load(file.path(path, "cnProbes.rda")) |
| 9 |
+ cnProbes <- get("cnProbes")
|
|
| 5 | 10 |
cnps <- rownames(cnProbes) |
| 6 | 11 |
cnps <- cnps[cnps %in% featureNames(object)] |
| 7 |
- featureNames(object)[match(cnps, featureNames(object))] |
|
| 8 |
- ##featureNames(object)[grep("CN_", featureNames(object))])
|
|
| 12 |
+ index <- match(cnps, featureNames(object), nomatch=0) |
|
| 13 |
+ index <- index[index != 0] |
|
| 14 |
+ featureNames(object)[index] |
|
| 9 | 15 |
}) |
| 10 |
-setMethod("snpIndex", "eSet", function(object, cdfName) match(snpNames(object, cdfName), featureNames(object)))
|
|
| 11 |
-setMethod("snpNames", "eSet", function(object, cdfName){
|
|
| 12 |
- path <- system.file("extdata", package=paste(cdfName, "Crlmm", sep=""))
|
|
| 16 |
+ |
|
| 17 |
+ |
|
| 18 |
+setMethod("snpIndex", "eSet", function(object){
|
|
| 19 |
+ index <- match(snpNames(object), featureNames(object), nomatch=0) |
|
| 20 |
+ index[index != 0] |
|
| 21 |
+}) |
|
| 22 |
+ |
|
| 23 |
+setMethod("snpNames", "eSet", function(object){
|
|
| 24 |
+ path <- system.file("extdata", package=paste(annotation(object), "Crlmm", sep=""))
|
|
| 13 | 25 |
load(file.path(path, "snpProbes.rda")) |
| 26 |
+ snpProbes <- get("snpProbes")
|
|
| 14 | 27 |
snps <- rownames(snpProbes) |
| 15 | 28 |
snps <- snps[snps %in% featureNames(object)] |
| 16 |
- featureNames(object)[match(snps, featureNames(object))] |
|
| 17 |
- }) |
|
| 29 |
+ index <- match(snps, featureNames(object), nomatch=0) |
|
| 30 |
+ index <- index[index != 0] |
|
| 31 |
+ featureNames(object)[index] |
|
| 32 |
+}) |
|
| 33 |
+ |
|
| 18 | 34 |
setMethod("chromosome", "eSet", function(object) fData(object)$chromosome)
|
| 19 | 35 |
setMethod("position", "eSet", function(object) fData(object)$position)
|
| 36 |
+ |
|
| 20 | 37 |
##setMethod("combine", signature=signature(x="eSet", y="eSet"),
|
| 21 | 38 |
## function(x, y, ...){
|
| 22 | 39 |
## ##Check that both x and y are valid objects |
| ... | ... |
@@ -142,4 +142,19 @@ loader <- function(theFile, envir, pkgname){
|
| 142 | 142 |
load(theFile, envir=envir) |
| 143 | 143 |
} |
| 144 | 144 |
|
| 145 |
+celfileDate <- function(filename) {
|
|
| 146 |
+ h <- affyio::read.celfile.header(filename, info="full") |
|
| 147 |
+ date <- grep("/", strsplit(h$DatHeader, " ")[[1]], value=TRUE)
|
|
| 148 |
+ if(length(date) < 1){
|
|
| 149 |
+ ##try something else |
|
| 150 |
+ results <- h$ScanDate |
|
| 151 |
+ } else{
|
|
| 152 |
+ date <- strsplit(date, split="/")[[1]] |
|
| 153 |
+ CC <- ifelse(substr(date[3],1,1)=="9", "19", "20") |
|
| 154 |
+ results <- as.character(as.Date(paste(paste(CC, date[3], sep=""), date[1], |
|
| 155 |
+ date[2], sep="-"))) |
|
| 156 |
+ } |
|
| 157 |
+ results |
|
| 158 |
+} |
|
| 159 |
+ |
|
| 145 | 160 |
|
| 146 | 161 |
deleted file mode 100644 |
| ... | ... |
@@ -1,113 +0,0 @@ |
| 1 |
-\name{.computeCopynumber}
|
|
| 2 |
-\alias{.computeCopynumber}
|
|
| 3 |
-\title{Internal function for computing copy number}
|
|
| 4 |
-\description{
|
|
| 5 |
- |
|
| 6 |
- This function is not meant to be called directly by the user and is |
|
| 7 |
- not exported in the package namespace. Arguments to this function can |
|
| 8 |
- be specified in the 'computeCopynumber' (no '.' preceding the name) |
|
| 9 |
- function. |
|
| 10 |
- |
|
| 11 |
-} |
|
| 12 |
-\usage{
|
|
| 13 |
-.computeCopynumber(chrom, A, B, calls, conf, NP, plate, MIN.OBS=1, |
|
| 14 |
-envir, P, DF.PRIOR = 50, CONF.THR = 0.99, bias.adj=FALSE, |
|
| 15 |
-priorProb, gender=NULL, SNR, SNRmin, seed=123, cdfName="genomewidesnp6", |
|
| 16 |
-verbose=TRUE, ...) |
|
| 17 |
-} |
|
| 18 |
-\arguments{
|
|
| 19 |
- \item{chrom}{Chromosome (an integer). Use 23 for X and 24 for Y.}
|
|
| 20 |
- \item{A}{The A allele intensities from \code{snprma}}
|
|
| 21 |
- \item{B}{The B allele intensities from \code{snprma}}
|
|
| 22 |
- \item{calls}{The genotype calls from \code{crlmm}}
|
|
| 23 |
- \item{conf}{The genotype confidence scores from \code{crlmm}}
|
|
| 24 |
- \item{NP}{The quantile normalized intensities of the nonpolymorphic probes}
|
|
| 25 |
- \item{plate}{The bach variable. Should be the same length as the
|
|
| 26 |
- number of columns in A} |
|
| 27 |
- \item{MIN.OBS}{Integer: The minimum number of observations in a genotype
|
|
| 28 |
- cluster for which a SNP is deemed complete.} |
|
| 29 |
- \item{envir}{An environment to save intermediate objects}
|
|
| 30 |
- \item{P}{Mainly for debugging a particular plate/batch.}
|
|
| 31 |
- \item{DF.PRIOR}{The degrees of freedom for the prior. Higher numbers
|
|
| 32 |
- will shrink the variance and correlation more.} |
|
| 33 |
- \item{CONF.THR}{A threshold for the genotype confidence scores.
|
|
| 34 |
- Genotypes with scores below the threshold are ignored when computing |
|
| 35 |
- SNP-specific within-genotype estimates of location and scale.} |
|
| 36 |
- \item{bias.adj}{Logical: whether to adjust the location and scale
|
|
| 37 |
- parameters to account for biases due to common copy number |
|
| 38 |
- variants. This is a SNP-specific adjustment. Parameters for |
|
| 39 |
- background and slope must have already been estimated and available |
|
| 40 |
- from the environment variable.} |
|
| 41 |
- \item{priorProb}{Numerical vector of length 4. The prior probability
|
|
| 42 |
- of each copy number state (0, 1, 2, 3, and 4). The default is a |
|
| 43 |
- uniform prior. Ignored if bias.adj=FALSE} |
|
| 44 |
- \item{gender}{Gender of subjects. If not specified, we predict the
|
|
| 45 |
- gender from the X chromosome.} |
|
| 46 |
- \item{SNR}{Signal to noise ratio from crlmm.}
|
|
| 47 |
- \item{SNRmin}{The minimum value for the SNR -- we suggest 5. Samples
|
|
| 48 |
- with SNR below SNRmin are excluded.} |
|
| 49 |
- \item{seed}{Seed used for random samples}
|
|
| 50 |
- \item{cdfName}{Annotation package }
|
|
| 51 |
- \item{verbose}{Logical: verbose output}
|
|
| 52 |
- \item{\dots}{Currently ignored}
|
|
| 53 |
-} |
|
| 54 |
- |
|
| 55 |
-\details{
|
|
| 56 |
- |
|
| 57 |
- This function transforms the intensities to a copy number scale. We |
|
| 58 |
- assume that the median within-SNP intensity across samples is 2. We |
|
| 59 |
- make no assumption about the chromosomal copy number. This function is |
|
| 60 |
- useful for detecting rare variants (e.g., variants that would not |
|
| 61 |
- affect the SNP-specific quantile-based estimators of location and |
|
| 62 |
- scale). A correction for more common variants is coming, as well as |
|
| 63 |
- improved estimates of the uncertainty. |
|
| 64 |
- |
|
| 65 |
-} |
|
| 66 |
- |
|
| 67 |
-\value{
|
|
| 68 |
- |
|
| 69 |
-All objects created by this function are stored in the environment |
|
| 70 |
- passed to this function. In addition, each of the elements are |
|
| 71 |
- specific to the chromosome(s) specified by the argument \code{chrom}.
|
|
| 72 |
- For instance the element \code{A} is the matrix of quantile-normalized
|
|
| 73 |
- intensities for the A-allele on chromosome(s) \code{chrom}. The
|
|
| 74 |
- element of this environment are as follows |
|
| 75 |
- |
|
| 76 |
- \item{A}{Matrix of quantile-normalized intensities for the A-allele}
|
|
| 77 |
- \item{B}{Matrix of quantile-normalized intensities for the A-allele}
|
|
| 78 |
- \item{CA}{Copy number estimate for the A-allele (x 100)}
|
|
| 79 |
- \item{CB}{Copy number estimate for the B-allele (x 100)}
|
|
| 80 |
- \item{calls}{CRLMM genotype calls (AA=1, AB=2, BB=3)}
|
|
| 81 |
- \item{chrom}{Integer(s) indicating the chromosome(s)}
|
|
| 82 |
- \item{cnvs}{Names of the nonpolymorphic probes. These are the
|
|
| 83 |
- rownames of \code{NP} and \code{CT}.}
|
|
| 84 |
- \item{conf}{CRLMM confidence scores for the genotypes: 'round(-1000*log2(1-p))'}
|
|
| 85 |
- \item{corr}{Correlation of the A and B alleles for genotypes AB}
|
|
| 86 |
- \item{corrA.BB}{Correlation of A and B alleles for genotypes BB}
|
|
| 87 |
- \item{corrB.AA}{Correlation of A and B alleles for genotypes AA}
|
|
| 88 |
- \item{CT}{Copy number estimates for nonpolymorphic probe. See
|
|
| 89 |
- \code{cnvs} for the rownames.}
|
|
| 90 |
- \item{CT.sds}{Standard deviation estimates for \code{CT}}
|
|
| 91 |
- \item{npflags}{Flags for the nonpolymorphic probes.}
|
|
| 92 |
- \item{Ns}{The number of observations for each genotype/plate}
|
|
| 93 |
- \item{nuA}{Background/cross-hyb for the A allele (plate- and locus-specific)}
|
|
| 94 |
- \item{nuB}{Background/cross-hyb for the B allele (plate- and locus-specific)}
|
|
| 95 |
- \item{nuT}{Background for the nonpolymorphic probes (plate- and locus-specific)}
|
|
| 96 |
- \item{phiA}{Slope for the A allele (plate- and locus-specific)}
|
|
| 97 |
- \item{phiB}{Slope for the B allele (plate- and locus-specific)}
|
|
| 98 |
- \item{phiT}{Slope for the nonpolymorphic probes (plate- and locus-specific)}
|
|
| 99 |
- \item{plate}{Factor indicating batch (same length as number of cel files)}
|
|
| 100 |
- \item{sig2A}{Variance estimate for the A-allele signal (plate- and locus-specific)}
|
|
| 101 |
- \item{sig2B}{Variance estimate for the B-allele signal (plate- and locus-specific)}
|
|
| 102 |
- \item{sig2T}{Variance estimate for the nonpolymorphic signal (plate- and locus-specific)}
|
|
| 103 |
- \item{snpflags}{Flags for polymorphic probes}
|
|
| 104 |
- \item{snps}{Rownames for \code{A}, \code{B}, \code{CA}, \code{CB}, ...}
|
|
| 105 |
- \item{sns}{ Sample names -- the column names for \code{A}, \code{B}, ...}
|
|
| 106 |
- \item{steps}{Steps completed. For internal use.}
|
|
| 107 |
- \item{tau2A}{Variance estimate for the B-allele background/cross-hyb (plate- and locus-specific)}
|
|
| 108 |
- \item{tau2B}{Variance estimate for the B-allele background/cross-hyb (plate- and locus-specific)}
|
|
| 109 |
-} |
|
| 110 |
-\references{Nothing yet.}
|
|
| 111 |
-\author{Rob Scharpf}
|
|
| 112 |
-\keyword{manip}
|
|
| 113 |
- |
| 114 | 0 |
deleted file mode 100644 |
| ... | ... |
@@ -1,55 +0,0 @@ |
| 1 |
-\name{ABset-class}
|
|
| 2 |
-\Rdversion{1.1}
|
|
| 3 |
-\docType{class}
|
|
| 4 |
-\alias{ABset-class}
|
|
| 5 |
-\alias{A}
|
|
| 6 |
-\alias{B}
|
|
| 7 |
-\alias{A,ABset-method}
|
|
| 8 |
-\alias{B,ABset-method}
|
|
| 9 |
-\alias{A<-,ABset,matrix-method}
|
|
| 10 |
-\alias{B<-,ABset,matrix-method}
|
|
| 11 |
-\title{Class "ABset"}
|
|
| 12 |
-\description{Containter for quantile-normalized intensities}
|
|
| 13 |
-\section{Objects from the Class}{
|
|
| 14 |
-Objects can be created by calls of the form \code{new("ABset", assayData, phenoData, featureData, experimentData, annotation, protocolData, ...)}.
|
|
| 15 |
-} |
|
| 16 |
-\section{Slots}{
|
|
| 17 |
- \describe{
|
|
| 18 |
- \item{\code{assayData}:}{Object of class \code{"AssayData"} }
|
|
| 19 |
- \item{\code{phenoData}:}{Object of class \code{"AnnotatedDataFrame"} }
|
|
| 20 |
- \item{\code{featureData}:}{Object of class \code{"AnnotatedDataFrame"} }
|
|
| 21 |
- \item{\code{experimentData}:}{Object of class \code{"MIAME"} }
|
|
| 22 |
- \item{\code{annotation}:}{Object of class \code{"character"} }
|
|
| 23 |
- \item{\code{protocolData}:}{Object of class \code{"AnnotatedDataFrame"} }
|
|
| 24 |
- \item{\code{.__classVersion__}:}{Object of class \code{"Versions"} }
|
|
| 25 |
- } |
|
| 26 |
-} |
|
| 27 |
-\section{Extends}{
|
|
| 28 |
-Class \code{\link[Biobase:class.eSet]{eSet}}, directly.
|
|
| 29 |
-Class \code{\link[Biobase:class.VersionedBiobase]{VersionedBiobase}}, by class "eSet", distance 2.
|
|
| 30 |
-Class \code{\link[Biobase:class.Versioned]{Versioned}}, by class "eSet", distance 3.
|
|
| 31 |
-} |
|
| 32 |
-\section{Methods}{
|
|
| 33 |
- \describe{
|
|
| 34 |
- \item{A}{\code{signature(object="ABset")}: accessor for the
|
|
| 35 |
- quantile-normalized intensities of allele A for polymorphic probes |
|
| 36 |
- and the quantile normalized intensities for the copy number |
|
| 37 |
- probes.} |
|
| 38 |
- |
|
| 39 |
- \item{"A<-"}{\code{signature(object="ABset", value="matrix")}:
|
|
| 40 |
- replacement method for the A allele intensities.} |
|
| 41 |
- |
|
| 42 |
- \item{"B<-"}{\code{signature(object="ABset", value="matrix")}:
|
|
| 43 |
- replacement method for the B allele intensities.} |
|
| 44 |
- |
|
| 45 |
- \item{B}{\code{signature(object="ABset")}: accessor for the
|
|
| 46 |
- quantile-normalized intensities of allele B for polymorphic probes.} |
|
| 47 |
- \item{plot}{\code{signature(x = "ABset", y = "CopyNumberSet")}:
|
|
| 48 |
-... }} |
|
| 49 |
- |
|
| 50 |
-} |
|
| 51 |
-\author{R. Scharpf}
|
|
| 52 |
-\examples{
|
|
| 53 |
-showClass("ABset")
|
|
| 54 |
-} |
|
| 55 |
-\keyword{classes}
|
| 56 | 0 |
deleted file mode 100644 |
| ... | ... |
@@ -1,194 +0,0 @@ |
| 1 |
-\name{CrlmmSetList-class}
|
|
| 2 |
-\Rdversion{1.1}
|
|
| 3 |
-\docType{class}
|
|
| 4 |
-\alias{CrlmmSetList-class}
|
|
| 5 |
-\alias{[,CrlmmSetList-method}
|
|
| 6 |
-\alias{$,CrlmmSetList-method}
|
|
| 7 |
-\alias{addFeatureAnnotation}
|
|
| 8 |
-\alias{addFeatureAnnotation,CrlmmSetList-method}
|
|
| 9 |
-\alias{A<-}
|
|
| 10 |
-\alias{B<-}
|
|
| 11 |
-\alias{A,CrlmmSetList-method}
|
|
| 12 |
-\alias{A<-,CrlmmSetList,matrix-method}
|
|
| 13 |
-\alias{B,CrlmmSetList-method}
|
|
| 14 |
-\alias{B<-,CrlmmSetList,matrix-method}
|
|
| 15 |
-\alias{batch,CrlmmSetList-method}
|
|
| 16 |
-\alias{CA,CrlmmSetList-method}
|
|
| 17 |
-\alias{CB,CrlmmSetList-method}
|
|
| 18 |
-\alias{CA<-,CrlmmSetList,matrix-method}
|
|
| 19 |
-\alias{CB<-,CrlmmSetList,matrix-method}
|
|
| 20 |
-\alias{calls,CrlmmSetList-method}
|
|
| 21 |
-\alias{chromosome,CrlmmSetList-method}
|
|
| 22 |
-\alias{cnIndex,CrlmmSetList-method}
|
|
| 23 |
-\alias{combine,CrlmmSetList,CrlmmSetList-method}
|
|
| 24 |
-\alias{confs,CrlmmSetList-method}
|
|
| 25 |
-\alias{copyNumber,CrlmmSetList-method}
|
|
| 26 |
-\alias{nrow,CrlmmSetList-method}
|
|
| 27 |
-\alias{ncol,CrlmmSetList-method}
|
|
| 28 |
-\alias{plot,CrlmmSetList,ANY-method}
|
|
| 29 |
-\alias{points,CrlmmSetList-method}
|
|
| 30 |
-\alias{position,CrlmmSetList-method}
|
|
| 31 |
-\alias{sampleNames,CrlmmSetList-method}
|
|
| 32 |
-\alias{show,CrlmmSetList-method}
|
|
| 33 |
-\alias{snpIndex,CrlmmSetList-method}
|
|
| 34 |
-\alias{update,CrlmmSetList-method}
|
|
| 35 |
-\alias{update,character-method}
|
|
| 36 |
- |
|
| 37 |
- |
|
| 38 |
-\title{Class "CrlmmSetList"}
|
|
| 39 |
-\description{Container for quantile normalized intensities and genotype calls.}
|
|
| 40 |
-\section{Objects from the Class}{
|
|
| 41 |
- Objects from the class are created by calls to the \code{crlmmWrapper}
|
|
| 42 |
- function. |
|
| 43 |
-} |
|
| 44 |
-\section{Slots}{
|
|
| 45 |
- \describe{
|
|
| 46 |
- \item{\code{.Data}:}{Object of class \code{"list"} ~~ }
|
|
| 47 |
- } |
|
| 48 |
-% \item{\code{ABset}:}{Object of class \code{"ABset"}.
|
|
| 49 |
-% Contains quantile normalized A and B intensities. Quantile |
|
| 50 |
-% normalized intensities for nonpolymorphic probes are stored |
|
| 51 |
-% in the A slot.} |
|
| 52 |
-% \item{\code{crlmmResult}:}{Object of class
|
|
| 53 |
-% \code{"SnpSet"}. Contains genotype calls and the
|
|
| 54 |
-% corresponding confidence estimates.} |
|
| 55 |
-} |
|
| 56 |
- |
|
| 57 |
-\section{Details}{
|
|
| 58 |
- |
|
| 59 |
- Instances of \code{CrlmmSetList} are a list with two elements:
|
|
| 60 |
- |
|
| 61 |
- 1. an object of class \code{ABset}
|
|
| 62 |
- |
|
| 63 |
- 2. an object of class \code{SnpSet}
|
|
| 64 |
- |
|
| 65 |
- The \code{featureNames} and \code{sampleNames} of both objects are
|
|
| 66 |
- required to be identical. |
|
| 67 |
- |
|
| 68 |
- Quantile-normalized intensities for the copy number probes are stored |
|
| 69 |
- in the assay data element A of the ABset object. Corresponding rows |
|
| 70 |
- for the B allele are recorded as NAs. Genotype calls for copy number |
|
| 71 |
- probes in the SnpSet object are recorded as NAs. |
|
| 72 |
- |
|
| 73 |
-} |
|
| 74 |
- |
|
| 75 |
-\section{Extends}{
|
|
| 76 |
- Class \code{"\linkS4class{list}"}, from data part.
|
|
| 77 |
- Class \code{"\linkS4class{vector}"}, by class "list", distance 2.
|
|
| 78 |
- Class \code{\link[Biobase:assayData]{assayData}}, by class
|
|
| 79 |
- "list", distance 2. |
|
| 80 |
-} |
|
| 81 |
-\section{Methods}{
|
|
| 82 |
- \describe{
|
|
| 83 |
- \item{"["}{\code{signature(x = "CrlmmSetList")}: subsets both the
|
|
| 84 |
- features and samples of each element in th elist} |
|
| 85 |
- |
|
| 86 |
- \item{"$"}{\code{signature(x = "CrlmmSetList")}: used to access
|
|
| 87 |
- element from the \code{featureData} of the \code{CopyNumberSet}
|
|
| 88 |
- element. In particular, useful for extract the SNP- and |
|
| 89 |
- batch-specific parameters used to estimate copy number that are |
|
| 90 |
- currently stored in the \code{featureData} of \code{CopyNumberSet}
|
|
| 91 |
- objects. } |
|
| 92 |
- |
|
| 93 |
- \item{"addFeatureAnnotation"}{\code{signature(object =
|
|
| 94 |
- "CrlmmSetList")}: Creates an object of class AnnotatedDataFrame from |
|
| 95 |
- the CrlmmSetList object. The new annotation object contains |
|
| 96 |
- information on chromosome and physical position and has the same |
|
| 97 |
- ordering of rows as the CrlmmSetList object. } |
|
| 98 |
- |
|
| 99 |
- \item{"A"}{\code{signature(object = "CrlmmSetList")}: extracts the
|
|
| 100 |
- quantile normalized intensities for the A allele in the |
|
| 101 |
- \code{ABset} element.}
|
|
| 102 |
- |
|
| 103 |
- \item{"A<-"}{\code{signature(object = "CrlmmSetList", value =
|
|
| 104 |
- "matrix")}: Replacement method for the A intensities.} |
|
| 105 |
- |
|
| 106 |
- \item{"B"}{\code{signature(object = "CrlmmSetList")}: extracts the
|
|
| 107 |
- quantile normalized intensities for the B allele in the |
|
| 108 |
- \code{ABset} element.}
|
|
| 109 |
- |
|
| 110 |
- \item{"B<-"}{\code{signature(object = "CrlmmSetList", value =
|
|
| 111 |
- "matrix")}: Replacement method for the B intensities.} |
|
| 112 |
- |
|
| 113 |
- \item{"batch"}{\code{signature(object = "CrlmmSetList")}: extracts the
|
|
| 114 |
- batch information used to estimate copy number.} |
|
| 115 |
- |
|
| 116 |
- \item{"CA"}{\code{signature(object = "CrlmmSetList")}: extracts the
|
|
| 117 |
- copy number for allele A at polymorphic loci. For nonpolymorphic |
|
| 118 |
- probes, CA returns the total copy number. } |
|
| 119 |
- |
|
| 120 |
- \item{"CB"}{\code{signature(object = "CrlmmSetList")}: extracts the
|
|
| 121 |
- copy number for allele B at polymorphic loci. For nonpolymorphic |
|
| 122 |
- probes, CB returns 'NA'.} |
|
| 123 |
- |
|
| 124 |