Dear Gunter, and all,

Please may I ask for an advice regarding the use of cn.MOPS on WES data containing TUMOR-GERMLINE paired samples. I wrote the following piece of R code (below) -- please would you let me know if it looks fine, or if it needs any modifications.

#####

segments <- read.table("hg38_refFlat.longest.isoform",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))

TUMOR <- "zzz2.bam"       
GERMLINE <- "zzz1.bam"

tumor <- getSegmentReadCountsFromBAM(TUMOR, GR=gr, parallel=1)
normal <- getSegmentReadCountsFromBAM(GERMLINE, GR=gr, parallel=1)

#####

X <- tumor
values(X) <- cbind(values(tumor),values(normal))
head(X)

#####

resCNMOPS <- exomecn.mops(X)
resCNMOPS <- calcIntegerCopyNumbers(resCNMOPS)
 

#############
 

-- bogdan

Hello Bogdan,

The code is ok until "exomecn.mops". In a tumor vs. reference setting, I recommend to use "referencecn.mops", which also works for exome sequencing data.

Regards,

Günter

Dear Gunter, thank you very much for your suggestions. If I may add, would the following R code be suitable for WES analysis on TUMOR-GERMLINE pairs ? I am including below the full R code that I have  put recently together :

library(cn.mops)
library(parallel)
library(snow)
library(RUnit)
library(BiocParallel)

# registered()
# register(MulticoreParam(workers = 1))
#### args <- commandArgs(TRUE)

TUMOR <- "zzz2.bam"        
GERMLINE <- "zzz1.bam"

#######################################
#######################################

segments <- read.table("hg38_refFlat.longest.isoform",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1], IRanges(segments[,2], segments[,3]))

###############

tumor <- getSegmentReadCountsFromBAM(TUMOR, GR=gr, parallel=1)
normal <- getSegmentReadCountsFromBAM(GERMLINE, GR=gr, parallel=1)

####################
####################

X <- tumor
values(X) <- cbind(values(tumor),values(normal))

#### NORMALIZE the DATA .. do we use POISSON normalization for exome sequencing ?

X.mode <- normalizeGenome(X, normType="poisson")

### Parameter settings for tumor:
### - norm=0, because we already have normalized
### - integer copy numbers higher than 8 allowed
### - DNAcopy as segmentation algorithm.

## the position 1 is for TUMOR, the position 2 is for NORMAL.

ref_analysis_norm0 <- referencecn.mops(X.mode[,1], X.mode[,2],
    norm=0,
    I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64),
    classes = c("CN0", "CN1", "CN2", "CN3", "CN4", "CN5", "CN6", "CN7","CN8","CN16","CN32","CN64","CN128"),
    segAlgorithm="DNAcopy")
 
resCNMOPS <- calcIntegerCopyNumbers(ref_analysis_norm0)
resCNMOPS <- cn.mops:::.replaceNames(resCNMOPS, basename(TUMOR) ,"TUMOR")

##############
name_file_display <- paste("z.cn.mops.seqplot", basename(TUMOR), basename(GERMLINE), basename(CHR), "png", sep=".")

png(name_file_display)
# segplot(resCNMOPS)
segplot(resCNMOPS, zlcol="blue", segcol="red")
dev.off()
 
###############

## cnvs(resCNMOPS) # look at individual CNV regions
## cnvr(resCNMOPS) # look at REGIONS

############### SAVING THE RESULTS :

results.segm  <- as.data.frame(segmentation(resCNMOPS))
results.CNVs  <- as.data.frame(cnvs(resCNMOPS))
results.CNVRegions  <- as.data.frame(cnvr(resCNMOPS))

############### NAMING THE FILES :

name_file_results_segm <- paste("z.cn.mops.segments.normalization.POISSON", basename(TUMOR), basename(GERMLINE), "txt", sep=".")
name_file_results.CNVs <- paste("z.cn.mops.cnvs_shorts.normalization.POISSON", basename(TUMOR), basename(GERMLINE), "txt", sep=".")
name_file_results.CNVRegions <- paste("z.cn.mops.cnvs_regions.normalization.POISSON", basename(TUMOR), basename(GERMLINE), "txt", sep=".")

############### WRITING THE TABLES :

write.table(results.segm, file=name_file_results_segm, sep="\t", row.names=F)
write.table(results.CNVs, file=name_file_results.CNVs, sep="\t", row.names=F)
write.table(results.CNVRegions, file=name_file_results.CNVRegions, sep="\t", row.names=F)