Hello,

I have some mouse RNA seq data I am trying to process and so far its been one issue after another.

I am using SystemPipeR to complete the workflow. Most processes are done  on a school computer cluster

Alignment with Tophat2 & Bowtie2  will not give me alignment percentage greater than  96%:

Some of the concordant pair alignment rate of these these files would be ~68%, while some are ~90%

Example Tophat execution:  /tophat2/2.1.0/bin/tophat -p 16 -o /results/022722_C3_1.fastq.gz.tophat /data/genome /data/022722_C3_1.fastq.gz /data/022722_C3_2.fastq.gz

Here is the fastq report:https://drive.google.com/open?id=0B6NkdsrSEj_wUExMSnlXOC14Unc 

Any ideas as to why I get such low alignment or how I can improve it?

Okay, besides the low rates, everything else goes fine with alignment. Next step is to get counts. I use SummarizeOverlaps to get counts from bam file using txdb made from mm10

Here is where I am stuck so far: 

Trying to do it with parallelization on the Cluster gives an error. I was able to get it to work just using the local computer processor on 2 files. But when I try it on 20 files no luck: I run the following

f <- function(x) {
    library(systemPipeR); library(BiocParallel); library(GenomicFeatures); library(AnnotationDbi)
    txdb <- loadDb("./data/genesmrna.sqlite")
    eByg <- exonsBy(txdb, by=c("gene"))
    args <- systemArgs(sysma="param/tophat.param", mytargets="targetsPE.txt")
    bfl <- BamFileList(outpaths(args), yieldSize=50000, index=character())
    summarizeOverlaps(eByg, bfl[x], mode="Union", ignore.strand=TRUE, inter.feature=TRUE, singleEnd=FALSE)
countDFeByg <- sapply(seq(along=counteByg), function(x) assays(counteByg[[x]])$counts)
}
funs <- makeClusterFunctionsTorque("torque.tmpl")
param <- BatchJobsParam(length(args), resources=list(walltime="20:00:00", nodes="1:ppn=8", memory="16gb"), cluster.functions=funs)
register(param)
counteByg <- bplapply(seq(along=args), f) 

I get the following error:

Warning in .make_GAlignmentPairs_from_GAlignments(gal, strandMode = strandMode,  :
    255004 alignments with ambiguous pairing were dumped.
    Use 'getDumpedAlignments()' to retrieve them from the dump environment.
Warning in .make_GAlignmentPairs_from_GAlignments(gal, strandMode = strandMode,  :
    516780 alignments with ambiguous pairing were dumped.
    Use 'getDumpedAlignments()' to retrieve them from the dump environment.
Result:
List of 2
 $ message: chr "object 'counteByg' not found"
 $ call   : language seq(along = counteByg)
 - attr(*, "class")= chr [1:3] "remote_error" "bperror" "condition"
 - attr(*, "traceback")= chr [1:29] "17: BatchJobs:::doJob(reg = loadRegistry(\"//BiocParallel_tmp_6c923953a7b5\"), " "        ids = c(1L), multiple.result.files = FALSE, disable.mail = FALSE, " "        first = 2L, last = 1L, array.id = NA)" "16: try(applyJobFunction(reg, job, cache), silent = TRUE)" ...

No idea what: "alignments with ambiguous pairing were dumped." means.

Any help would be greatly appreciated!

There may be subtleties here that I am missing, in which case Thomas Girke will probably be along with a better answer. In the mean time, a 75% alignment isn't particularly bad, and expecting better than 96% alignment is probably a bit optimistic on your part, so I wouldn't be particularly worried about those statistics.

Also note that you aren't getting an Error, you are getting a Warning. These are not the same thing. An error means that something went horribly wrong, and R refused to continue, sputtering in indignation. A warning is just a polite notification that something happened that you might not have expected. In this case, you are just being warned that you had a few hundred thousand alignments dumped, because the pairings were ambiguous.

Unfortunately, the warning is a bit, um, ambiguous itself, and I believe it comes from readGAlignmentPairs in the GenomicAlignments package. If you look at ?readGAlignmentPairs you can see what exactly is being tested to say if a pairing is ambiguous, which has to do with lots of 'inside baseball' about the SAM specification and whatnot. The basic idea behind 'proper pairing' is that we expect (for Illumina data, at least) for paired ends to be on opposite strands, facing inward, with an expected insert size of like 450 bp or so (these are all flags that you send to the aligner - I am giving some idea of the usual defaults). Anything that doesn't fulfill those basic criteria will be marked as not properly paired by the aligner.

Anyway, you have like 25M - 75M reads per sample, so having half a million or so dumped off is inconsequential.
 

Dear Abass, 

What makes you believe that an alignment rate around 90% is low? Especially, for PE reads this is reasonable. Even 70% is not too bad. You could certainly try to diagnose the source of this drop for some of your samples (e.g. lower quality, adapter contaminations, PCR artifacts, etc). For comparison you also may want to take a look what alignment frequencies are reported by Tophat directly. These are stored in the Tophat output directories of this naming structure: results/*.fastq.gz.tophat/align_summary.txt. Please also note, you might get better alignment frequencies with newer aligners such as the new version of Tophat called HISAT2  (param file for systemPipeR is here) or rsubread. Those aligners are more sensitive and also much more time efficient. 

As for the error in your read counting step on your computer cluster, one issue is the last line in your f() function. This line doesn't belong there. It needs to be run after the job submission with bplapply (please see vignette for this). Also, in case your computer nodes run several R versions, please make sure that your cluster template file (here torque.tmpl) loads the same R version you are using to submit the jobs (preferentially a recent version). The first messages you are seeing are just warnings from summarizeOverlaps. To understand their meaning it can be helpful to read the help file of the summarizeOverlaps function and/or the vignette of the GenomicAlignments package where this function is defined. 

Thomas