##' Read FastQ input files
##'
##' Uses the global config to find input files
##' @return Reads as list of ShortRead objects
##' @export
##' @keywords internal
readInputFiles <- function() {
loginfo("io.R/readInputFiles: reading FASTQ file(s)...")
## get config parameters
paired_ends <- getConfig.logical("paired_ends")
filename1 <- getConfig("input_file")
if (paired_ends) filename2 <- getConfig("input_file2")
## check file existence
if (!file.exists(filename1)) stop("cannot open file indicated by config parameter input_file=", filename1)
if (paired_ends) if (!file.exists(filename2)) stop("cannot open file indicated by config parameter input_file2=", filename2)
## read files
lreads <- list(readFastq(filename1))
loginfo(paste("io.R/readInputFiles: read", length(lreads[[1]]), "reads from input_file=", filename1))
if (paired_ends) {
lreads <- c(lreads, list(readFastq(filename2)))
loginfo(paste("io.R/readInputFiles: read", length(lreads[[2]]), "reads from input_file2=", filename2))
}
loginfo("io.R/readInputFiles: done")
return(lreads)
}
##' Write reads to file
##'
##' @param lreads List of reads as ShortRead objects
##' @param dir Save directory
##' @param filename1 Name of file 1
##' @param filename2 Name of file 2
##' @return Named list of filepaths
##' @export
##' @keywords internal
##' @importMethodsFrom ShortRead writeFastq
writeFastQFiles <- function (lreads, dir, filename1, filename2) {
## get config parameters
paired_ends <- getConfig.logical("paired_ends")
## write forward reads
filepaths <- list()
filepaths$fastq_for_aligner_1 <- file.path(dir, filename1)
loginfo(paste("io.R/writeFastQFiles: writing filename=", filepaths$fastq_for_aligner_1))
if (length(lreads[[1]])>0) writeFastq(lreads[[1]], file=filepaths$fastq_for_aligner_1, lane="", compress=FALSE)
else cat("", file=filepaths$fastq_for_aligner_1) ## empty FASTQ file
## write reverse reads
if (paired_ends) {
filepaths$fastq_for_aligner_2 <- file.path(dir, filename2)
loginfo(paste("io.R/writeFastQFiles: writing filename=", filepaths$fastq_for_aligner_2))
if (length(lreads[[2]])>0) writeFastq(lreads[[2]], file=filepaths$fastq_for_aligner_2, lane="", compress=FALSE)
else cat("", file=filepaths$fastq_for_aligner_2) ## empty FASTQ file
}
return(filepaths)
}
##' Count reads in Fastq file
##'
##' @title Count reads in Fastq file
##' @param filename Name of FastQ file
##' @return Number of reads
##' @author Gregoire Pau
##' @export
##' @importMethodsFrom ShortRead FastqStreamer
##' @keywords internal
getNumberOfReadsInFASTQFile <- function(filename) {
if (length(grep("\\.gz$", filename))>0) con <- gzfile(filename)
else con <- file(filename)
fqs <- FastqStreamer(con, n=1e6)
n <- 0
repeat {
nreads <- length(safe.yield(fqs))
if (nreads==0) break
else n <- n + nreads
}
close(fqs)
n
}
##' Open a streaming connection to a FastQ file
##'
##' Only one FastQStreamer object can be open at any time.
##'
##' @title Open a streaming connection to a FastQ file
##' @param input_file Path to a FastQ file
##' @param input_file2 Optional path to a FastQ file. Default is NULL.
##' @param chunk_size Number of reads per chunk
##' @param subsample_nbreads Optional number of reads to subsample (deterministic) from the input files. Default is NULL.
##' @param max_nbchunks Optional maximal number of chunks to read
##' @seealso FastQStreamer.getReads
##' @return Nothing.
##' @author Gregoire Pau
##' @importFrom stats runif
##' @export
##' @keywords internal
FastQStreamer.init <- function(input_file, input_file2=NULL, chunk_size, subsample_nbreads=NULL, max_nbchunks=NULL) {
## check file existences
if (!file.exists(input_file)) stop("cannot open file indicated by parameter input_file=", input_file)
if (!is.null(input_file2)) {
if (!file.exists(input_file2)) stop("cannot open file indicated by parameter input_file2=", input_file2)
}
## build subsampling_filter if subsample_nbreads is not NULL
subsampling_filter <- NULL
if (!is.null(subsample_nbreads)) {
loginfo(paste("io.R/FastQStreamer.init: counting number of reads in file=", input_file))
totnbreads <- getNumberOfReadsInFASTQFile(input_file)
if (subsample_nbreads<totnbreads) {
## save and set fixed random seed, to allow reproducible behaviors
runif(1) ; originalseed <- .Random.seed ; set.seed(1)
## build subsampling filter
subsampling_filter <- rep(FALSE, totnbreads)
subsampling_filter[sample(seq_len(totnbreads), subsample_nbreads)] <- TRUE
## restore seed
set.seed(originalseed)
loginfo(paste("io.R/FastQStreamer.init: subsampling_filter set (subsampled reads=", subsample_nbreads,
", totnbreads=", totnbreads, ")"))
} else {
loginfo(paste("io.R/FastQStreamer.init: the requested number of subsampled reads=", subsample_nbreads,
"is higher than the total number of reads=", totnbreads, ": I won't subsample"))
}
}
## initialise streamers
filenames <- list(input_file)
if (!is.null(input_file2)) filenames <- c(filenames, input_file2)
lfqs <- lapply(filenames, function (filename) {
if (length(grep("\\.gz$", filename))>0) con1 <- gzfile(filename)
else con1 <- file(filename)
fqs <- FastqStreamer(con1, n=chunk_size)
loginfo(paste("io.R/FastQStreamer.init: initialised FastQ streamer for filename=", filename))
fqs
})
FastQStreamer.lfqs <<- lfqs
FastQStreamer.chunkid <<- 0
FastQStreamer.subsampling_filter <<- subsampling_filter
if (is.null(max_nbchunks)) FastQStreamer.max_nbchunks <<- Inf
else FastQStreamer.max_nbchunks <<- as.integer(max_nbchunks)
}
##' Get FastQ reads from the FastQ streamer
##'
##' @title Get FastQ reads from the FastQ streamer
##' @return A list of ShortRead object containing reads. NULL if there are no more reads to read.
##' @author Gregoire Pau
##' @seealso FastQStreamer.init
##' @export
##' @keywords internal
FastQStreamer.getReads <- function() {
## check
if (!exists("FastQStreamer.lfqs")) {
FastQStreamer.lfqs <- NULL ## to avoid the 'no visible binding' R CMD check warning
FastQStreamer.chunkid <- NULL ## to avoid the 'no visible binding' R CMD check warning
FastQStreamer.max_nbchunks <- NULL ## to avoid the 'no visible binding' R CMD check warning
FastQStreamer.subsampling_filter <- NULL ## to avoid the 'no visible binding' R CMD check warning
stop("io.R/FastQStreamer.getReads: FastQStreamer.init() has not been called")
}
## get reads
lfqs <- FastQStreamer.lfqs
paired_ends <- length(lfqs)==2
lreads <- list(safe.yield(lfqs[[1]]))
if (paired_ends) {
lreads <- c(lreads, list(safe.yield(lfqs[[2]])))
if (length(lreads[[1]])!=length(lreads[[2]]))
stop("io.R/FastQStreamer.getReads: input files must have the same number of reads")
}
## increment chunkid
FastQStreamer.chunkid <<- FastQStreamer.chunkid + 1
if (length(lreads[[1]])==0 || FastQStreamer.chunkid>FastQStreamer.max_nbchunks) lreads <- NULL
else {
## subsample, if required
if (!is.null(FastQStreamer.subsampling_filter)) {
len <- length(lreads[[1]])
if (len>0) {
z <- FastQStreamer.subsampling_filter[1:len]
## if (sum(z)>0) {
lreads <- lapply(lreads, function(reads) reads[z])
FastQStreamer.subsampling_filter <<- FastQStreamer.subsampling_filter[-(1:len)]
## } else {
## stop("io.R/FastQStreamer.getReads: (known bug) the chunk is empty after subsampling. I won't process an empty chunk! Please disable or increase 'subsample_nbreads', or 'increase chunk_size'")
##}
}
}
}
lreads
}
##' Close the FastQStreamer
##'
##' @title Close the FastQStreamer
##' @return Nothing
##' @seealso FastQStreamer.init
##' @author Gregoire Pau
##' @export
##' @keywords internal
FastQStreamer.release <- function() {
## check
if (!exists("FastQStreamer.lfqs")) {
FastQStreamer.lfqs <- NULL ## to avoid the 'no visible binding' R CMD check warning
stop("io.R/FastQStreamer.release: FastQStreamer.init() has not been called")
}
invisible(lapply(FastQStreamer.lfqs, close))
}
##' Save an R object
##'
##' Exists so objects can be serialized and reloaded with the a unique
##' identifier in the symbol. Stores the data object with a new nane
## 'created by appending 'orig_name' and 'id'
##' @param data The data to store
##' @param orig_name The original name of the data
##' @param id A meaningful id the is prepended to the stored objects name
##' @param save_dir The directory where the data should be saved in
##' @param compress Save the data compressed or not
##' @param format Choice of 'RData' or 'tab'(ular)
##' @return Name of the stored file
##' @importFrom utils write.table
##' @export
##' @keywords internal
saveWithID <- function(data, orig_name, id, save_dir, compress=TRUE, format="RData") {
new_name <- paste(id, orig_name, sep=".")
if (format=="RData") {
## RData format
assign(new_name, data)
filename <- paste(new_name, ".RData", sep="")
filename <- file.path(save_dir, filename)
loginfo(paste("io.R/saveWithID: saving file=", filename))
save(list=new_name, file=filename, compress=compress)
} else if (format=="tab"){
## tab-separated format
filename <- paste(new_name, ".tab", sep="")
filename <- file.path(save_dir, filename)
loginfo(paste("io.R/saveWithID: saving file=", filename))
write.table(data, filename, row.names=FALSE, sep="\t", quote=FALSE)
}
else{
stop("saveWithID: unknown format: ", format)
}
return(filename)
}
##' Make a directory after performing an existence check
##'
##' Throws an exception if file or directory with same name exist
##' and overwrite is TRUE.
##' @param dir Name of directory to create
##' @param overwrite A character string: never (default), erase, overwrite
##' @return Path to created directory
##' @export
##' @keywords internal
makeDir <- function(dir, overwrite="never"){
if (file.exists(dir)) {
if (overwrite=="never") stop("io.R/makeDir: I won't overwrite data present in dir=", dir)
if (overwrite=="erase") safeUnlink(dir)
}
dir.create(dir, recursive=TRUE, showWarnings=FALSE)
return(dir)
}
##' Create a random directory with prefix in R temp dir
##'
##' Especially for testing code it is very helpful to have
##' a temp directory with a defined prefix, so one knows which test
##' produced which directory.
##' @param prefix A string that will preceed the directory name
##' @param dir Directory where the random dir will be created under.
##' Defaults to tempdir()
##' @return Name of temporary directory
##' @keywords internal
createTmpDir <- function(prefix=NULL, dir=tempdir()) {
tmpname <- NULL
while(is.null(tmpname)){
tmpname <- try(
{
tmp <- file.path(dir, paste(c(prefix, sample(letters,8)), collapse=""))
makeDir(tmp, overwrite="never")
return(tmp)
}, silent=TRUE)
if (class(tmpname)=="try-error") {
tmpname <- NULL ## restarts while
}
}
return(tmpname)
}
##' Creates a directory with all needed subdirectories for pipeline
##' outputs
##'
##' @title Create output directory and subdirectories for sequencing pipeline analysis outputs
##' @param save_dir path to the directory that will contain all
##' needed subdirectories
##' @param overwrite A character string: never (default), erase, overwrite
##' @return Nothing. Called for its side effects
##' @author Cory Barr, Jens Reeder
##' @export
##' @keywords internal
setUpDirs <- function(save_dir, overwrite="never") {
makeDir(save_dir, overwrite=overwrite)
dirs <- c("bams", "logs", "results", "reports", "reports/images")
dirs <- file.path(save_dir, dirs)
invisible(sapply(dirs, makeDir, overwrite=overwrite))
}
##' Safely load a R data file
##'
##' Attempts to load a file given by object_name.
##' Bails out if none or more than one files match the object name.
##'
##' @param dir_path Save dir of a pipeline run
##' @param object_name object name, can be a regexp
##' @return loaded object
##' @export
##' @keywords internal
safeGetObject <- function(dir_path, object_name){
filename <- getObjectFilename(file.path(dir_path,"results"), object_name)
obj <- try({
get(load(filename))
}, silent=TRUE)
if (class(obj)=="try-error") {
stop("error: cannot load ", object_name, "in", dir_path, "\n", sep="")
}
return(obj)
}
##' Load tabular data from the NGS pipeline result directory
##' @param save_dir A character string containing an NGS pipeline output directory.
##' @param object_name A character string ontaining the regular expression matching a filename in dir_path
##' @return A data frame.
##' @export
getTabDataFromFile <- function(save_dir, object_name) {
filename <- getObjectFilename(file.path(save_dir, "results"), paste(".", object_name, "\\.tab$", sep=""))
read.table(filename, sep="\t", header=TRUE, stringsAsFactors=FALSE)
}
##' Load data as numerical values
##'
##' @param dir_path Save dir of a pipeline run
##' @param object_name Object name
##' @return loaded data as table of numbers
##' @author Jens Reeder
##' @export
##' @keywords internal
getNumericVectorDataFromFile <- function(dir_path, object_name) {
df <- getTabDataFromFile(dir_path, object_name)
v <- as.numeric(df$value)
names(v) <- df$name
return(v)
}
##' Detect quality protocol from a FASTQ file
##'
##' @param filename Path to a FASTQ or gzipped-FASTQ file
##' @param nreads Number of reads to test quality on. Default is 5000.
##' @return A character vector containing the compatible qualities. NULL if none.
##' @author Jens Reeder
##' @export
##' @importMethodsFrom IRanges as.vector
##' @keywords internal
detectQualityInFASTQFile <- function(filename, nreads=5000) {
cquals <- try({
if (length(grep("\\.gz$", filename))>0) con <- gzfile(filename)
else con <- file(filename)
## read qualities
fqs <- FastqStreamer(con, n=nreads)
reads <- safe.yield(fqs)
z <- paste(as.vector(quality(quality(reads))), collapse="")
quals <- as.numeric(charToRaw(z))
close(con)
## check
if (length(quals)==0) stop()
## build a dictionary of known qualities
## "illumina1.5" upper bound is now 105 instead of 104, to accomodate Phred-qualities of 41 (instead of 40)
## some bam files from the TCGA project come with rescaled phred qualities from 1-50
## to accomodate for them we invent a new quality range "GATK-rescaled"
knownquals <- data.frame(qual=c("solexa", "illumina1.3", "illumina1.5", "illumina1.8", "GATK-rescaled", "sanger"),
min=c(59, 64, 66, 33, 33, 33),
max=c(104, 104, 105, 74, 83, 73), stringsAsFactors=FALSE)
## check compatible qualities
rquals <- range(quals)
unlist(lapply(1:nrow(knownquals), function(i) {
if (rquals[1]>=knownquals$min[i] && rquals[2]<=knownquals$max[i]) knownquals$qual[i]
else NULL
}))
}, silent=TRUE)
if (class(cquals)=="try-error") stop("io.R/detectQualityInFASTQFile: cannot read FASTQ reads in filename=", filename, cquals)
else cquals
}
##' Symlink-safe file/directory delete function
##'
##' Unlike unlink(), safeUnlink() does not follow symlink directories for deletion.
##'
##' @title safeUnlink
##' @param path A character string indicating which file/directory to delete.
##' @return Nothing
##' @author Gregoire Pau
##' @export
##' @keywords internal
safeUnlink <- function(path) {
cmd <- paste("rm -rf", path, "2>&1")
z <- suppressWarnings(system(cmd, intern=TRUE))
if (length(z)>0) stop("io.R/safeUnlink: cannot delete file=", path, "; ", z)
}
##' Get a filename given a directory and the object name
##'
##' @title Get a filename given a directory and the object name
##' @param dir_path A character string containing a dir path
##' @param object_name A character string containing the regular expression matching a filename in dir_path
##' @return A character vector containing an existing filename, stops if 0 or more than 1
##' @author Gregoire Pau
##' @export
##' @keywords internal
getObjectFilename <- function(dir_path, object_name) {
files <- dir(dir_path, full.names=TRUE)
filename <- grep(paste(object_name, sep=""), files, value=TRUE)
if (length(filename)!=1) {
stop("io.R/getObjectFilename: found non-unique (0 or multiple) filename matching '", object_name,
"' in directory ", dir_path)
}
return(filename)
}
##' Get the filename of the variant file
##'
##' Depending on the variant caller used and the version of
##' VariantAnotation used to create the file a file might have the
##' ending vcf.gz, vcf.bgz. To function hides all this mess.
##' @title Get a vcf filename given a HTSeqGenie directory
##' @param dir_path A character string containing a dir path
##' @return A character vector containing an existing filename, stops if 0 or more than 1
##' @author Jens Reeder
##' @export
##' @keywords internal
findVariantFile <- function(save_dir){
vcf <- getObjectFilename(file.path(save_dir, "results"),
"variants.vcf.b?gz$")
return(vcf)
}
## this path is ciontructed over and over again, so we factor it out here
## similar functions for other files will follow
getAnalyzedBamFile <- function(){
file <- file.path(getConfig('save_dir'), 'bams',
paste(getConfig('prepend_str'),
'analyzed.bam', sep="."))
return(file)
}
##' Overloaded yield(...) method catching truncated exceptions for FastqStreamer
##'
##' @title Overloaded yield(...) method catching truncated exceptions for FastqStreamer
##' @param fqs An instance from the FastqSampler or FastqStreamer class.
##' @return Same as FastqStreamer::yield
##' @author Gregoire Pau
##' @export
##' @keywords internal
##' @importMethodsFrom ShortRead yield
safe.yield <- function(fqs) {
tryCatch(yield(fqs), IncompleteFinalRecord=function(x) stop("io.R/safe.yield: input gzipped file is corrupted/truncated. Aborting.\n"))
}
## file.rename only warns, but we like it to fail instead
safe.file.rename <- function(from, to){
success <- file.rename(from, to)
if(!all(success)){
stop(paste("Renaming file", from, "to", to, "failed"))
}
}
##' Read single/paired end BAM files with requested columns from the BAM
##'
##' @title Read single/paired End Bam Files
##' @param filename Path to a bam file
##' @param paired_ends A logical indicating whether the reads are paired
##' @param remove.strandness A logical indicating whether read strands should be set to "*".
##' @return GRangesList
##' @author Cory Barr
##' @export
##' @keywords internal
##' @importMethodsFrom Rsamtools ScanBamParam
##' @importFrom Rsamtools scanBamFlag
##' @importMethodsFrom GenomicAlignments readGAlignments grglist
##' @importFrom Rsamtools BamFile
readRNASeqEnds <- function(filename, paired_ends, remove.strandness=TRUE) {
## bam file
scf <- scanBamFlag(isNotPassingQualityControls=FALSE, isDuplicate=FALSE)
sbp <- ScanBamParam(flag=scf, what="groupid")
if (paired_ends) bamfile <- BamFile(filename, asMates=TRUE)
else bamfile <- BamFile(filename)
## read bam
reads <- readGAlignments(bamfile, param=sbp, use.names=FALSE, with.which_label=FALSE)
if (remove.strandness) strand(reads) <- "*"
groupid <- values(reads)$groupid
reads <- grglist(reads) ## use cigar information to split gapped reads into granges
## pair reads
if (paired_ends) {
z <- rep(groupid, elementNROWS(reads))
reads <- unlist(reads)
reads <- split(reads, z)
}
## return reads
return(reads)
}
