[BioC] MEDIPS.createSet error

[Vining, Kelly]

Hi,
Thanks for the advice thus far! To confirm what is in my BSgenome variable, I did this:

> class(BSgenome)
[1] "BSgenome"
attr(,"package")
[1] "BSgenome"
> names(BSgenome)
 [1] "Chr01" "Chr02" "Chr03" "Chr04" "Chr05" "Chr06" "Chr07" "Chr08" "Chr09"
[10] "Chr10" "Chr11" "Chr12" "Chr13" "Chr14" "Chr15" "Chr16" "Chr17" "Chr18"
[19] "Chr19" "scaf"

And then:
> print(BSgenome)
Black cottonwood genome
|
| organism: Populus trichocarpa (Black cottonwood)
| provider: Phytozome (JGI)
| provider version: 3.0
| release date: January 2010
| release name: Populus trichocarpa v3.0
|
| single sequences (see '?seqnames'):
|   Chr01  Chr02  Chr03  Chr04  Chr05  Chr06  Chr07  Chr08  Chr09  Chr10  Chr11
|   Chr12  Chr13  Chr14  Chr15  Chr16  Chr17  Chr18  Chr19
|
| multiple sequences (see '?mseqnames'):
|   scaf
|
| (use the '$' or '[[' operator to access a given sequence)


So that looks ok. Interestingly, when I followed the vignette and did the equivalent of 
BSgenome="BSgenome.Hsapiens.UCSC.hg19"

That didn't work for me. It only worked without quotes. If I included quotes, it just assigned a character vector to that variable.

Then, following your advice: 

> si.bsg <- seqinfo(BSgenome.Ptrichocarpa.Phytozome.v3)
> si.bam <- seqinfo(BamFile("FallBud_brep1_aln_sorted.bam"))
Error in seqinfo(BamFile("FallBud_brep1_aln_sorted.bam")) :
  error in evaluating the argument 'x' in selecting a method for function 'seqinfo': Error: could not find function "BamFile"

Hmmm....

Continuing down the rabbit hole...

--Kelly
________________________________________
From: mailinglist.honeypot at gmail.com [mailinglist.honeypot at gmail.com] on behalf of Steve Lianoglou [lianoglou.steve at gene.com]
Sent: Monday, March 31, 2014 12:20 PM
To: Vining, Kelly
Cc: Lukas Chavez; bioconductor at r-project.org
Subject: Re: [BioC] MEDIPS.createSet error

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