Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
# test if there is at least one argument: if not, return an error
if (length(args)==0) {
stop("Please provide prefix for output", call.=FALSE)
} else if (length(args)==1) {
## default bin size
args[2] = 500
}
### R code from vignette source 'QDNAseq.Rnw'
###################################################
### code chunk number 1: style-Sweave
###################################################
## BiocStyle::latex()
future::plan("multiprocess")
options(mc.cores = 8L)
###################################################
### code chunk number 2: QDNAseq.Rnw:23-24
###################################################
library(QDNAseq)
library(methods)
###################################################
### code chunk number 3: QDNAseq.Rnw:27-29
###################################################
options("QDNAseq::verbose"=NA)
options("QDNAseq::cache"=TRUE)
options(width=40)
nameprefix <- args[1]
cachedir <- "cache"
dir.create(cachedir, showWarnings=F)
R.cache::setCacheRootPath(path="./cache")
binsize <- as.numeric(args[2]) # It's one kb size
binsfile <- file.path(cachedir,paste0(nameprefix,"bins.", binsize))
bins <- NULL
if (!file_test("-f", binsfile)) {
bins <- getBinAnnotations(binSize = binsize)
saveRDS(bins, file= binsfile)
} else {
bins <- readRDS(binsfile)
}
###################################################
### code chunk number 4: QDNAseq.Rnw:66-78 (eval = FALSE)
###################################################
## readCounts <- binReadCounts(bins)
## # all files ending in .bam from the current working directory
##
## # or
##
## readCounts <- binReadCounts(bins, bamfiles='tumor.bam')
## # file 'tumor.bam' from the current working directory
##
## # or
##
readCounts <- binReadCounts(bins, path='data')
## # all files ending in .bam from the subdirectory 'tumors'
###################################################
### code chunk number 5: QDNAseq.Rnw:89-92
###################################################
## data(LGG150)
## readCounts <- LGG150
readCounts
###################################################
### code chunk number 6: QDNAseq.Rnw:98-99
###################################################
png(paste0(nameprefix,"rawprofile_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 7: rawprofile
###################################################
plot(readCounts, logTransform=FALSE, ylim=c(-50, 200))
highlightFilters(readCounts, logTransform=FALSE,
residual=TRUE, blacklist=TRUE)
###################################################
### code chunk number 8: QDNAseq.Rnw:106-107
###################################################
dev.off()
###################################################
### code chunk number 9: QDNAseq.Rnw:123-125
###################################################
readCountsFiltered <- applyFilters(readCounts,
residual=TRUE, blacklist=TRUE)
###################################################
### code chunk number 10: QDNAseq.Rnw:127-128
###################################################
png(paste0(nameprefix,"isobar_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 11: isobar
###################################################
isobarPlot(readCountsFiltered)
###################################################
### code chunk number 12: QDNAseq.Rnw:133-134
###################################################
dev.off()
###################################################
### code chunk number 13: QDNAseq.Rnw:152-153
###################################################
readCountsFiltered <- estimateCorrection(readCountsFiltered)
###################################################
### code chunk number 14: QDNAseq.Rnw:155-156
###################################################
png(paste0(nameprefix,"noise_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 15: noise
###################################################
noisePlot(readCountsFiltered)
###################################################
### code chunk number 16: QDNAseq.Rnw:161-162
###################################################
dev.off()
###################################################
### code chunk number 17: QDNAseq.Rnw:176-180
###################################################
## readCountsFilteredwSex <- applyFilters(readCountsFiltered, chromosomes=NA)
readCountsFilteredwSex <- applyFilters(readCountsFiltered)
copyNumbers <- correctBins(readCountsFilteredwSex)
copyNumbers
copyNumbersNormalized <- normalizeBins(copyNumbers)
copyNumbersSmooth <- smoothOutlierBins(copyNumbersNormalized)
###################################################
### code chunk number 18: QDNAseq.Rnw:182-183
###################################################
png(paste0(nameprefix,"profile_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 19: profile
###################################################
plot(copyNumbersSmooth)
###################################################
### code chunk number 20: QDNAseq.Rnw:188-189
###################################################
dev.off()
###################################################
### code chunk number 21: QDNAseq.Rnw:203-206 (eval = FALSE)
###################################################
## exportBins(copyNumbersSmooth, file="LGG150.txt")
## exportBins(copyNumbersSmooth, file="LGG150.igv", format="igv")
exportBins(copyNumbersSmooth, file=paste0(nameprefix,"%s_bin_", binsize, ".bed"), format="bed")
exportBins(copyNumbersSmooth, file=paste0(nameprefix,"smooth_bin_", binsize, ".igv"), format="igv")
###################################################
### code chunk number 22: QDNAseq.Rnw:219-221
###################################################
copyNumbersSegmented <- segmentBins(copyNumbersSmooth, transformFun="sqrt")
copyNumbersSegmented <- normalizeSegmentedBins(copyNumbersSegmented)
###################################################
### code chunk number 23: QDNAseq.Rnw:223-224
###################################################
png(paste0(nameprefix,"segment_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 24: segments
###################################################
plot(copyNumbersSegmented)
###################################################
### code chunk number 25: QDNAseq.Rnw:229-230
###################################################
dev.off()
###################################################
### code chunk number 26: QDNAseq.Rnw:243-244
###################################################
copyNumbersCalled <- callBins(copyNumbersSegmented, ncpus=6)
###################################################
### code chunk number 27: QDNAseq.Rnw:246-247
###################################################
png(paste0(nameprefix,"calls_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
###################################################
### code chunk number 28: calls
###################################################
plot(copyNumbersCalled)
###################################################
### code chunk number 29: QDNAseq.Rnw:252-253
###################################################
dev.off()
exportBins(copyNumbersCalled, file=paste0(nameprefix, "Calls_bin_", binsize, "_%s.bed"), format="bed")
exportBins(copyNumbersCalled, file=paste0(nameprefix, "Calls_bin_", binsize, ".igv"), format="igv")
###################################################
### code chunk number 30: QDNAseq.Rnw:278-280
###################################################
cgh <- makeCgh(copyNumbersCalled)
exportBins(cgh, file=paste0(nameprefix, "cgh_bin_", binsize, "_%s.bed"), format="bed")
png(paste0(nameprefix,"cgh_bin_",binsize, "_%03d.png"), width = 1920, height = 1080)
plot(cgh)
dev.off()
###################################################
### code chunk number 31: QDNAseq.Rnw:315-316 (eval = FALSE)
###################################################
## copyNumbers <- callBins(..., ncpus=4L)
###################################################
### code chunk number 32: QDNAseq.Rnw:324-325 (eval = FALSE)
###################################################
## future::plan("eager")
###################################################
### code chunk number 33: QDNAseq.Rnw:333-334 (eval = FALSE)
###################################################
## future::plan("multiprocess")
###################################################
### code chunk number 34: QDNAseq.Rnw:348-349 (eval = FALSE)
###################################################
## options(mc.cores=4L)
###################################################
### code chunk number 35: QDNAseq.Rnw:357-359 (eval = FALSE)
###################################################
## cl <- parallel::makeCluster(...)
## future::plan("cluster", cluster=cl)
###################################################
### code chunk number 36: QDNAseq.Rnw:379-384 (eval = FALSE)
###################################################
## readCounts <- binReadCounts(getBinAnnotations(15))
## readCounts <- applyFilters(readCounts)
## readCounts <- estimateCorrection(readCounts)
## readCounts <- applyFilters(readCounts, chromosomes=NA)
## copyNumbers <- correctBins(readCounts)
###################################################
### code chunk number 37: QDNAseq.Rnw:399-401 (eval = FALSE)
###################################################
## readCounts <- estimateCorrection(readCounts,
## control=loess.control(surface="direct"))
###################################################
### code chunk number 38: QDNAseq.Rnw:415-425 (eval = FALSE)
###################################################
## # load required packages for human reference genome build hg19
## library(QDNAseq)
## library(Biobase)
## library(BSgenome.Hsapiens.UCSC.hg19)
##
## # set the bin size
## binSize <- 15
##
## # create bins from the reference genome
## bins <- createBins(bsgenome=BSgenome.Hsapiens.UCSC.hg19, binSize=binSize)
###################################################
### code chunk number 39: QDNAseq.Rnw:442-446 (eval = FALSE)
###################################################
## # calculate mappabilites per bin from ENCODE mapability tracks
## bins$mappability <- calculateMappability(bins,
## bigWigFile='/path/to/wgEncodeCrgMapabilityAlign50mer.bigWig',
## bigWigAverageOverBed='/path/to/bigWigAverageOverBed')
###################################################
### code chunk number 40: QDNAseq.Rnw:458-462 (eval = FALSE)
###################################################
## # calculate overlap with ENCODE blacklisted regions
## bins$blacklist <- calculateBlacklist(bins,
## bedFiles=c('/path/to/wgEncodeDacMapabilityConsensusExcludable.bed',
## '/path/to/wgEncodeDukeMapabilityRegionsExcludable.bed'))
###################################################
### code chunk number 41: QDNAseq.Rnw:470-476 (eval = FALSE)
###################################################
## # load data for the 1000 Genomes (or similar) data set, and generate residuals
## ctrl <- binReadCounts(bins,
## path='/path/to/control-set/bam/files')
## ctrl <- applyFilters(ctrl, residual=FALSE, blacklist=FALSE,
## mappability=FALSE, bases=FALSE)
## bins$residual <- iterateResiduals(ctrl)
###################################################
### code chunk number 42: QDNAseq.Rnw:485-488 (eval = FALSE)
###################################################
## # by default, use all autosomal bins that have a reference sequence
## # (i.e. not only N's)
## bins$use <- bins$chromosome %in% as.character(1:22) & bins$bases > 0
###################################################
### code chunk number 43: QDNAseq.Rnw:494-507 (eval = FALSE)
###################################################
## # convert to AnnotatedDataFrame and add metadata
## bins <- AnnotatedDataFrame(bins,
## varMetadata=data.frame(labelDescription=c(
## 'Chromosome name',
## 'Base pair start position',
## 'Base pair end position',
## 'Percentage of non-N nucleotides (of full bin size)',
## 'Percentage of C and G nucleotides (of non-N nucleotides)',
## 'Average mappability of 50mers with a maximum of 2 mismatches',
## 'Percent overlap with ENCODE blacklisted regions',
## 'Median loess residual from 1000 Genomes (50mers)',
## 'Whether the bin should be used in subsequent analysis steps'),
## row.names=colnames(bins)))
###################################################
### code chunk number 44: QDNAseq.Rnw:513-521 (eval = FALSE)
###################################################
## attr(bins, "QDNAseq") <- list(
## author="Ilari Scheinin",
## date=Sys.time(),
## organism="Hsapiens",
## build="hg19",
## version=packageVersion("QDNAseq"),
## md5=digest::digest(bins@data),
## sessionInfo=sessionInfo())
###################################################
### code chunk number 45: QDNAseq.Rnw:531-564 (eval = FALSE)
###################################################
## # download table of samples
## server <- "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/"
## g1k <- read.table(paste0(nameprefix,server, "sequence.index"),
## header=TRUE, sep="\t", as.is=TRUE, fill=TRUE)
##
## # keep cases that are Illumina, low coverage, single-read, and not withdrawn
## g1k <- g1k[g1k$INSTRUMENT_PLATFORM == 'ILLUMINA', ]
## g1k <- g1k[g1k$ANALYSIS_GROUP == 'low coverage', ]
## g1k <- g1k[g1k$LIBRARY_LAYOUT == 'SINGLE', ]
## g1k <- g1k[g1k$WITHDRAWN == 0, ]
##
## # keep cases with read lengths of at least 50 bp
## g1k <- g1k[!g1k$BASE_COUNT %in% c("not available", ""), ]
## g1k$BASE_COUNT <- as.numeric(g1k$BASE_COUNT)
## g1k$READ_COUNT <- as.integer(g1k$READ_COUNT)
## g1k$readLength <- g1k$BASE_COUNT / g1k$READ_COUNT
## g1k <- g1k[g1k$readLength > 50, ]
##
## # keep samples with a minimum of one million reads
## readCountPerSample <- aggregate(g1k$READ_COUNT,
## by=list(sample=g1k$SAMPLE_NAME), sum)
## g1k <- g1k[g1k$SAMPLE_NAME %in%
## readCountPerSample$sample[readCountPerSample$x >= 1e6], ]
##
## g1k$fileName <- basename(g1k$FASTQ_FILE)
##
## # download fastq files
## for (i in rownames(g1k)) {
## sourceFile <- paste0(nameprefix,server, g1k[i, "FASTQ_FILE"])
## destFile <- g1k[i, "fileName"]
## if (!file.exists(destFile))
## download.file(sourceFile, destFile)
## }
###################################################
### code chunk number 46: QDNAseq.Rnw:579-580
###################################################
## toLatex(sessionInfo())
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.