Hi to the experts,

First of all, I would like to thank you for the efforts you put in to the Bioconductor team. Those efforts make the lives of struggling bioinformatics students like me much easy.

This is a very naive inquiry about an issue I face with counting features using summarizeOverlaps() function.

I have an analysis on stranded, paired end RNA-Seq data of transgenic Arabidopsis. A snippet of my code used is below.

> metadata                <- as.data.frame(read.table("bamFiles/metadata.txt",
+                                                     header=TRUE))
> metadata$Code           <- paste(metadata$Tissue,
+                                  metadata$Treatment, metadata$Genotype,
+                                  sep="_")
> metadata$ExpGroup       <- paste(substr(metadata$Tissue, 1, 1),
+                                  metadata$Genotype,
+                                  substr(metadata$Treatment, 1, 1),
+                                  substr(metadata$Time, 1, 1), sep='-')
> metadata$biolrep        <- c(1,1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9,9,
+                              10,10,11,11,12,12,13,13,14,14,15,15,
+                              16,16,17,17,18,18,19,19,20,20,21,21,
+                              22,22,23,23,24,24,25,25,26,26,27,27,
+                              28,28,29,29,30,30,31,31)
> #extract out the names of the BAM files
> #to the bamFls variable (character variable)
> bamFls                  <- list.files(path="bamFiles/",
+                                       pattern="bam$", full=TRUE)
> #ordering according to the order of metadata
> bamFls_ordered          <- bamFls[match(metadata$Filename,
+                                         (sub("\\..*", "", basename(bamFls))))]
> #Assigning names to each data in bamFls_ordred using a
> #section of its basename
> names(bamFls_ordered)   <- sub("\\..*", "",
+                                basename(bamFls_ordered))

> source("https://bioconductor.org/biocLite.R")
> biocLite(c("GenomicAlignments", “GenomicFeatures”))
> library("GenomicAlignments")
> library("GenomicFeatures")

> txdb                    <- makeTxDbFromGFF(
+                            file="TAIR10_GFF3_genes.gff",
+                            format="gff",
+                            dataSource="TAIR",
+                            organism="Arabidopsis thaliana")
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK

This is where the issue begins. I managed to count the reads for each gene the first time without an issue using the following code.

> gnModel               <- transcriptsBy(txdb, "gene")
> options(mc.cores=1)
> ine_counts            <- summarizeOverlaps(gnModel,
+                          bamFls_ordered, ignore.strand=FALSE,
+                          singleEnd=FALSE,
+                          mode="Union")

This seem to function properly. However, the first attempt to create the count matrix gave me  an “acceptable” count matrix (completed in December 2015) with reads spanning up to ~1200 for the overexpressed transgene (attached the counts barplot as 1.png).

As I had a little time on my hand until the second dataset arrives, I tried to reproduce the data (in Jan 2016) using the same codes I used previously. However, now the count matrix looked very different compared to the first attempt (has very low numbers). I have attached the barplot for the counts of my overexpressed transgene from this count matrix as 2.png.

I tried the following to see whether it is a problem with my gff file, bam files, metadata file, etc.

> counter <- function(fl, gnModel) {
+     aln           <- readGAlignments(fl)
+     strand(aln)   <- "*" # for strand-blind sample prep protocol
+     hits          <- countOverlaps(aln, gnModel)
+     counts        <- countOverlaps(gnModel, aln[hits==1])
+     names(counts) <- names(gnModel)
+     counts
+ }
> counts            <- sapply(bamFls_ordered, counter, gnModel)

And the resulting count matrix was much similar to the initial count matrix I produced (attached the barplot for counts for the overexpressed transgene as 3.png). Of course there are higher numbers as the double counting is allowed.

Therefore, I concluded that (I might be wrong) it has got to do something with the summarizeOverlaps() function or the GenomicFetures() package.  I tried looking at the information from updates. However, I was unable to locate anything related to this. I prefer to use the summarizeOverlaps() function, as it lets me define whether my data is stranded/unstranded and paired/single end. And also the double counting is avoided.

If you find a free time, please advise me on where I might be going wrong, or please direct me to a place where I can find a solution.  I  pondered on it for nearly 2 weeks now, and cannot seem to find a way out. I’m sure it is a very naïve issue that has just slipped out of my sight. I would really appreciate your help.

Thanks you so much for your time.

Yours sincerely,

Shani A.

Hi Shari,

Sorry for the delayed response - this one slipped by me.

The first thing I notice is that in the 'ine_counts' example ignore.strand is FALSE and in the counter() function you change strand to '*' which is equivalent to ignore.strand=TRUE. When ignore.strand=TRUE you will (generally) see higher counts.

Another potential issue is treatment of single vs paired end reads. If your bam files are paired end you can use summarizeOverlaps() with 'single.end=FALSE'. For a rough comparison readGAlignmentPairs() followed by countOverlaps(). readGAlignments() reads data in as single-end reads. For example, if you read in paired-end data with readGAlignements() you will get roughly double the number of records than if you read the data in with readGAlignmentPairs(). When performing overlap operations, the GAlignmentPairs container only counts one hit per pair not one hit for each read in the pair.

Also there's no need to set mc.cores=1; instead use BPPARAM=MulticoreParam(1). The BPPARAM argument is passed to bplapply() called from summarizedOverlaps(). I see this isn't well documented on the man page and I will fix that.

I would recommend testing with a single bam files for starters. See if you can't get your expected counts with the single file and then go ahead and process the list.

Valerie