@@ -5,6 +5,7 @@

     - added calculateGCContentByInterval

     - added calculatePowerDetectSomatic

     - added callAlterations

+    - added callLOH

     - added getDiploid

     - added getSexFromVcf

@@ -1,48 +1,48 @@

-export(runAbsoluteCN)

-export(filterVcfBasic)

-export(bootstrapResults)

-export(filterVcfMuTect)

-export(setPriorVcf)

-export(getDiploid)

 export(autoCurateResults)

-export(segmentationCBS)

-export(segmentationPSCBS)

-export(poolCoverage)

-export(readCoverageGatk)

-export(correctCoverageBias)

-export(callAlterations)

-export(findFocal)

-export(readCurationFile)

-export(calculatePowerDetectSomatic)

+export(bootstrapResults)

 export(calculateBamCoverageByInterval)

 export(calculateGCContentByInterval)

+export(calculatePowerDetectSomatic)

+export(callAlterations)

+export(callLOH)

+export(correctCoverageBias)

 export(createCurationFile)

+export(createExonWeightFile)

 export(createNormalDatabase)

 export(createSNPBlacklist)

-export(createExonWeightFile)

+export(filterVcfBasic)

+export(filterVcfMuTect)

 export(findBestNormal)

-export(plotBestNormal)

-export(predictSomatic)

-export(plotAbs)

+export(findFocal)

+export(getDiploid)

 export(getSexFromCoverage)

 export(getSexFromVcf)

-import(VariantAnnotation)

+export(plotAbs)

+export(plotBestNormal)

+export(poolCoverage)

+export(predictSomatic)

+export(readCoverageGatk)

+export(readCurationFile)

+export(runAbsoluteCN)

+export(segmentationCBS)

+export(segmentationPSCBS)

+export(setPriorVcf)

 import(DNAcopy)

-importFrom("data.table", "data.table")

-importFrom("RColorBrewer", "brewer.pal")

+import(IRanges)

+import(VariantAnnotation)

+importFrom("Biostrings", "letterFrequency")

 importFrom("GenomeInfoDb", "seqnames")

 importFrom("GenomicRanges", "GRanges")

-importFrom("SummarizedExperiment", "rowRanges")

+importFrom("RColorBrewer", "brewer.pal")

 importFrom("Rsamtools", "ScanBamParam", "scanBamFlag", "scanBam", "scanFa")

-import(IRanges)

 importFrom("S4Vectors", "queryHits", "subjectHits")

+importFrom("SummarizedExperiment", "rowRanges")

+importFrom("data.table", "data.table")

 importFrom("grDevices", "colorRampPalette", "dev.off", "png", "adjustcolor")

-importFrom("graphics", "abline", "axis", "boxplot", "contour", "hist",

-    "image", "legend", "lines", "par", "plot", "text", "mtext", "polygon",

-    "points")

-importFrom("stats", "C", "complete.cases", "dbeta", "dist", "dnorm",

-    "dunif", "loess", "pbeta", "prcomp", "predict", "runif", "t.test",

-    "weighted.mean", "dbinom", "fisher.test", "hclust", "cutree")

-importFrom("utils", "data", "head", "read.csv", "read.delim",

-    "read.table", "write.csv", "write.table", "tail", "packageVersion")

-importFrom("Biostrings", "letterFrequency")

+importFrom("graphics", "abline", "axis", "boxplot", "contour", "hist","image", 

+    "legend", "lines", "par", "plot", "text", "mtext", "polygon","points")

+importFrom("stats", "C", "complete.cases", "dbeta", "dist", "dnorm", "dunif", 

+    "loess", "pbeta", "prcomp", "predict", "runif", "t.test", "weighted.mean", 

+    "dbinom", "fisher.test", "hclust", "cutree")

+importFrom("utils", "data", "head", "read.csv", "read.delim", "read.table", 

+    "write.csv", "write.table", "tail", "packageVersion")

new file mode 100644

@@ -0,0 +1,91 @@

+callLOH <- structure(function(# Get regions of LOH

+### This function provides detailed LOH information by region.

+res, 

+### Return object of the runAbsoluteCN() function.

+id=1, 

+### Candidate solution to extract LOH from. id=1 will  

+### use the maximum likelihood solution.

+arm.cutoff=0.9

+### Min fraction LOH on a chromosome arm to call 

+### whole arm events.

+) {

+    if (is.null(res$input$vcf)) {

+        .stopUserError("runAbsoluteCN was run without a VCF file.")

+    }

+    chr.hash <- res$input$chr.hash

+    armLocations <- .getArmLocations(res, 1)

+    loh <- res$results[[id]]$SNV.posterior$beta.model$loh$output

+    if (!nrow(armLocations)) {

+        warning("Centromere positions not available or matching.")

+        return(loh)

+    }

+    armLocationsGR <- GRanges(seqnames=armLocations$chrom,

+IRanges(start=armLocations$start, end=armLocations$end))

+    seg <- res$results[[id]]$seg

+

+    minorChrNumber <- cbind(

+        res$results[[id]]$SNV.posterior$beta.model$segment.ids,

+        res$results[[id]]$SNV.posterior$beta.model$posteriors$ML.M.Segment

+    )

+    minorChrNumber <- minorChrNumber[!duplicated(minorChrNumber[,1]),]

+    seg$M[minorChrNumber[,1]] <- minorChrNumber[,2]

+    seg <- seg[complete.cases(seg),]

+    seg$chrom <- .add.chr.name(seg$chrom, chr.hash)

+    segGR <- GRanges(seqnames=seg$chrom, IRanges(start=seg$loc.start,

+    end=seg$loc.end))

+

+    ov <- findOverlaps(armLocationsGR, segGR)

+    segLOH <-  cbind(seg[subjectHits(ov),], armLocations[queryHits(ov),])

+    segLOH$loc.start <- apply(segLOH[,c("loc.start", "start")],1,max)

+    segLOH$loc.end <- apply(segLOH[,c("loc.end", "end")],1,min)

+    segLOH$chrom <- paste(segLOH$chrom, segLOH$arm, sep="")

+    segLOH$size <- segLOH$loc.end-segLOH$loc.start+1

+    segLOH$fraction.arm <- round(segLOH$size/armLocations$size[match(segLOH$chrom,

+        paste(armLocations$chrom, armLocations$arm, sep=""))], digits=2)

+    segLOH$type <- ""

+    segLOH$type[segLOH$C %% 2 == 0 & segLOH$C > 1 & 

+        segLOH$M == 0] <- "COPY-NEUTRAL LOH"

+    segLOH$type[segLOH$type=="" & segLOH$M == 0] <- "LOH"

+    idx <- segLOH$fraction.arm > arm.cutoff & segLOH$type != ""

+    segLOH$type[idx] <- paste("WHOLE ARM",

+        segLOH$type)[idx]

+    rownames(segLOH) <- NULL

+    segLOH <- segLOH[, c("chrom", "loc.start", "loc.end", "C", "M", "type")]

+    # standardize colnames

+    colnames(segLOH)[1:3] <- c("chr", "start", "end")

+    segLOH

+### Returns data.frame with LOH regions.    

+}, ex=function() {

+data(purecn.example.output)

+head(callLOH(purecn.example.output))

+})

+

+.getArmLocations <- function(res, id) {

+    chr.hash <- res$input$chr.hash

+    centromeres <- res$input$centromeres

+

+    chromCoords <- t(vapply(split(

+        res$results[[id]]$SNV.posterior$beta.model$loh$data$maploc,

+        res$results[[id]]$SNV.posterior$beta.model$loh$data$chrom), function(x)

+        c(min(x), max(x)), c(min=double(1), max=double(1))))

+

+    centromeres <- cbind(centromeres, chromCoords[match(centromeres$chrom,

+        rownames(chromCoords)),])

+    centromeres <- centromeres[complete.cases(centromeres),]

+

+    pArms <- centromeres[centromeres$min < centromeres$chromStart,

+        c("chrom", "min", "chromStart")]

+    qArms <- centromeres[centromeres$max > centromeres$chromEnd,

+        c("chrom", "chromEnd", "max")]

+    

+    colnames(pArms) <- c("chrom", "start", "end")

+    colnames(qArms) <- colnames(pArms)

+    pArms$arm <- "p"

+    qArms$arm <- "q"

+    armLocations <- rbind(pArms, qArms)

+    armLocations <- armLocations[order(match(armLocations$chrom,

+        chr.hash[,1])),]

+    armLocations$size <- armLocations$end-armLocations$start+1

+    armLocations

+}

+

new file mode 100755

@@ -0,0 +1,6 @@

+test_callLOH <- function() {

+    data(purecn.example.output)

+    ret <- callLOH(purecn.example.output)

+    checkEquals("data.frame", class(ret))

+    checkEqualsNumeric(6, ncol(ret))

+}    

new file mode 100644

@@ -0,0 +1,24 @@

+\name{callLOH}

+\alias{callLOH}

+\title{Get regions of LOH}

+\description{This function provides detailed LOH information by region.}

+\usage{callLOH(res, id = 1, arm.cutoff = 0.9)}

+\arguments{

+  \item{res}{Return object of the runAbsoluteCN() function.}

+  \item{id}{Candidate solution to extract LOH from. id=1 will  

+use the maximum likelihood solution.}

+  \item{arm.cutoff}{Min fraction LOH on a chromosome arm to call 

+whole arm events.}

+}

+

+\value{Returns data.frame with LOH regions.    }

+

+\author{Markus Riester}

+

+

+

+

+\examples{

+data(purecn.example.output)

+head(callLOH(purecn.example.output))

+}

@@ -615,6 +615,17 @@ gene.calls <- callAlterations(ret)

 head(gene.calls)

 @

+\section{Find genomic regions in LOH}

+

+The \Rcode{gene.calls} data.frame described above provides gene-level LOH

+information. To find the corresponding regions in LOH, the \Rfunction{getLOH}

+provides a convenient way:

+

+<<loh>>=

+loh <- callLOH(ret)

+head(loh)

+@

+  

 \section{Power to detect somatic mutations}

 As final quality control step, we can test if coverage and tumor purity are