R/findFocal.R

In PureCN: Copy number calling and SNV classification using targeted short read sequencing

#' Find focal amplifications
#' 
#' Function to find focal amplifications in segmented data.  This is
#' automatically called in \code{\link{runAbsoluteCN}}.
#' 
#' 
#' @param seg Segmentation data.
#' @param max.size Cutoff for focal in base pairs.
#' @param cn.diff Minimum copy number delta between neighboring segments.
#' @param min.amp.cn Minimum amplification integer copy number. Segments with
#' lower copy number are not tested.
#' @return \code{logical(n)}, indicating for all n segments whether they are
#' focally amplified or not.
#' @author Markus Riester
#' @seealso \code{\link{runAbsoluteCN}}
#' @examples
#' 
#' normal.coverage.file <- system.file("extdata", "example_normal_tiny.txt", 
#'     package="PureCN")
#' tumor.coverage.file <- system.file("extdata", "example_tumor_tiny.txt", 
#'     package="PureCN")
#' vcf.file <- system.file("extdata", "example.vcf.gz", 
#'     package="PureCN")
#' interval.file <- system.file("extdata", "example_intervals_tiny.txt", 
#'     package="PureCN")
#' 
#' # The max.candidate.solutions, max.ploidy and test.purity parameters are set to
#' # non-default values to speed-up this example.  This is not a good idea for real
#' # samples.
#' ret <-runAbsoluteCN(normal.coverage.file=normal.coverage.file, 
#'     tumor.coverage.file=tumor.coverage.file, vcf.file=vcf.file, genome="hg19", 
#'     sampleid="Sample1", interval.file=interval.file,
#'     max.candidate.solutions=1, max.ploidy=4, test.purity=seq(0.3,0.7,by=0.05), 
#'     args.focal=list(max.size = 2e+06), fun.focal=findFocal)
#' 
#' @export findFocal
findFocal <- function(seg, max.size = 3000000, cn.diff = 2, min.amp.cn = 5) {
    focal <- rep(FALSE, nrow(seg))
    for (i in seq_len(nrow(seg))) {
        if (seg$C[i] < min.amp.cn) next
        if (seg$size[i] > max.size) next
        size <- seg$size[i]
        if (i>1) {
            for (j in (i-1):1) {
                if (seg$C[j] < seg$C[i]-cn.diff) {
                    break
                }
                size <- size + seg$size[j]    
            }        
        }    
        if (i<nrow(seg)) {
            for (j in (i+1):nrow(seg)) {
                if (seg$C[j] < seg$C[i]-cn.diff) {
                    break
                }
                size <- size + seg$size[j]    
            }        
        }    
        focal[i] <- size < max.size        
    }
    focal
}