Hi Guys,

I am using cn.mops for CNV calling for whole exome data with one normal and three tumor samples. I am running individually for each tumor sample (code below) and getting an output csv file containing seq name, start, end, width, strand, tumor, gene-ids (output.jpeg; attached). I am not sure what CN2, CN4, CN32 values are ?. Are these values corresponding to fold change copy number as mentioned in code [I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64), classes = paste0("CN", c(0:8, 16, 32, 64, 128)). So, if there is CN4: that means copy number fold change is 2 and CN128 corresponds to fold change 64. Am I understanding correctly?

BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene") library("TxDb.Hsapiens.UCSC.hg19.knownGene")
suppressPackageStartupMessages({ library(cn.mops) })

1. Load the input data

bamFiles=c("200.bam","836.bam")

2. We can bin and count the reads

reads_gr <- getReadCountsFromBAM(BAMFiles = bamFiles, sampleNames = c("tumor", "normal"),WL = 10000, parallel = 20)
write.csv(reads_gr,'reads-200.csv')
X <- normalizeGenome(reads_gr)

4. Detect cnv's

ref_analysis <- referencecn.mops(X[,1], X[,2], norm=0, I = c(0.025, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 8, 16, 32, 64), classes = paste0("CN", c(0:8, 16, 32, 64, 128)), segAlgorithm="DNAcopy")
resCNMOPS <-calcIntegerCopyNumbers(ref_analysis)

5. Plot the results

pdf("Cnmopsplots-200.pdf", height = 9, width = 18)
segplot(resCNMOPS).
dev.off()

6. add gene information to your segments

human_cn <- cnvr(resCNMOPS)
txdb <-keepSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, c("chr1","chr2","chr3","chr4","chr5","chr6","chr7","chr8","chr9","chr10","chr11","chr12","chr13","chr14","chr15","chr16","chr17","chr18","chr19","chr20","chr21","chr22","chrX","chrY")).
head(seqlevels(txdb))#chromosomes that are active
genes0 <-genes(txdb).

'unlist' so that each range is associated with a single gene identifier

idx <- rep(seq_along(genes0), elementNROWS(genes0$gene_id)).
genes <- granges(genes0)[idx]
genes$gene_id = unlist(genes0$gene_id).
olaps <- findOverlaps(genes, human_cn, type="within")
idx <-factor(subjectHits(olaps), levels=seq_len(subjectLength(olaps))).
human_cn$gene_ids <- splitAsList(genes$gene_id[queryHits(olaps)], idx).
human_cn

7. Save output in a csv file

write.csv(human_cn,'tumor_cn-200.csv')

I will really appreciate it if you respond to my query.

Thanks