using cn.MOPS on WES data GERMLINE-TUMOR
2
0
Entering edit mode
Bogdan ▴ 620
@bogdan-2367
Last seen 21 days ago
Palo Alto, CA, USA

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

cn.mops • 757 views
ADD COMMENT
3
Entering edit mode
@gunter-klambauer-5426
Last seen 8 months ago
Austria

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

ADD COMMENT
1
Entering edit mode
Bogdan ▴ 620
@bogdan-2367
Last seen 21 days ago
Palo Alto, CA, USA

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)
ADD COMMENT

Login before adding your answer.

Traffic: 240 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6