sequenza (version 3.0.0)

find.breaks: Segmentation of sequencing data using an allele-specific copy number algorithm

Description

This function uses aspcf or pcf from the package copynumber to segment depth ratio and B-allele frequency obtained from sequencing data.

Usage

find.breaks(seqz.baf, gamma = 80, kmin = 10, baf.thres = c(0, 0.5),
            verbose = FALSE, seg.algo = "aspcf", ...)
segment.breaks(seqz.tab, breaks, min.reads.baf = 1, weighted.mean = TRUE)

Arguments

seqz.baf

an seqz file containing only the heterozygous positions.

seqz.tab

a complete seqz file.

gamma, kmin, baf.thres, verbose

arguments passed to the segmentation algorithm.

breaks

breaks as output by find.breaks.

min.reads.baf

threshold on the depth of the positions included to calculate the average BAF for segment.

weighted.mean

boolean to select if the segments have to calculated using the read depth as a weights to calculate depth ratio and B-allele frequency means.

seg.algo

Selects the algorithm used for the segmentation. Available options are aspcf of pcf.

...

additional arguments passed to aspcf.

Details

copynumber is a package to perform efficient segmentation of SNP-array data. The function find.breaks uses the algorithms from the copynumber package to find break points, where the default parameters have been optimized for sequencing data, but a careful choice of an optimal gamma value is advised.

Examples

# NOT RUN {
  
# }
# NOT RUN {
data.file <-  system.file("extdata", "example.seqz.txt.gz", package = "sequenza")
# read all the chromosomes:
seqz.data  <- read.seqz(data.file)
# Gather genome wide GC-stats from raw file:
gc.stats <- gc.sample.stats(data.file)
gc.vect  <- setNames(gc.stats$raw.mean, gc.stats$gc.values)
# Read only one chromosome:
seqz.data  <- read.seqz(data.file, chr.name = 12)

# Correct the coverage of the loaded chromosome:
seqz.data$adjusted.ratio <- seqz.data$depth.ratio /
                           gc.vect[as.character(seqz.data$GC.percent)]
# Select the heterozygous positions
seqz.hom  <- seqz.data$zygosity.normal == 'hom'
seqz.het  <- seqz.data[!seqz.hom, ]
# Detect breakpoints
breaks <- find.breaks(seqz.het, gamma = 80, kmin = 10, baf.thres = c(0, 0.5))
# use heterozygous and homozygous position to measure segment values
segment.breaks(seqz.data, breaks = breaks)
   
# }