cnmops.R 2.01 KB
Newer Older
Wai Yi Leung's avatar
Wai Yi Leung committed
1 2 3
#!/usr/bin/env Rscript
suppressPackageStartupMessages(library('cn.mops'))
suppressPackageStartupMessages(library('optparse'))
Wai Yi Leung's avatar
Wai Yi Leung committed
4 5 6 7

# Script from  https://git.lumc.nl/lgtc-bioinformatics/gapss3/blob/master/src/CNV/makeCnmops.sh
# modified to take arguments

Wai Yi Leung's avatar
Wai Yi Leung committed
8 9 10 11
option_list <- list(
    make_option(c("--rawoutput"), dest="rawoutput"),
    make_option(c("--cnv"), dest="cnv"),
    make_option(c("--cnr"), dest="cnr"),
Wai Yi Leung's avatar
Wai Yi Leung committed
12 13
    make_option(c("--chr"), dest="chr"),
    make_option(c("--threads"), dest="threads", default=8, type="integer")
Wai Yi Leung's avatar
Wai Yi Leung committed
14
    )
Wai Yi Leung's avatar
Wai Yi Leung committed
15

Wai Yi Leung's avatar
Wai Yi Leung committed
16 17 18 19 20 21 22 23 24
parser <- OptionParser(usage = "%prog [options] file", option_list=option_list)
arguments = parse_args(parser, positional_arguments=TRUE)
opt = arguments$options
args = arguments$args

chromosome <- opt$chr
CNVoutput <- opt$cnv
CNRoutput <- opt$cnr
bamFile <- args
Wai Yi Leung's avatar
Wai Yi Leung committed
25 26

BAMFiles <- c(bamFile)
Wai Yi Leung's avatar
Wai Yi Leung committed
27
bamDataRanges <- getReadCountsFromBAM(BAMFiles, mode="paired", refSeqName=chromosome, WL=1000, parallel=opt$threads)
Wai Yi Leung's avatar
Wai Yi Leung committed
28 29 30

write.table(as.data.frame( bamDataRanges ), quote = FALSE, opt$rawoutput, row.names=FALSE)

Wai Yi Leung's avatar
Wai Yi Leung committed
31 32 33 34 35 36
res <- cn.mops(bamDataRanges)
res <- calcIntegerCopyNumbers(res)

write.table(as.data.frame(cnvs(res)), quote = FALSE, CNVoutput, row.names=FALSE)
write.table(as.data.frame(cnvr(res)), quote = FALSE, CNRoutput, row.names=FALSE)

Wai Yi Leung's avatar
Wai Yi Leung committed
37 38 39 40 41

ppi <- 300
plot_margins <- c(3,4,1,2)+0.1
label_positions <- c(2,0.5,0)

Wai Yi Leung's avatar
Wai Yi Leung committed
42
dir.create(chromosome, showWarnings=FALSE, recursive=TRUE, mode="0744")
Wai Yi Leung's avatar
Wai Yi Leung committed
43

Wai Yi Leung's avatar
Wai Yi Leung committed
44 45
# Plot chromosome per sample.
for ( i in 1:length(BAMFiles)){
Wai Yi Leung's avatar
Wai Yi Leung committed
46 47 48 49 50 51 52
  png(file=paste(chromosome,"/",chromosome,"-segplot-",i,".png", sep=""),
  width = 16 * ppi, height = 10 * ppi,
    res=ppi, bg = "white"
  )
    par(mfrow = c(1,1))
    par(mar=plot_margins)
    par(mgp=label_positions)
Wai Yi Leung's avatar
Wai Yi Leung committed
53 54 55 56 57 58
  segplot(res,sampleIdx=i)
  dev.off()
}

# Plot cnvr regions.
for ( i in 1:nrow(as.data.frame(cnvr(res)))) {
Wai Yi Leung's avatar
Wai Yi Leung committed
59 60 61 62 63 64
  png(file=paste(chromosome,"/",chromosome,"-cnv-",i,".png",sep=""),
  width = 16 * ppi, height = 10 * ppi,
    res=ppi, bg = "white")
    par(mfrow = c(1,1))
    par(mar=plot_margins)
    par(mgp=label_positions)
Wai Yi Leung's avatar
Wai Yi Leung committed
65 66 67 68
  plot(res,which=i,toFile=TRUE)
  dev.off()
}