From f9441d377be5f385a5235856f4b1619bc60565d5 Mon Sep 17 00:00:00 2001 From: fferrari Date: Tue, 2 Mar 2021 15:46:21 +0100 Subject: [PATCH] Update with alterantive read.bam.tags function The alternative function will randomly pick only one of the reads in the pair when both are aligned, thus resulting in a taglist object with characteristics similar to a single end BAM file. however, I am not sure this was the intended purpose of the original code handling BAM files with paired end reads. --- R/zroutines.R | 66 +++++++++++++++++++++++++++++++++++++++++---------- 1 file changed, 53 insertions(+), 13 deletions(-) diff --git a/R/zroutines.R b/R/zroutines.R index 6e46b37..7e0b036 100644 --- a/R/zroutines.R +++ b/R/zroutines.R @@ -112,31 +112,70 @@ read.bowtie.tags <- function(filename,read.tag.names=F,fix.chromosome.names=F) { } } -read.bam.tags <- function(filename,read.tag.names=F,fix.chromosome.names=F) { +## this function reads BAM files into a taglist object of SPP +## this function code is derived from the original "read.bam.tags" from spp package by Peter Kharchenko +## the version included here has been revised to hanlde BAM files containing paired end reads +## the argument read.tag.names is here should be set to TRUE to properly handle paired end reads in BAM files +f_read.bam.tags <- function(filename,read.tag.names=F,fix.chromosome.names=F) { #require(Rsamtools) if(!is.element("Rsamtools", installed.packages()[, 1])) { stop("Rsamtools Bioconductor package is now required for BAM file support. Please install") } - + +## this is setting the list of fileds to be extracted from the BAM file +## note that "pos" is the mapping position on the reference sequence (chromosome name stored in "rname") +## whereas the "qname" field contains the ID of each sequencing reads (or read pairs for paired end reads in the bam file) ww <- c("flag","rname","pos","isize","strand","mapq","qwidth"); if(read.tag.names) { ww <- c(ww,"qname") }; bam <- Rsamtools::scanBam(filename,param=Rsamtools::ScanBamParam(what=ww,flag=Rsamtools::scanBamFlag(isUnmappedQuery=FALSE)))[[1]]; + +## this is returning an empty tagglist object if the BAM file contains no valid alignment if(is.null(bam$pos) || length(bam$pos)==0) { return(list(tags=c(),quality=c())) } + +## this is just creating a 1/0 vector for postiive/negative strand mapped reads strm <- as.integer(bam$strand=="+") - if(any(bitwAnd(bam$flag,0x1))) { # paired-end data - # use only positive strand mappings - posvi <- which(strm==1); - rl <- list(tags=tapply(posvi,bam$rname[posvi],function(ii) as.numeric(bam$pos[ii])), - flen=tapply(posvi,bam$rname[posvi],function(ii) as.numeric(abs(bam$isize[ii])))) - # alternatively, handle reads with NA isize (unpaired?) just like single-ended reads - #pos <- tapply(1:length(bam$pos),bam$rname,function(ii) ifelse(is.na(bam$isize[ii]), bam$pos[ii]*strm[ii] - (1-strm[ii])*(bam$pos[ii]+bam$qwidth[ii]), strm[ii]*bam$pos[ii] - (1-strm[ii])*(bam$pos[ii]+bam$isize[ii]))) + +## this is checking if the BAM file is actually containing paired end reads + if(any(bitwAnd(bam$flag,0x1))) { + if (!read.tag.names) { + stop("read.tag.names must be set to TRUE to handle paired end reads BAM files") + } + # paired-end data + ## for paired end data, we can select one (random) read out of the pair, so as to have equally represented both the positive and the negative strand mapped reads + ## we must design the code so as to take 1 random read for each read ID (qname) so that we get 1 even if we have only one read mapped in the pair + oneSelectedInPair<-unlist(tapply(X=1:length(bam$pos), INDEX=bam$qname, FUN=function(ii) { + if (length(ii)>1) { + return(sample(ii, size=1)) + } else if (length(ii)==1) { + return(ii) + } else { + stop("unexpected BAM file content format") + } + })) # return only one for each pair (or one for each group if more than 2 alignments are present) + ## the selection of indexes (oneSelectedInPair) is performed on the full vectors, thus we can use these indexes to perfrm subselections on the full BAM vectors + rl <- list(tags=tapply(X=oneSelectedInPair, INDEX=bam$rname[oneSelectedInPair],function(ii) bam$pos[ii]*strm[ii] - (1-strm[ii])*(bam$pos[ii]+bam$qwidth[ii]))) + rl <- c(rl,list(quality=tapply(X=oneSelectedInPair,INDEX=bam$rname[oneSelectedInPair],function(ii) bam$mapq[ii]))) + + ## return also the read IDS (query name = "qname") if required + if(read.tag.names) { + rl <- c(rl,list(names=tapply(X=oneSelectedInPair,INDEX=bam$rname[oneSelectedInPair],function(ii) bam$qname[ii]))) + } + } else { + ## this is the "standard" workflow in case the BAM file contains only single end reads + ## most of the ChIP-seq peaks calling alogorithms expect single end reads, disributed on both positive and negative strand + ## this line of code is traversing the BAM file content (1:length(bam$pos)), chromosome by chromosome (bam$rname) + ## and keeping the annotated position as 5'end of positive strand mapped reads (bam$pos[ii]*strm[ii] ) + ## or the "-" 3'end (i.e. the 5'-end of the reads mapped on hte negative strand (- (1-strm[ii])*(bam$pos[ii]+bam$qwidth[ii]))) + ## this "ifelse" condition for positive and negative strand reads is actually managed by the 1/0 vectors for strand, + ## as it will change to "zero" either the first or the second element in the subtraction below rl <- list(tags=tapply(1:length(bam$pos),bam$rname,function(ii) bam$pos[ii]*strm[ii] - (1-strm[ii])*(bam$pos[ii]+bam$qwidth[ii]))) + rl <- c(rl,list(quality=tapply(1:length(bam$pos),bam$rname,function(ii) bam$mapq[ii]))) + ## return also the read IDS (query name = "qname") if required + if(read.tag.names) { + rl <- c(rl,list(names=tapply(1:length(bam$pos),bam$rname,function(ii) bam$qname[ii]))) + } } - rl <- c(rl,list(quality=tapply(1:length(bam$pos),bam$rname,function(ii) bam$mapq[ii]))) - if(read.tag.names) { - rl <- c(rl,list(names=tapply(1:length(bam$pos),bam$rname,function(ii) bam$qname[ii]))) - } if(fix.chromosome.names) { # remove ".fa" names(rl) <- gsub("\\.fa","",names(rl)) @@ -145,6 +184,7 @@ read.bam.tags <- function(filename,read.tag.names=F,fix.chromosome.names=F) { } + read.helicos.tags <- function(filename,read.tag.names=F,fix.chromosome.names=F,include.length.info=T) { if(read.tag.names) { rtn <- as.integer(1); } else { rtn <- as.integer(0); }; tl <- lapply(.Call("read_helicostabf",path.expand(filename),rtn),function(d) {