Question: using cn.MOPS on WES data GERMLINE-TUMOR
0
gravatar for Bogdan
2.5 years ago by
Bogdan580
Palo Alto, CA, USA
Bogdan580 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 • 538 views
ADD COMMENTlink modified 2.5 years ago • written 2.5 years ago by Bogdan580
Answer: using cn.MOPS on WES data GERMLINE-TUMOR
3
gravatar for Günter Klambauer
2.5 years 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 2.5 years ago by Günter Klambauer540
Answer: using cn.MOPS on WES data GERMLINE-TUMOR
1
gravatar for Bogdan
2.5 years ago by
Bogdan580
Palo Alto, CA, USA
Bogdan580 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 2.5 years ago by Bogdan580
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: 272 users visited in the last hour