[BioC] MEDIPS.createSet error
lianoglou.steve at gene.com
Mon Mar 31 21:20:09 CEST 2014
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:
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> 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.
More information about the Bioconductor