Update with alterantive read.bam.tags function

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) {