cnmops_wes.R 2.32 KB
Newer Older
Wai Yi Leung's avatar
Wai Yi Leung committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28
#!/usr/bin/env Rscript
suppressPackageStartupMessages(library('cn.mops'))
suppressPackageStartupMessages(library('optparse'))

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

option_list <- list(
    make_option(c("--rawoutput"), dest="rawoutput"),
    make_option(c("--cnv"), dest="cnv"),
    make_option(c("--cnr"), dest="cnr"),
    make_option(c("--chr"), dest="chr"),
    make_option(c("--targetBed"), dest="targetBed"),
    make_option(c("--threads"), dest="threads", default=8, type="integer")
    )

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

BAMFiles <- c(bamFile)

29
### WES Specific code
Wai Yi Leung's avatar
Wai Yi Leung committed
30 31 32 33
segments <- read.table(opt$targetBed, sep="\t", as.is=TRUE)
# filter the segments by the requested chromosome
segments <- segments[ segments[,1] == chromosome, ]
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
34 35
### END WES Specific code

Wai Yi Leung's avatar
Wai Yi Leung committed
36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
bamDataRanges <- getSegmentReadCountsFromBAM(BAMFiles, GR=gr, mode="paired", parallel=opt$threads)

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

res <- exomecn.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)


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

dir.create(chromosome, showWarnings=FALSE, recursive=TRUE, mode="0744")

# Plot chromosome per sample.
for ( i in 1:length(BAMFiles)){
  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)
  segplot(res,sampleIdx=i)
  dev.off()
}

# Plot cnvr regions.
for ( i in 1:nrow(as.data.frame(cnvr(res)))) {
  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)
  plot(res,which=i,toFile=TRUE)
  dev.off()
}