[BioC] MEDIPS.createSet error

[Steve Lianoglou]

Hi,

On Mon, Mar 31, 2014 at 12:01 PM, Vining, Kelly
<Kelly.Vining at oregonstate.edu> wrote:
> Hi Lukas,
> Well, I'm basically back to the original error I had when I started this work back in December. At that point, I'd decided that my bam files must have had chromosome or scaffold names that didn't exist in the reference genome. So I went back to the original data set and redid all of the alignments to make sure the version of the reference genome I was using was the same. It seems that I something is still not matching up.
>
> How can I compare the set of chromosome names between the bam file and the reference genome?
>
>> budFall1_MEDIP = MEDIPS.createSet(file="FallBud_brep1_aln_sorted.bam", BSgenome=BSgenome, extend=300, uniq=TRUE, shift=0, window_size=1000)
> Reading bam alignment FallBud_brep1_aln_sorted.bam
> Total number of imported short reads: 21038273
> Extending reads...
> Creating GRange Object...
> Extract unique regions...
> Number of unique short reads: 17855868
> Error in as.environment(pos) :
>   no item called "paste("package:", BSgenome, sep = "")" on the search list
> In addition: Warning message:
> In ls(paste("package:", BSgenome, sep = "")) :
>   āpaste("package:", BSgenome, sep = "")ā converted to character string

I'm not so sure that's what this error is telling you.

First check if "BSgenome" is a real thing in your workspace -- you are
passing some variable named "BSgenome" as the BSgenome parameter to
the createSet function.

Skimming the vignette, it looks your BSgenoome variable should be set
to the name of a BSgenome package to use/load.

So, in the same R session where this call is failing, what is the
output of `print(BSgenome)`.

Note how the vignette has this line:

BSgenome="BSgenome.Hsapiens.UCSC.hg19"

Which uses the hg19 genome as the reference.

Anyway, if BSgenome is set correctly for you, try to load it first to
make sure it is installed, eg:

R> library(BSgenome, character.only=TRUE)

Finally, the BSgenome that you load and the BAM file you are parsing
need concordant names for their chromosomes. You can check the
chromosome info in each by invoking the `seqinfo` method.

If I was using the "BSgenome.Hsapiens.UCSC.hg19" package, I'd do:

R> library(Rsamtools)
R> library(BSgenome.Hsapiens.UCSC.hg19)
R> si.bsg <- seqinfo(BSgenome.Hsapiens.UCSC.hg19)
R> si.bam <- seqinfo(BamFile("FallBud_brep1_aln_sorted.bam"))

Now look to see that chromosomes (seqnames) in si.bsg and si.bam are concordant.

HTH,
-steve

-- 
Steve Lianoglou
Computational Biologist
Genentech