ExomeCopy - gsudre/autodenovo GitHub Wiki

11/02/2017

Decided to give a try to ExomeCopy, since it's one of the tools tested in the CNVscan paper (https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-2374-2#Sec8). It seems to operate in a similar fashion to XHMM. So, I think we'll eventually need to run it on the entire cohort to get the normalizations correct. For now, let's just make sure it works by running on a trio.

Following directions from http://bioconductor.org/packages/release/bioc/vignettes/exomeCopy/inst/doc/exomeCopy.pdf in BW:

library(exomeCopy)
target.file <- "/data/NCR_SBRB/simplex/SeqCapEZ_Exome_v3.0_Design_Annotation_files/SeqCap_EZ_Exome_v3_hg19_capture_targets.bed"
sample.names <- c("child", "father", "mother")
reference.file <- '/fdb/igenomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa'
target.df <- read.delim(target.file, header = FALSE)
target <- GRanges(seqname = target.df[, 1], IRanges(start = target.df[,2] + 1, end = target.df[, 3]))
counts <- target

bam.files <- c("/data/NCR_SBRB/simplex/BAM/CLIA_400216/CLIA_400216_sorted_RG_markduplicate_recalibrated.bam","/data/NCR_SBRB/simplex/BAM/CLIA_400176/CLIA_400176_sorted_RG_markduplicate_recalibrated.bam","/data/NCR_SBRB/simplex/BAM/CLIA_400177/CLIA_400177_sorted_RG_markduplicate_recalibrated.bam")
for (i in 1:length(bam.files)) {
   mcols(counts)[sample.names[i](/gsudre/autodenovo/wiki/sample.names[i)] <- countBamInGRanges(bam.files[i],target)
}
counts$GC <- getGCcontent(target, reference.file)
counts$GC.sq <- counts$GC^2
counts$bg <- generateBackground(sample.names, counts, median)
counts$log.bg <- log(counts$bg + 0.1)
counts$width <- width(counts)