@@ -2,16 +2,15 @@ Package: PureCN

 Type: Package

 Title: Estimating tumor purity, ploidy, LOH, and SNV status using

         hybrid capture NGS data

-Version: 0.99.6

+Version: 0.99.7

 Date: 2016-03-30

 Author: Markus Riester

 Maintainer: Markus Riester <[email protected]>

 Description: This package estimates tumor purity, copy number, loss of

-        heterozygosity (LOH), and status of short nucleotide variants

-        (SNVs). PureCN is designed for hybrid capture next generation

-        sequencing (NGS) data, integrates well with standard somatic

-        variant detection pipelines, and has support for tumor samples

-        without matching normal samples.

+    heterozygosity (LOH), and status of short nucleotide variants (SNVs). 

+    PureCN is designed for hybrid capture next generation sequencing (NGS) 

+    data, integrates well with standard somatic variant detection pipelines, 

+    and has support for tumor samples without matching normal samples.

 Depends: R (>= 3.3), DNAcopy, VariantAnnotation (>= 1.14.1)

 Imports: GenomicRanges (>= 1.20.3), IRanges (>= 2.2.1), RColorBrewer,

         S4Vectors, data.table, grDevices, graphics, stats, utils,

@@ -22,3 +21,4 @@ License: Artistic-2.0

 biocViews: CopyNumberVariation, Software, Sequencing,

         VariantAnnotation, VariantDetection

 NeedsCompilation: no

+Packaged: 2016-04-14 21:52:04 UTC; riestma1

@@ -17,6 +17,7 @@ export(findBestNormal)

 export(plotBestNormal)

 export(predictSomatic)

 export(plotAbs)

+export(getSexFromCoverage)

 import(VariantAnnotation)

 import(DNAcopy)

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

@@ -86,11 +86,6 @@ max.exon.ratio) {

                     Mi + g * (1 - p))/(p * lr.C[j] + 2 * (1 - p)), sd = sd.ar, log = TRUE) - 

                     log(lr.C[j] + 1), -Inf))))

             } else {

-                # p.ar <- lapply(c(0,1), function(g) lapply(1:length(ar_all), function(j)

-                # sapply(test.num.copy, function(Mi) ifelse(Mi <= lr.C[j], dbeta(x = (p * Mi + g

-                # * (1 - p))/(p * lr.C[j] + 2 * (1 - p)), shape1=ar_all[j]*dp_all[j]+1,

-                # shape2=(1-ar_all[j])*dp_all[j]+1, log = TRUE)-log(lr.C[j]+1+haploid.penalty) ,

-                # -Inf))))

                 p.ar <- lapply(c(0, 1), function(g) lapply(1:length(ar_all), function(j) sapply(test.num.copy, 

                   function(Mi) ifelse(Mi <= lr.C[j], dbeta(x = (p * Mi + g * (1 - 

                     p))/(p * lr.C[j] + 2 * (1 - p)), shape1 = ar_all[j] * dp_all[j] + 

@@ -265,7 +260,7 @@ max.exon.ratio) {

     return(result)

 }

-.flagResults <- function(results, max.non.clonal = 0.2, max.logr.sdev, logr.sdev, 

+.flagResults <- function(results, max.non.clonal = 0.2, max.logr.sdev, logr.sdev, max.segments,

     flag = NA, flag_comment = NA) {

     results <- lapply(results, .flagResult, max.non.clonal = max.non.clonal)

@@ -273,6 +268,7 @@ max.exon.ratio) {

     # results[[2]]$total.log.likelihood ) / abs( results[[1]]$total.log.likelihood )

     # if (ldiff < 0.1) { results[[1]]$flag <- TRUE results[[1]]$flag_comment <-

     # 'SIMILAR TO SECOND MOST LIKELY OPTIMUM' return(results) }

+    number.segments <- nrow(results[[1]]$seg)

     if (logr.sdev > max.logr.sdev) {

         for (i in 1:length(results)) {

@@ -281,6 +277,14 @@ max.exon.ratio) {

                 "NOISY LOG-RATIO")

         }

     }

+

+    if (number.segments > max.segments) {

+        for (i in 1:length(results)) {

+            results[[i]]$flag <- TRUE

+            results[[i]]$flag_comment <- .appendComment(results[[i]]$flag_comment, 

+                "NOISY SEGMENTATION")

+        }

+    }

     # some global flags

     if (!is.na(flag) && flag) {

@@ -519,4 +523,10 @@ max.exon.ratio) {

     data(chr.hash, envir = environment())

     chr.hash[as.character(ls), 2]

 }

+.add.chr.name <- function(ls) {

+    chr.hash <- NULL

+    data(chr.hash, envir = environment())

+    as.character(chr.hash$chr[match(ls, chr.hash$number)])

+}

+    

@@ -1,10 +1,13 @@

 createNormalDatabase <- structure(function(

-### Function to create a database of normal samples, used to find a good match for tumor copy number normalization.

+### Function to create a database of normal samples, used to find 

+### a good match for tumor copy number normalization.

 gatk.normal.files,

-### Vector with file names pointing to GATK coverage files of normal samples. 

+### Vector with file names pointing to GATK coverage files 

+### of normal samples. 

 ...

 ### Arguments passed to the prcomp function.

 ) {

+    gatk.normal.files <- normalizePath(gatk.normal.files)

     normals <- lapply(gatk.normal.files, readCoverageGatk)

     normals.m <- do.call(cbind, lapply(normals, function(x) x$average.coverage))

     idx <- complete.cases(normals.m)

@@ -15,7 +18,8 @@ gatk.normal.files,

         exons.used=idx, 

         coverage=apply(normals.m, 2, mean, na.rm=TRUE), 

         exon.median.coverage=apply(normals.m,1,median, na.rm=TRUE),

-        exon.log2.sd.coverage=apply(log2(normals.m+1),1,sd, na.rm=TRUE)

+        exon.log2.sd.coverage=apply(log2(normals.m+1),1,sd, na.rm=TRUE),

+        sex=sapply(normals,getSexFromCoverage)

     )

 ### A normal database that can be used in the findBestNormal function to 

 ### retrieve good matching normal samples for a given tumor sample.

@@ -27,11 +31,13 @@ normalDB <- createNormalDatabase(gatk.normal.files)

 })    

 createExonWeightFile <- structure(function(# Calculate exon weights

-### Creates an exon weight file useful for segmentation, by down-weighting unreliable exons.

+### Creates an exon weight file useful for segmentation, by 

+### down-weighting unreliable exons.

 gatk.tumor.files, 

 ### A small number (1-3) of GATK tumor coverage samples.

 gatk.normal.files,

-### A large number of GATK normal coverage samples (>20) to estimate exon log-ratio standard deviations.

+### A large number of GATK normal coverage samples (>20) 

+### to estimate exon log-ratio standard deviations.

 exon.weight.file

 ### Output filename.

 ) {

@@ -51,7 +57,7 @@ exon.weight.file

     ret <- data.frame(Target=tumor.coverage[[1]][,1], Weights=zz)

     write.table(ret, file=exon.weight.file,row.names=FALSE, quote=FALSE, sep="\t")

-    ret

+    invisible(ret)

 ###A data.frame with exon weights.

 }, ex=function() {

 exon.weight.file <- "exon_weights.txt"

@@ -7,13 +7,16 @@ normalDB,

 pcs=1:3,

 ### Principal components to use for distance calculation.

 num.normals=1,

-## Return the num.normals best normals.

+### Return the num.normals best normals.

+ignore.sex=FALSE,

+### If FALSE, detects sex of sample returns a best normal 

+### with matching sex.

 verbose=TRUE

+### Verbose output.

 ) {

     if (is.character(gatk.tumor.file)) {

         tumor  <- readCoverageGatk(gatk.tumor.file)

     } else {

-        if (verbose) message("gatk.tumor.file does not appear to be a filename, assuming it is valid GATK coverage data.")

         tumor <- gatk.tumor.file

     }    

     x <- t(tumor[normalDB$exons.used,"average.coverage", drop=FALSE])

@@ -22,8 +25,22 @@ verbose=TRUE

     idx.pcs <- pcs

     idx.pcs <- idx.pcs[idx.pcs %in% 1:ncol(normalDB$pca$x)]

-    best.match <- order(sapply(1:length(normalDB$gatk.normal.files), function(i) dist(  rbind(predict(normalDB$pca, x)[1,idx.pcs], predict(normalDB$pca)[i,idx.pcs]))[1]))

-    normalDB$gatk.normal.files[head(best.match, num.normals)]

+    idx.normals <- seq_along(normalDB$gatk.normal.files)

+

+    if (!ignore.sex && !is.null(normalDB$sex) && sum(!is.na(normalDB$sex))>0) {

+        sex <- getSexFromCoverage(tumor, verbose=FALSE)

+        if (verbose) message(paste("Sex of sample:", sex))

+        if (!is.na(sex)) {

+            idx.normals <- which(normalDB$sex == sex)

+        }

+        if (length(idx.normals) < 2) { 

+            warning(paste("Not enough samples of sex", sex, "in database. Ignoring sex."))

+            idx.normals <- seq_along(normalDB$gatk.normal.files)

+        }    

+    }    

+

+    best.match <- order(sapply(idx.normals, function(i) dist(  rbind(predict(normalDB$pca, x)[1,idx.pcs], predict(normalDB$pca)[i,idx.pcs]))[1]))

+    normalDB$gatk.normal.files[idx.normals][head(best.match, num.normals)]

 ### Return the filename of the best matching normal    

 },ex=function() {

 gatk.normal.file <- system.file("extdata", "example_normal.txt", package="PureCN")

new file mode 100755

@@ -0,0 +1,42 @@

+getSexFromCoverage <- function(# Get sample sex from coverage

+### This function determines the sex of a sample by the coverage 

+### ratio of chrX and chrY.

+gatk.coverage, 

+### GATK coverage file or data read with readCoverageGatk.

+min.ratio=20,

+### Min chrX/chrY coverage ratio to call sample as female.

+verbose=TRUE

+### Verbose output.

+) {

+    if (is.character(gatk.coverage)) {

+        x <- readCoverageGatk(gatk.coverage)

+    } else {

+        x <- gatk.coverage

+    }

+

+    sex.chr <- .getSexChr(x)

+    xx <- split(x$average.coverage, x$chr)

+    avg.coverage <- sapply(xx, mean, na.rm=TRUE)

+    if (is.na(avg.coverage[sex.chr[1]]) || is.na(avg.coverage[sex.chr[2]]) ) {

+        if (verbose) message("Allosome coverage appears to be missing, cannot determine sex.")

+        return(NA)

+    }    

+

+    autosome.ratio <- mean(avg.coverage[-match(sex.chr, names(avg.coverage))], na.rm=TRUE)/(avg.coverage[sex.chr[1]]+0.0001)

+    if (autosome.ratio > 5) { 

+        if (verbose) message("Allosome coverage very low, cannot determine sex.")

+        return(NA)

+    }

+    XY.ratio <- avg.coverage[sex.chr[1]]/ (avg.coverage[sex.chr[2]]+ 0.0001)

+    if (XY.ratio > min.ratio) return("F")

+    return("M")    

+### Returns "M" for male, "F" for female, or NULL if unknown.    

+}

+

+.getSexChr <- function(gatk.coverage) {

+    if ("chrX" %in% gatk.coverage$chr) {

+        return(c("chrX", "chrY"))

+    }

+    return(as.character(23:24))    

+}

+       

@@ -1,35 +1,67 @@

 runAbsoluteCN <-

-structure(function(# Run our implementation of ABSOLUTE

-### This function takes as input tumor and normal control coverage and allelic fractions of germline variants and somatic mutations.

-### Coverage data is provied in GATK DepthOfCoverage format, allelic fraction in VCF format (e.g. obtained by MuTect). 

-### Normal control does not need to be matched (from the same patient). In case VCF does not contain somatic status, it should contain

-### dbSNP and optionally COSMIC annotation.

+structure(function(# Run PureCN implementation of ABSOLUTE

+### This function takes as input tumor and normal control coverage and 

+### allelic fractions of germline variants and somatic mutations.

+### Coverage data is provied in GATK DepthOfCoverage format, allelic fraction 

+### in VCF format (e.g. obtained by MuTect). Normal control does not need to 

+### be matched (from the same patient). In case VCF does not contain somatic 

+### status, it should contain dbSNP and optionally COSMIC annotation.

 gatk.normal.file=NULL, 

-### GATK coverage file of normal control (optional if log.ratio is provided - then it will be only used to filter low coverage exons). Should be already GC-normalized.

+### GATK coverage file of normal control (optional if 

+### log.ratio is provided - then it will be only used to filter low coverage 

+### exons). Should be already GC-normalized. Needs to be either a file name 

+### or data read with the readCoverageGatk function.

 gatk.tumor.file, 

-### GATK coverage file of tumor. Should be already GC-normalized.

+### GATK coverage file of tumor. Should be already 

+### GC-normalized. Needs to be either a file name or data read with the 

+### readCoverageGatk function.

 log.ratio=NULL, 

-### Copy number log-ratios for all exons in the coverage files. If NULL, calculated based on coverage files.

+### Copy number log-ratios for all exons in the coverage files. 

+### If NULL, calculated based on coverage files.

 seg.file=NULL,

-### Segmented data. Optional, to support matched SNP6 data. If null, use coverage files or log.ratio to segment the data.  

+### Segmented data. Optional, to support matched SNP6 data. 

+### If null, use coverage files or log.ratio to segment the data.  

 seg.file.sdev=0.4,

-### If seg.file provided, the log-ratio standard deviation, used to model likelihood of sub-clonal copy number events.

+### If seg.file provided, the log-ratio standard deviation, 

+### used to model likelihood of sub-clonal copy number events.

 vcf.file=NULL, 

-### VCF file, tested with MuTect output files.  Optional, but typically needed to select between local optima of similar likelihood. Can also be a CollapsedVCF, read with the readVcf function. Requires a DB info flag for dbSNP membership. The default fun.setPriorVcf function will also look for a Cosmic.CNT slot, containing the hits in the COSMIC database. Again, do not expect very useful results without a VCF file.

+### VCF file, tested with MuTect output files.  Optional, but

+### typically needed to select between local optima of similar likelihood. Can 

+### also be a CollapsedVCF, read with the readVcf function. Requires a DB info 

+### flag for dbSNP membership. The default fun.setPriorVcf function will also

+### look for a Cosmic.CNT slot, containing the hits in the COSMIC database.

+### Again, do not expect very useful results without a VCF file.

 genome="hg19",

 ### Genome version, required for the readVcf function.

+sex=c("?","F","M"),

+### Sex of sample. If ?, detect.

 fun.filterVcf=filterVcfMuTect, 

-### Function for filtering variants. Expected output is a list with elements vcf (CollapsedVCF), flag (TRUE/FALSE) and flag_comment (string). The flags will be added to the output data and can be used to warn users, for example when samples look too noisy. Default filter will remove variants flagged by MuTect, but will keep germline variants. If ran in matched normal mode, it will by default use somatic status of variants and filter non-somatic calls with allelic fraction significantly different from 0.5 in normal. 

+### Function for filtering variants. Expected output is a 

+### list with elements vcf (CollapsedVCF), flag (TRUE/FALSE) and flag_comment 

+### (string). The flags will be added to the output data and can be used to 

+### warn users, for example when samples look too noisy. Default filter will 

+### remove variants flagged by MuTect, but will keep germline variants. If 

+### ran in matched normal mode, it will by default use somatic status of 

+### variants and filter non-somatic calls with allelic fraction significantly 

+### different from 0.5 in normal. 

 args.filterVcf=list(),

-### Arguments for variant filtering function. Arguments vcf, tumor.id.in.vcf, coverage.cutoff and verbose are required in the filter function and are automatically set (do NOT set them here again).

+### Arguments for variant filtering function. Arguments 

+### vcf, tumor.id.in.vcf, coverage.cutoff and verbose are required in the 

+### filter function and are automatically set (do NOT set them here again).

 fun.setPriorVcf=setPriorVcf,

-### Function to set prior for somatic status for each variant in the VCF.

+### Function to set prior for somatic status for each 

+### variant in the VCF.

 args.setPriorVcf=list(),

 ### Arguments for somatic prior function.

 fun.segmentation=segmentationCBS, 

-### Function for segmenting the copy number log-ratios. Expected return value is a list with elements seg (the segmentation) and size (the size in bp for all segments).

+### Function for segmenting the copy number log-ratios. 

+### Expected return value is a list with elements seg (the segmentation) and 

+### size (the size in bp for all segments).

 args.segmentation=list(),

-### Arguments for segmentation function. Arguments normal, tumor, log.ratio, plot.cnv, coverage.cutoff, sampleid, vcf, tumor.id.in.vcf, verbose are required in the segmentation function and automatically set (do NOT set them here again).

+### Arguments for segmentation function. Arguments 

+### normal, tumor, log.ratio, plot.cnv, coverage.cutoff, sampleid, vcf, 

+### tumor.id.in.vcf, verbose are required in the segmentation function and 

+### automatically set (do NOT set them here again).

 fun.focal=findFocal,

 ### Function for identifying focal amplifications.

 args.focal=list(),

@@ -41,37 +73,72 @@ min.ploidy=1,

 max.ploidy=6, 

 ### Maximum ploidy to be considered.

 test.num.copy=0:7, 

-### Copy numbers tested in the grid search. Note that focal amplifications can have much higher copy numbers, but they will be labeled as subclonal (because they do not fit the integer copy numbers).

+### Copy numbers tested in the grid search. Note that focal 

+### amplifications can have much higher copy numbers, but they will be labeled

+### as subclonal (because they do not fit the integer copy numbers).

 test.purity=seq(0.05,0.95,by=0.01), 

 ### Considered tumor purity values. 

 prior.purity=rep(1,length(test.purity))/length(test.purity), 

-### Priors for purity if they are available. Only change when you know what you are doing.

+### Priors for purity if they are available. Only change 

+### when you know what you are doing.

 max.candidate.solutions=15, 

-### Number of local optima considered in optimization and variant fitting steps. If there are too many local optima, it will use specified number of top candidate solutions, but will also include all optima close to diploid, because silent genomes have often lots of local optima.

+### Number of local optima considered in optimization 

+### and variant fitting steps. If there are too many local optima, it will use

+### specified number of top candidate solutions, but will also include all 

+### optima close to diploid, because silent genomes have often lots of local 

+### optima.

 candidates=NULL, 

-### Candidates to optimize from a previous run (return.object$candidates). If NULL, do 2D grid search and find local optima. 

+### Candidates to optimize from a previous run 

+### (return.object$candidates). 

+### If NULL, do 2D grid search and find local optima. 

 coverage.cutoff=15, 

-### Minimum exon coverage in both normal and tumor. Exons with lower coverage are ingored. The cutoff choice depends on the expected purity and overall coverage. High purity samples might need a lower cutoff to call homozygous deletions. If an exon.weigh.file (below) is NOT specified, it is recommended to set a higher cutoff (e.g. 20) to remove noise from unreliable exon measurements. 

+### Minimum exon coverage in both normal and tumor. Exons

+### with lower coverage are ingored. The cutoff choice depends on the expected

+### purity and overall coverage. High purity samples might need a lower cutoff

+### to call homozygous deletions. If an exon.weigh.file (below) is NOT 

+### specified, it is recommended to set a higher cutoff (e.g. 20) to remove 

+### noise from unreliable exon measurements. 

 max.non.clonal=0.2, 

-### Maximum genomic fraction assigned to a subclonal copy number state.

+### Maximum genomic fraction assigned to a subclonal copy 

+### number state.

+max.homozygous.loss=0.1,

+### Maximum genomic fraction assigned to homozygous loss. 

+### This is set to a fairly high default value to not exclude correct

+### solutions, especially in noisy segmentations. 

 iterations=30, 

-### Maximum number of iterations in the Simulated Annealing copy number fit optimization.

+### Maximum number of iterations in the Simulated Annealing copy 

+### number fit optimization.

 log.ratio.calibration=0.25,

-### re-calibrate log-ratios in the window sd(log.ratio)*log.ratio.calibration.

+### re-calibrate log-ratios in the window 

+### sd(log.ratio)*log.ratio.calibration.

 gc.gene.file=NULL, 

-### GC gene file, which assigns GC content and gene symbols to each exon in the coverage files. Used for generating gene level calls. First column in format CHR:START-END. Second column GC content (0 to 1). Third column gene symbol.

+### GC gene file, which assigns GC content and gene symbols 

+### to each exon in the coverage files. Used for generating gene level calls. 

+### First column in format CHR:START-END. Second column GC content (0 to 1). 

+### Third column gene symbol.

 filter.lowhigh.gc.exons=0.001,

-### Quantile q (defines lower q and upper 1-q) for removing exons with outlier GC profile. Assuming that GC correction might not have been worked on those. Requires gc.gene.file.

+### Quantile q (defines lower q and upper 1-q) 

+### for removing exons with outlier GC profile. Assuming that GC correction 

+### might not have been worked on those. Requires gc.gene.file.

 filter.targeted.base=4,

-### Exclude exons with targeted base (size) smaller than this cutoff.

+### Exclude exons with targeted base (size) smaller 

+### than this cutoff.

 max.logr.sdev=0.75,

-### Flag noisy samples with segment log-ratio standard deviation larger than this. Assay specific and needs to be calibrated.

+### Flag noisy samples with segment log-ratio standard deviation 

+### larger than this. Assay specific and needs to be calibrated.

+max.segments=200,

+### Flag noisy samples with a large number of segments. Assay 

+### specific and needs to be calibrated.

 plot.cnv=TRUE, 

 ### Generate segmentation plots.

 verbose=TRUE, 

 ### Verbose output.

 post.optimize=FALSE,

-### Optimize purity using final SCNA-fit and SNVs. This might take a long time when lots of SNVs need to be fitted, but will result in a slightly more accurate purity, especially for rather silent genomes or very low purities. Otherwise, it will just use the purity determined via the SCNA-fit.

+### Optimize purity using final SCNA-fit and SNVs. This 

+### might take a long time when lots of SNVs need to be fitted, but will 

+### typically result in a slightly more accurate purity, especially for rather

+### silent genomes or very low purities. Otherwise, it will just use the 

+### purity determined via the SCNA-fit.

 ... 

 ### Additional parameters passed to the segmentation function.

 ) {

@@ -132,6 +199,29 @@ post.optimize=FALSE,

         if (is.null(gatk.normal.file)) normal <- tumor

         if (!is.null(seg.file)) stop("Provide either log.ratio or seg.file, not both.") 

     }        

+    sex <- match.arg(sex)

+    sex.chr <- .getSexChr(tumor)

+    if (sex =="?") {

+        sex.tumor <- getSexFromCoverage(tumor, verbose=FALSE)

+        sex.normal <- getSexFromCoverage(normal, verbose=FALSE)

+        if (is.na(sex.tumor)) {

+            tumor <- .removeChr(tumor, remove.chrs=sex.chr)

+        }    

+        if (is.na(sex.normal)) {

+            normal <- .removeChr(normal, remove.chrs=sex.chr)

+        }    

+        if (!identical(sex.tumor, sex.normal)) warning(paste("Sex tumor/normal mismatch: tumor =", sex.tumor, "normal =", sex.normal))  

+        sex <- sex.tumor    

+        if (is.na(sex)) sex = "?"

+    } 

+    if (sex=="M") {

+        tumor <- .removeChr(tumor, remove.chrs=sex.chr)

+    } else if (sex=="F") {

+        tumor <- .removeChr(tumor, remove.chrs=sex.chr[2])

+    }       

+          

+    if (verbose) message(paste("Sex of sample:", sex))

+          

     # NA's in log.ratio confuse the CBS function

     idx <- !is.na(log.ratio) & !is.infinite(log.ratio)

     log.ratio <- log.ratio[idx]    

@@ -171,6 +261,7 @@ post.optimize=FALSE,

     vcf.germline <- NULL

     tumor.id.in.vcf <- NULL

     prior.somatic <- NULL

+    vcf.filtering <-list(flag=FALSE, flag_comment="")

     if (!is.null(vcf.file)) {

         if (verbose) message("Loading VCF...")

@@ -358,6 +449,10 @@ post.optimize=FALSE,

                 # set probability to zero if ploidy is not within requested range 

                 log.prior.ploidy <- log(ifelse(ploidy<min.ploidy | ploidy>max.ploidy,0,1))

                 if (iter > 1) p.rij <- p.rij+log.prior.ploidy

+                

+                frac.homozygous.loss <- sapply(test.num.copy, function(Ci) (sum(li[-i] * ifelse(C[-i]==0,1,0)) + li[i] * ifelse(Ci==0,1,0))/sum(li) )

+                log.prior.homozygous.loss <- log(ifelse(frac.homozygous.loss > max.homozygous.loss,0,1))

+                if (iter > 1) p.rij <- p.rij+log.prior.homozygous.loss

                 # model sub clonal state with a uniform distribution

                 p.rij <- c(p.rij,  .calcLlikSegmentSubClonal( exon.lrs[[i]]+log.ratio.offset[i], max.exon.ratio))

@@ -443,7 +538,7 @@ post.optimize=FALSE,

             SNV.posterior <- res.snvllik[[idx]]

         }

-        list(log.likelihood=llik, purity=p, ploidy=weighted.mean(C,li), total.ploidy=total.ploidy, seg=seg.adjusted, C.posterior=data.frame(C.posterior/rowSums(C.posterior), ML.C=C, ML.Subclonal=subclonal), SNV.posterior=SNV.posterior, fraction.subclonal=subclonal.f, gene.calls=gene.calls, log.ratio.offset=log.ratio.offset, failed=FALSE)

+        list(log.likelihood=llik, purity=p, ploidy=weighted.mean(C,li), total.ploidy=total.ploidy, seg=seg.adjusted, C.posterior=data.frame(C.posterior/rowSums(C.posterior), ML.C=C, ML.Subclonal=subclonal), SNV.posterior=SNV.posterior, fraction.subclonal=subclonal.f, fraction.homozygous.loss=sum(li[which(C<0.01)])/sum(li), gene.calls=gene.calls, log.ratio.offset=log.ratio.offset, failed=FALSE)

     }

     results <- lapply(1:nrow(candidate.solutions$candidates),.optimizeSolution)

@@ -451,7 +546,7 @@ post.optimize=FALSE,

     results <- .rankResults(results)

     results <- .filterDuplicatedResults(results)

-    results <- .flagResults(results, max.non.clonal=max.non.clonal, max.logr.sdev=max.logr.sdev, logr.sdev=sd.seg, flag=vcf.filtering$flag, flag_comment=vcf.filtering$flag_comment)  

+    results <- .flagResults(results, max.non.clonal=max.non.clonal, max.logr.sdev=max.logr.sdev, logr.sdev=sd.seg, max.segments=max.segments, flag=vcf.filtering$flag, flag_comment=vcf.filtering$flag_comment)  

     if (!is.null(gc.gene.file)) {

         # Add LOH calls to gene level calls 

         for (i in 1:length(results)) {

@@ -459,30 +554,37 @@ post.optimize=FALSE,

         }

     }    

-    for (i in 1:length(results)) {

-        results[[i]]$SNV.posterior$beta.model$posteriors <- .getVariantCallsLOH( results[[i]] )

+    if (!is.null(vcf.file)) {

+        for (i in 1:length(results)) {

+            results[[i]]$SNV.posterior$beta.model$posteriors <- .getVariantCallsLOH( results[[i]] )

+        }

     }

     ##value<< A list with elements

     list(

         candidates=candidate.solutions, ##<< Results of the grid search.

         results=results, ##<< All local optima, sorted by final rank.

-        input=list(tumor=gatk.tumor.file, normal=gatk.normal.file, log.ratio=data.frame(probe=normal[,1], log.ratio=log.ratio), log.ratio.sdev=sd.seg, vcf=vcf, sampleid=sampleid) ##<< The input data.

+        input=list(tumor=gatk.tumor.file, normal=gatk.normal.file, log.ratio=data.frame(probe=normal[,1], log.ratio=log.ratio), log.ratio.sdev=sd.seg, vcf=vcf, sampleid=sampleid, sex=sex) ##<< The input data.

         )

 ##end<<

 },ex=function(){

-gatk.normal.file <- system.file("extdata", "example_normal.txt", package="PureCN")

-gatk.tumor.file <- system.file("extdata", "example_tumor.txt", package="PureCN")

-vcf.file <- system.file("extdata", "example_vcf.vcf", package="PureCN")

-gc.gene.file <- system.file("extdata", "example_gc.gene.file.txt", package="PureCN")

-

-# Speed-up the runAbsoluteCN by using the stored grid-search 

+gatk.normal.file <- system.file("extdata", "example_normal.txt", 

+    package="PureCN")

+gatk.tumor.file <- system.file("extdata", "example_tumor.txt", 

+    package="PureCN")

+vcf.file <- system.file("extdata", "example_vcf.vcf", 

+    package="PureCN")

+gc.gene.file <- system.file("extdata", "example_gc.gene.file.txt", 

+    package="PureCN")

+

+# Speed-up the runAbsoluteCN call by using the stored grid-search 

 # (purecn.example.output$candidates).

 data(purecn.example.output)

 # The max.candidate.solutions parameter is set to a very low value only to

 # speed-up this example.  This is not a good idea for real samples.

-ret <-runAbsoluteCN(gatk.normal.file=gatk.normal.file, gatk.tumor.file=gatk.tumor.file, 

-   candidates=purecn.example.output$candidates, max.candidate.solutions=2,

+ret <-runAbsoluteCN(gatk.normal.file=gatk.normal.file, 

+    gatk.tumor.file=gatk.tumor.file, 

+    candidates=purecn.example.output$candidates, max.candidate.solutions=2,

     vcf.file=vcf.file, sampleid='Sample1', gc.gene.file=gc.gene.file)

 })    

@@ -69,7 +69,7 @@ ret <-runAbsoluteCN(gatk.normal.file=gatk.normal.file, gatk.tumor.file=gatk.tumo

 .pruneByVCF <- function(x, vcf, tumor.id.in.vcf, min.size=5, max.pval=0.00001, iterations=3, debug=FALSE) {

     seg <- segments.p(x$cna)

     for (iter in 1:iterations) {

-        seg.gr <- GRanges(seqnames=paste("chr", seg$chrom,sep=""), IRanges(start=seg$loc.start, end=seg$loc.end))

+        seg.gr <- GRanges(seqnames=.add.chr.name(seg$chrom), IRanges(start=seg$loc.start, end=seg$loc.end))

         ov <- findOverlaps(seg.gr, vcf)

         ar <- sapply(geno(vcf)$FA[,tumor.id.in.vcf], function(x) x[1])

         ar.r <- ifelse(ar>0.5, 1-ar, ar)

deleted file mode 100755

@@ -1,119 +0,0 @@

-<h1 id="quick-start">Quick Start</h1>

-<h2 id="basic-input-files">Basic Input Files</h2>

-<p>You will need coverage files in GATK DepthOfCoverage format for tumor and a control sample:</p>

-<pre><code>Target  total_coverage  average_coverage    Sample1_total_cvg   Sample1_mean_cvg    Sample1_granular_Q1 Sample1_granular_median Sample1_granular_Q3 Sample1_._above_15

-chr1:69091-70009    0   0   0   0   1   1   1   0

-chr1:367659-368598  6358    9.25121153819759    6358    9.25121153819759    1   7   13  11.6

-chr1:621096-622035  6294    9.16910019318401    6294    9.16910019318401    1   7   12  9.5</code></pre>

-<p>PureCN will look for the columns &quot;Target&quot;, &quot;total_coverage&quot;, and &quot;average_coverage&quot;. The algorithm works best when the coverage data is already GC normalized. The steps below will NOT GC normalize.</p>

-<p>You will also need a VCF file with both somatic and germline variants, for example as obtained by MuTect.</p>

-<h2 id="first-run">First Run</h2>

-<p>In R:</p>

-<pre><code>library(PureCN)

-ret <-runAbsoluteCN(gatk.normal.file=gatk.normal.file, 

-                    gatk.tumor.file=gatk.tumor.file, 

-                    vcf.file=vcf.file)</code></pre>

-<p>To obtain gene level copy number calls, we need an exon-to-gene mapping file. This is provided together with GC content via the gc.gene.file argument. The expected format:</p>

-<pre><code>Target  gc_bias Gene

-chr1:69091-70009    0.427638737758433   OR4F5

-chr1:367659-368598  0.459574468085106   OR4F29

-chr1:621096-622035  0.459574468085106   OR4F3</code></pre>

-<h2 id="pool-of-normals">Pool of Normals</h2>

-<p>If a pool of normal samples is available, we can use the coverage data to estimate the expected variance in coverage per exon. This will help the segmentation algorithm to filter noise:</p>

-<pre><code>gatk.normal.files &lt;- dir(&quot;poolofnormals&quot;, 

-    pattern="coverage.sample_interval_summary", full.names=TRUE) 

-createExonWeightFile(gatk.normal.files[1:2], gatk.normal.files[-(1:2)], &quot;exon_weights.txt&quot;)</code></pre>

-<p>This will calculate log-ratios for all normals using the first two as tumor and the remaining ones as normal. This estimates then exon level coverage standard deviations.</p>

-<p>We can also use the pool of normal to find SNPs with biased allelic fractions (significantly different from 0.5 for heterozygous SNPs):</p>

-<pre><code>mutect.normal.files &lt;- dir(&quot;poolofnormals&quot;, pattern=&quot;vcf$&quot;, full.names=TRUE) 

-snp.blacklist <- createSNPBlacklist(mutect.normal.files)

-write.csv(snp.bl, file="SNP_blacklist.csv")

-write.csv(snp.bl, file=&quot;SNP_blacklist_segmented.csv&quot;, row.names=FALSE, quote=FALSE)</code></pre>

-<p>Finally, we can use the pool of normals to find normal samples for log-ratio calculation, especially when no matched normal sample is available:</p>

-<pre><code>if (!file.exists(&quot;normalDB.rds&quot;)) {

-    normalDB <- createNormalDatabase(gatk.normal.files)

-    saveRDS(normalDB, "normalDB.rds")

-} else {

-    normalDB <- readRDS("normalDBe4.rds")

-}

-

-# get the best 4 normals and average them

-gatk.normal.file <- findBestNormal(gatk.tumor.file, normalDB, num.normals=4)

-pool <- pool.coverage2(lapply(gatk.normal.file, read.coverage.gatk), remove.chrs=c('chrX', 'chrY'))

-

-# get the best normal

-gatk.normal.file &lt;- findBestNormal(gatk.tumor.file, normalDB)</code></pre>

-<h2 id="vcf-artifact-filtering">VCF Artifact Filtering</h2>

-<p>If MuTect was run in matched normal mode, then all germline variants are rejected, that means we cannot just filter by the KEEP/REJECT MuTect flags. The filterVcfMuTect function optionally reads the MuTect stats file and will keep germline variants, while removing potential artifacts. Otherwise it will use only the filters based on read depths as defined in filterVcfBasic.</p>

-<h2 id="recommended-run">Recommended Run</h2>

-<p>Finally, we can run PureCN with all that information:</p>

-<pre><code>  ret &lt;-runAbsoluteCN(gatk.normal.file=pool, gatk.tumor.file=gatk.tumor.file, 

-    vcf.file=vcf.file, sampleid='Sample1', gc.gene.file=gc.gene.file, 

-    args.filterVcf=list(snp.blacklist=snp.blacklist, stats.file=mutect.stats.file), 

-    args.segmentation=list(exon.weight.file=exon.weight.file), 

-    post.optimize=TRUE)</code></pre>

-<p>The post.optimize flag will increase the runtime by a lot, but might be worth it if the copy number profile is not very clean. This is recommended with lower coverage datasets or if there are significant capture biases.</p>

-<p>We also need to create a few output files:</p>

-<pre><code>file.rds &lt;- &#39;Sample1_PureCN.rds&#39;

-saveRDS(ret, file=file.rds)

-pdf('Sample1_PureCN.pdf'), width=10, height=12)

-plotAbs(ret, type='all')

-dev.off()</code></pre>

-<h2 id="manual-curation">Manual Curation</h2>

-<p>For prediction of SNV status (germline vs somatic, sub-clonal vs. clonal, homozyous vs heterozygous), it is important that both purity and ploidy are correct. We provide functionality for curating results:</p>

-<pre><code>createCurationFile(file.rds) </code></pre>

-<p>This will generate a CSV file, in which the correct purity and ploidy values can be manually entered. It also contains a column &quot;Curated&quot;, which should be set to TRUE, otherwise the file will be overwritten when re-run.</p>

-<p>Then in R, the correct solution (closest to the combination in the CSV file) can be loaded with the readCurationFile function:</p>

-<pre><code>ret &lt;- readCurationFile(file.rds)</code></pre>

-<h2 id="custom-segmentation">Custom Segmentation</h2>

-<p>By default, we will use DNAcopy to segment the log-ratio. You can easily change that to your favorite method and the segmentationCBS function can serve as an example.</p>

-<pre><code>segmentationCBS &lt;- function (normal, tumor, log.ratio, plot.cnv, coverage.cutoff, 

-    sampleid = sampleid, exon.weight.file = NULL, alpha = 0.005, 

-    vcf = NULL, tumor.id.in.vcf = 1, verbose = TRUE, ...) </code></pre>

-<p>It is also possible to provide segmented file, which we however only recommend when matched SNP6 data is available (otherwise it is better to customize the segmentation function as described above). The expected file format is:</p>

-<pre><code>ID  chrom   loc.start   loc.end num.mark    seg.mean

-Sample1   1   61723   5773942 2681    0.125406444072723

-Sample1   1   5774674 5785170 10  -0.756511807441712</code></pre>

-<h2 id="output">Output</h2>

-<p>The plotAbs() call above will generate the main plots shown in the manuscript. The R data file (file.rds) contains gene level copy number calls, SNV status and LOH calls. The purity/ploidy combinations are sorted by likelihood at stored in ret$results.</p>

-<pre><code>names(ret)

-head(ret$results[[1]]$gene.calls, 3)

-chr    start      end C seg.mean seg.id number.exons gene.mean   gene.min

-TP73 chr1  3598833  3649721 3   0.3434      1           14 0.1755614 -0.3861084

-CHD5 chr1  6161904  6240126 3   0.3434      1           44 0.3059672 -0.8812580

-MTOR chr1 11166594 11319542 3   0.3434      1           62 0.2360036 -0.1346406

-gene.max focal num.snps.loh.segment percentage.loh.in.loh.segment

-TP73 0.6437573 FALSE                    0                             0

-CHD5 1.5984124 FALSE                   49                             0

-MTOR 0.6705342 FALSE                   49                             0</code></pre>

-<p>This data.frame also contains gene level LOH information. The SNV posteriors:</p>

-<pre><code>head(ret$results[[1]]$SNV.posterior$beta.model$posteriors, 3)

-           seqnames   start     end SOMATIC.M0   SOMATIC.M1   SOMATIC.M2

-rs6696489      chr1 6162054 6162054          0 9.155037e-04 2.938312e-07

-rs59788818     chr1 6171992 6171992          0 2.364795e-05 3.397817e-07

-rs2843493      chr1 6184092 6184092          0 9.362524e-16 2.714454e-04

-             SOMATIC.M3   SOMATIC.M4    SOMATIC.M5    SOMATIC.M6 SOMATIC.M7

-rs6696489  6.592308e-21 2.743271e-76 2.002534e-167 3.562023e-275          0

-rs59788818 3.491208e-34 1.711087e-93 7.197813e-188 2.275374e-298          0

-rs2843493  1.037518e-11 5.167773e-67 3.034251e-158 3.725864e-266          0

-            GERMLINE.M0  GERMLINE.M1  GERMLINE.M2  GERMLINE.M3   GERMLINE.M4

-rs6696489  1.723074e-04 9.989057e-01 6.216658e-06 3.302211e-27  1.067221e-82

-rs59788818 5.836546e-23 9.999757e-01 2.683928e-07 1.233256e-51 2.340805e-111

-rs2843493  8.802336e-49 1.023147e-07 9.997285e-01 2.798399e-18  3.718923e-74

-             GERMLINE.M5   GERMLINE.M6 GERMLINE.M7 GERMLINE.CONTHIGH

-rs6696489  6.678185e-174 1.070845e-281           0      3.060961e-62

-rs59788818 5.523414e-206 1.183050e-316           0     8.253575e-127

-rs2843493  9.760218e-166 6.969108e-274           0      5.841728e-59

-           GERMLINE.CONTLOW ML.SOMATIC ML.M ML.C     ML.AR        AR

-rs6696489      2.105422e-22      FALSE    1    3 0.3571429 0.2803790

-rs59788818     3.106132e-86      FALSE    1    3 0.3571429 0.3957373

-rs2843493     4.644636e-138      FALSE    2    3 0.6428571 0.6259540

-           CN.Subclonal Log.Ratio Prior.Somatic Prior.Contamination ML.LOH

-rs6696489         FALSE 1.5984124  0.0004950495                0.01  FALSE

-rs59788818        FALSE 0.2600908  0.0004950495                0.01  FALSE

-rs2843493         FALSE 0.7559721  0.0004950495                0.01  FALSE

-           num.snps.loh.segment percentage.loh.in.loh.segment

-rs6696489                    49                             0

-rs59788818                   49                             0

-rs2843493                    49                             0</code></pre>

-<p>This lists all posterior probabilities for all possible SNV states.</p>

deleted file mode 100755

@@ -1,191 +0,0 @@

-# Quick Start

-

-## Basic Input Files 

-

-You will need coverage files in GATK DepthOfCoverage format for tumor and a control sample:

-

-    Target  total_coverage  average_coverage    Sample1_total_cvg   Sample1_mean_cvg    Sample1_granular_Q1 Sample1_granular_median Sample1_granular_Q3 Sample1_._above_15

-    chr1:69091-70009    0   0   0   0   1   1   1   0

-    chr1:367659-368598  6358    9.25121153819759    6358    9.25121153819759    1   7   13  11.6

-    chr1:621096-622035  6294    9.16910019318401    6294    9.16910019318401    1   7   12  9.5

-

-PureCN will look for the columns "Target", "total_coverage", and

-"average_coverage". The algorithm works best when the coverage data is already

-GC normalized. The steps below will NOT GC normalize.

-

-You will also need a VCF file with both somatic and germline variants, for

-example as obtained by MuTect.

-

-## First Run

-

-In R:

-

-    library(PureCN)

-    ret <-runAbsoluteCN(gatk.normal.file=gatk.normal.file, 

-                        gatk.tumor.file=gatk.tumor.file, 

-                        vcf.file=vcf.file)

-

-To obtain gene level copy number calls, we need an exon-to-gene mapping file.

-This is provided together with GC content via the gc.gene.file argument. The

-expected format:

-

-    Target	gc_bias	Gene

-    chr1:69091-70009	0.427638737758433	OR4F5

-    chr1:367659-368598	0.459574468085106	OR4F29

-    chr1:621096-622035	0.459574468085106	OR4F3

-

-## Pool of Normals

-

-If a pool of normal samples is available, we can use the coverage data to

-estimate the expected variance in coverage per exon. This will help the

-segmentation algorithm to filter noise:

-

-    gatk.normal.files <- dir("poolofnormals", 

-        pattern="coverage.sample_interval_summary", full.names=TRUE) 

-    createExonWeightFile(gatk.normal.files[1:2], gatk.normal.files[-(1:2)], "exon_weights.txt")

-

-This will calculate log-ratios for all normals using the first two as tumor and

-the remaining ones as normal. This estimates then exon level coverage standard

-deviations.

-

-We can also use the pool of normal to find SNPs with biased allelic fractions

-(significantly different from 0.5 for heterozygous SNPs):

-

-    mutect.normal.files <- dir("poolofnormals", pattern="vcf$", full.names=TRUE) 

-    snp.blacklist <- createSNPBlacklist(mutect.normal.files)

-    write.csv(snp.bl, file="SNP_blacklist.csv")

-    write.csv(snp.bl, file="SNP_blacklist_segmented.csv", row.names=FALSE, quote=FALSE)

-

-Finally, we can use the pool of normals to find normal samples for log-ratio

-calculation, especially when no matched normal sample is available:

-

-    if (!file.exists("normalDB.rds")) {

-        normalDB <- createNormalDatabase(gatk.normal.files)

-        saveRDS(normalDB, "normalDB.rds")

-    } else {

-        normalDB <- readRDS("normalDBe4.rds")

-    }

-    

-    # get the best 4 normals and average them

-    gatk.normal.file <- findBestNormal(gatk.tumor.file, normalDB, num.normals=4)

-    pool <- poolCoverage(lapply(gatk.normal.file, read.coverage.gatk), remove.chrs=c('chrX', 'chrY'))

-    

-    # get the best normal

-    gatk.normal.file <- findBestNormal(gatk.tumor.file, normalDB)

-

-## VCF Artifact Filtering

-

-If MuTect was run in matched normal mode, then all germline variants are

-rejected, that means we cannot just filter by the PASS/REJECT MuTect flags. The

-filterVcfMuTect function optionally reads the MuTect stats file and will keep

-germline variants, while removing potential artifacts. Otherwise it will use

-only the filters based on read depths as defined in filterVcfBasic.

-

-## Recommended Run

-

-Finally, we can run PureCN with all that information:

-

-      ret <-runAbsoluteCN(gatk.normal.file=pool, gatk.tumor.file=gatk.tumor.file, 

-        vcf.file=vcf.file, sampleid='Sample1', gc.gene.file=gc.gene.file, 

-        args.filterVcf=list(snp.blacklist=snp.blacklist, stats.file=mutect.stats.file), 

-        args.segmentation=list(exon.weight.file=exon.weight.file), 

-        post.optimize=TRUE)

-

-The post.optimize flag will increase the runtime by a lot, but might be worth

-it if the copy number profile is not very clean. This is recommended with lower

-coverage datasets or if there are significant capture biases. For high quality

-whole exome data, this is typically not necessary.

-

-We also need to create a few output files:

-

-    file.rds <- 'Sample1_PureCN.rds'

-    saveRDS(ret, file=file.rds)

-    pdf('Sample1_PureCN.pdf', width=10, height=12)

-    plotAbs(ret, type='all')

-    dev.off()

-

-## Manual Curation

-

-For prediction of SNV status (germline vs somatic, sub-clonal vs. clonal,

-homozyous vs heterozygous), it is important that both purity and ploidy are

-correct. We provide functionality for curating results:

-

-    createCurationFile(file.rds) 

-

-This will generate a CSV file, in which the correct purity and ploidy values

-can be manually entered. It also contains a column "Curated", which should be

-set to TRUE, otherwise the file will be overwritten when re-run.

-

-Then in R, the correct solution (closest to the combination in the CSV file)

-can be loaded with the readCurationFile function:

-

-    ret <- readCurationFile(file.rds)

-    

-## Custom Segmentation

-

-By default, we will use DNAcopy to segment the log-ratio. You can easily change

-that to your favorite method and the segmentationCBS function can serve as an

-example.

-

-    segmentationCBS <- function (normal, tumor, log.ratio, plot.cnv, coverage.cutoff, 

-        sampleid = sampleid, exon.weight.file = NULL, alpha = 0.005, 

-        vcf = NULL, tumor.id.in.vcf = 1, verbose = TRUE, ...) 

-

-It is also possible to provide segmented file, which we however only recommend

-when matched SNP6 data is available (otherwise it is better to customize the

-segmentation function as described above). The expected file format is:

-

-    ID  chrom   loc.start   loc.end num.mark    seg.mean

-    Sample1   1   61723   5773942 2681    0.125406444072723

-    Sample1   1   5774674 5785170 10  -0.756511807441712

-

-## Output

-

-The plotAbs() call above will generate the main plots shown in the manuscript. The R data file (file.rds) contains gene level copy number calls, SNV status and LOH calls.

-The purity/ploidy combinations are sorted by likelihood at stored in ret$results.

-

-    names(ret)

-    head(ret$results[[1]]$gene.calls, 3)

-    chr    start      end C seg.mean seg.id number.exons gene.mean   gene.min

-    TP73 chr1  3598833  3649721 3   0.3434      1           14 0.1755614 -0.3861084

-    CHD5 chr1  6161904  6240126 3   0.3434      1           44 0.3059672 -0.8812580

-    MTOR chr1 11166594 11319542 3   0.3434      1           62 0.2360036 -0.1346406

-    gene.max focal num.snps.loh.segment percentage.loh.in.loh.segment

-    TP73 0.6437573 FALSE                    0                             0

-    CHD5 1.5984124 FALSE                   49                             0

-    MTOR 0.6705342 FALSE                   49                             0

-

-This data.frame also contains gene level LOH information. The SNV posteriors:

-

-    head(ret$results[[1]]$SNV.posterior$beta.model$posteriors, 3)

-               seqnames   start     end SOMATIC.M0   SOMATIC.M1   SOMATIC.M2

-    rs6696489      chr1 6162054 6162054          0 9.155037e-04 2.938312e-07

-    rs59788818     chr1 6171992 6171992          0 2.364795e-05 3.397817e-07

-    rs2843493      chr1 6184092 6184092          0 9.362524e-16 2.714454e-04

-                 SOMATIC.M3   SOMATIC.M4    SOMATIC.M5    SOMATIC.M6 SOMATIC.M7

-    rs6696489  6.592308e-21 2.743271e-76 2.002534e-167 3.562023e-275          0

-    rs59788818 3.491208e-34 1.711087e-93 7.197813e-188 2.275374e-298          0

-    rs2843493  1.037518e-11 5.167773e-67 3.034251e-158 3.725864e-266          0

-                GERMLINE.M0  GERMLINE.M1  GERMLINE.M2  GERMLINE.M3   GERMLINE.M4

-    rs6696489  1.723074e-04 9.989057e-01 6.216658e-06 3.302211e-27  1.067221e-82

-    rs59788818 5.836546e-23 9.999757e-01 2.683928e-07 1.233256e-51 2.340805e-111

-    rs2843493  8.802336e-49 1.023147e-07 9.997285e-01 2.798399e-18  3.718923e-74

-                 GERMLINE.M5   GERMLINE.M6 GERMLINE.M7 GERMLINE.CONTHIGH

-    rs6696489  6.678185e-174 1.070845e-281           0      3.060961e-62

-    rs59788818 5.523414e-206 1.183050e-316           0     8.253575e-127

-    rs2843493  9.760218e-166 6.969108e-274           0      5.841728e-59

-               GERMLINE.CONTLOW ML.SOMATIC ML.M ML.C     ML.AR        AR

-    rs6696489      2.105422e-22      FALSE    1    3 0.3571429 0.2803790

-    rs59788818     3.106132e-86      FALSE    1    3 0.3571429 0.3957373

-    rs2843493     4.644636e-138      FALSE    2    3 0.6428571 0.6259540

-               CN.Subclonal Log.Ratio Prior.Somatic Prior.Contamination ML.LOH

-    rs6696489         FALSE 1.5984124  0.0004950495                0.01  FALSE

-    rs59788818        FALSE 0.2600908  0.0004950495                0.01  FALSE

-    rs2843493         FALSE 0.7559721  0.0004950495                0.01  FALSE

-               num.snps.loh.segment percentage.loh.in.loh.segment

-    rs6696489                    49                             0

-    rs59788818                   49                             0

-    rs2843493                    49                             0

-

-This lists all posterior probabilities for all possible SNV states.

-

@@ -1,12 +1,14 @@

 \name{createExonWeightFile}

 \alias{createExonWeightFile}

 \title{Calculate exon weights}

-\description{Creates an exon weight file useful for segmentation, by down-weighting unreliable exons.}

+\description{Creates an exon weight file useful for segmentation, by 

+down-weighting unreliable exons.}

 \usage{createExonWeightFile(gatk.tumor.files, gatk.normal.files, 

     exon.weight.file)}

 \arguments{

   \item{gatk.tumor.files}{A small number (1-3) of GATK tumor coverage samples.}

-  \item{gatk.normal.files}{A large number of GATK normal coverage samples (>20) to estimate exon log-ratio standard deviations.}

+  \item{gatk.normal.files}{A large number of GATK normal coverage samples (>20) 

+to estimate exon log-ratio standard deviations.}

   \item{exon.weight.file}{Output filename.}

 }

@@ -1,10 +1,12 @@

 \name{createNormalDatabase}

 \alias{createNormalDatabase}

 \title{createNormalDatabase}

-\description{Function to create a database of normal samples, used to find a good match for tumor copy number normalization.}

+\description{Function to create a database of normal samples, used to find 

+a good match for tumor copy number normalization.}

 \usage{createNormalDatabase(gatk.normal.files, ...)}

 \arguments{

-  \item{gatk.normal.files}{Vector with file names pointing to GATK coverage files of normal samples. }

+  \item{gatk.normal.files}{Vector with file names pointing to GATK coverage files 

+of normal samples. }

   \item{\dots}{Arguments passed to the prcomp function.}

 }

@@ -3,15 +3,15 @@

 \title{findBestNormal}

 \description{Function to find the best matching normal for a provided tumor sample.}