Leo -- thanks for your recent help with fullCoverage() in derfinder. After your updates, that function began working very well for me. However, I have discovered that analyzeChr() is still having issues with chromosomes > 22. When I give it a chromosome in the human range, it does fine, but once we get up into higher numbers, it breaks.

R 3.1.1
Bioconductor version 3.0 (BiocInstaller 1.16.0)
derfinder 1.0.4
GenomeInfoDb 1.2.2 (version including the Canis familiaris genome update)

My script:

library(derfinder)
library(GenomeInfoDb)
library(GenomicRanges)

args <- commandArgs(trailingOnly = TRUE)

perms = 200  # for alpha = 0.05
cores = 15   # too many -> memory problems
#perms = 1
#cores = 25
chroms = c(args[1])

files <- rawFiles(datadir="tophat", samplepatt="tophat-juncs-*")

# Set up our groups
tameIds <- read.table("/home/hekman2/Dropbox/UIUC/transcriptome/hypothalamus/notes/samples-hyp-tame.txt")
tameFileNames <- lapply(tameIds, function(id) { return(paste("tophat-juncs-", id, sep="")) } )
groups <- as.factor(sapply(names(files), function(filename) { ifelse(grepl(filename, tameFileNames), "tame", "aggr") }))

# Initial data analysis over all chromosomes:
print("Starting fullCoverage...")
fullCov <- fullCoverage(files, chrs = chroms, species = 'canis_familiaris')
sampleDepths <- sampleDepth(collapseFullCoverage(fullCov), 1)
models <- makeModels(sampleDepths, testvars = groups)

# Keep only bases with coverage > 2
print("Filtering...")
filteredCov <- lapply(fullCov, filterData, cutoff = 2)

originalWd <- getwd()
setwd(file.path(originalWd, 'derResults'))

# Analyze each chromosome individually (analysis will be parallelized
# and significance will be evaluated by permutations).
for (chrom in chroms) {
  print(chrom)
  currentChrom <- paste("chr", chrom, sep="")
  analyzeChr(chr = currentChrom, filteredCov[[currentChrom]], models, mc.cores=cores, groupInfo = groups, writeOutput = TRUE, cutoffFstat = 0.05, nPermute = perms, seeds = 19731107 + seq_len(perms))
}

print("Chromosome analysis completed.")

The error on chromosomes > 22:

2014-11-02 12:44:36 analyzeChr: Annotating regions
Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) :
  subscript contains NAs or out-of-bounds indices
Calls: analyzeChr ... extractROWS -> normalizeSingleBracketSubscript -> NSBS -> NSBS
In addition: Warning message:
In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
Execution halted

I also noticed this output message during execution of analyzeChr():

nearestgene: loading bumphunter hg19 transcript database

I'm guessing we should not be using hg19. In fact, this makes  me wonder if the other low-numbered chromosomes which completed correctly have valid results or if they need to be re-run against a different transcript database.

Any help you can give is much appreciated. I suspect that the only solution will involve you editing the code again, for which I'm sorry if it's true!

Jessica

Hi Jessica,

There are two different issues in this. First, I had not updated the BioC release version of `derfinder` to match the changes in 1.2.1 from GenomeInfoDb. That has been fixed already. Basically, I wasn't sure if they had merged the changes from GenomeInfoDb Bioc-devel version 1.3.2 to Bioc-release.

The second issue is that `derfinder` relies on `bumphunter` http://www.bioconductor.org/packages/release/bioc/html/bumphunter.html to annotate the candidate DERs. `derfinder` uses `bumphunter::annotateNearest()` and this function currently only supports Homo sapiens build hg19. You will need to generate your own annotation object for Canis familiaris and it should look similar to what you get running:

> library(bumphunter)
> data(TT, package = "bumphunter", envir = environment())
> TT

I highly recommend creating a new question with the `bumphunter` tag on how to proceed with it. Meanwhile, you can still use `derfinder` to find the DERs by using `analyzeChr(runAnnotation = FALSE)`. Note that your previous results for chrs 1 to 19 were annotated versus Homo sapiens because the defaults for `analyzeChr()` are `subject = 'hg19', runAnnotation = TRUE`, so the annotation results will not make sense. You don't have to re-run those chromosomes, simply remove the annotation information.

The following code  works and generates  Note that I used some artificial models and cutoffs for the toy example. The BAM and BAI file used are those you previously posted.

Used derfinder 1.0.6, derfinderHelper 1.0.3, and bumphunter 1.6.0

Best,
Leonardo
 

Great, thank you. I have posted a request for bumphunter help -- see Using bumphunter with non-human genomes.

I too was surprised to not get an email update, because I had set my options to request one -- but I'm getting email updates now so it's all good :)

I'm checking to make sure that the code you suggested actually does what it should on my end with the real BAMs. When it does, I'll accept this answer.

Thanks,

Jessica