Question: using cn.MOPS on WES data GERMLINE-TUMOR
0
gravatar for Bogdan
23 months ago by
Bogdan550
Palo Alto, CA, USA
Bogdan550 wrote:

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 • 444 views
ADD COMMENTlink modified 23 months ago • written 23 months ago by Bogdan550
Answer: using cn.MOPS on WES data GERMLINE-TUMOR
3
gravatar for Günter Klambauer
23 months ago by
Austria
Günter Klambauer540 wrote:

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 COMMENTlink written 23 months ago by Günter Klambauer540
Answer: using cn.MOPS on WES data GERMLINE-TUMOR
1
gravatar for Bogdan
23 months ago by
Bogdan550
Palo Alto, CA, USA
Bogdan550 wrote:

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 COMMENTlink written 23 months ago by Bogdan550
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 340 users visited in the last hour