Entering edit mode
Vining, Kelly
▴
220
@vining-kelly-5029
Last seen 10.2 years ago
Wow...this is worrisome if you guys are stumped! But...following your
most recent suggestion, I seem to have generated a new error which may
provide further insight. Getting chr_lengths seems to work, but...see
what happens when I then submit the MEDIPS.createSet command:
> library(BSgenome.Ptrichocarpa.Phytozome.v3)
> BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
> dataset=get(ls(paste("package:", BSgenome, sep="")))
> chromosomes=c("Chr01", "Chr02")
> chr_lengths=as.numeric(sapply(chromosomes,
function(x){as.numeric(length(dataset[[x]]))}))
> chr_lengths
[1] 50495391 25263035
> budFall1_MEDIP = MEDIPS.createSet(file="accepted_hits_F1.bam",
BSgenome=BSgenome, extend=300, uniq=TRUE, shift=0, window_size=1000)
Reading bam alignment accepted_hits_F1.bam
Total number of imported short reads: 25635579
Extending reads...
Creating GRange Object...
Extract unique regions...
Number of unique short reads: 19523016
Calculating genomic coordinates...Error in vector(length =
supersize_chr[length(chromosomes)], mode = "character") :
vector size cannot be NA/NaN
________________________________
From: Lukas Chavez [lukas.chavez.mailings@googlemail.com]
Sent: Friday, December 20, 2013 4:04 PM
To: Matthias Lienhard
Cc: Vining, Kelly; bioconductor@r-project.org
Subject: Re: [BioC] separate MEDIPS.createSet required for each
replicate? RE: using custom BSgenome with MEDIPS
Indeed, that also looks fine. Then, unfortunately, I do not have a
clue what's going on. What happens in the MEDIPS.createSet function is
the following:
(...)
dataset=get(ls(paste("package:", BSgenome, sep="")))
chr_lengths=as.numeric(seqlengths(dataset)[chromosomes])
(...)
where 'chromosomes' has been pre-calculated based on the set of
chromosomes in your bam files and/or by the chr.select parameter.
You could try if the following works:
library(BSgenome.Ptrichocarpa.Phytozome.v3)
BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
dataset=get(ls(paste("package:", BSgenome, sep="")))
chromosomes=c("Chr01", "Chr02")
chr_lengths=as.numeric(sapply(chromosomes,
function(x){as.numeric(length(dataset[[x]]))}))
If this works (i.e. you have the sizes of the chromosomes "Chr01",
"Chr02" at "chr_lengths"), then I do not see why there will be an
error in MEDIPS.
Can you send the whole script you execute up to the MEDIPS.createSet
function that causes the error?
Lukas
On Fri, Dec 20, 2013 at 4:01 PM, Matthias Lienhard
<lienhard@molgen.mpg.de<mailto:lienhard@molgen.mpg.de>> wrote:
Now I' stumbled. That was exactly the line that cased the error in
your function call. Does the error still occur if you call
MEDIPS.createSet in the same session where get(ls(paste ... worked?
Am 20.12.2013 15<tel:20.12.2013%2015>:53, schrieb Vining, Kelly:
OK, here's that output:
BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
library(BSgenome.Ptrichocarpa.Phytozome.v3)
get(ls(paste("package:",BSgenome, sep="")))
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)
________________________________________
From: bioconductor-bounces@r-project.org<mailto:bioconductor- bounces@r-project.org=""> [bioconductor-bounces@r-project.org<mailto :bioconductor-bounces@r-project.org="">] on behalf of Matthias Lienhard
[lienhard@molgen.mpg.de<mailto:lienhard@molgen.mpg.de>]
Sent: Friday, December 20, 2013 3:32 PM
To: bioconductor@r-project.org<mailto:bioconductor@r-project.org>
Subject: Re: [BioC] separate MEDIPS.createSet required for each
replicate? RE: using custom BSgenome with MEDIPS
OK, that looks alright. What happens when prompting
BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
library(BSgenome.Ptrichocarpa.Phytozome.v3)
get(ls(paste("package:",BSgenome, sep="")))
Am 20.12.2013 15<tel:20.12.2013%2015>:18, schrieb Vining, Kelly:
Hi Lukas,
Here's the output:
BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
ls("package:BSgenome.Ptrichocarpa.Phytozome.v3")
[1] "BSgenome.Ptrichocarpa.Phytozome.v3" "Ptrichocarpa"
get(ls("package:BSgenome.Ptrichocarpa.Phytozome.v3"))
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)
________________________________
From: Lukas Chavez [lukas.chavez.mailings@googlemail.com<mailto:lukas. chavez.mailings@googlemail.com="">]
Sent: Friday, December 20, 2013 2:11 PM
To: Vining, Kelly
Cc: Lukas Chavez;
bioconductor@r-project.org<mailto:bioconductor@r-project.org>
Subject: Re: [BioC] separate MEDIPS.createSet required for each
replicate? RE: using custom BSgenome with MEDIPS
Hi Kelly,
then there still something does not work with your custom genome.
What is your output for:
library(BSgenome.Ptrichocarpa.Phytozome.v3)
ls("package:BSgenome.Ptrichocarpa.Phytozome.v3")
get(ls("package:BSgenome.Ptrichocarpa.Phytozome.v3"))
?
Lukas
On Fri, Dec 20, 2013 at 1:52 PM, Vining, Kelly <kelly.vining@oregonsta te.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kelly.vining@oregon state.edu<mailto:kelly.vining@oregonstate.edu="">>> wrote:
Hi Lukas,
I have in my code:
library(BSgenome.Ptrichocarpa.Phytozome.v3)
BSgenome = "BSgenome.Ptrichocarpa.Phytozome.v3"
budFall2_MEDIP = MEDIPS.createSet(file="accepted_hits_F2.bam",
BSgenome=BSgenome, extend=300, uniq=TRUE, shift=0, window_size=1000)
Yes, I am seeing the same error for the first two .bam files I
submitted, so figured I'd better stop until I figured out what was
going on.
--Kelly
From: Lukas Chavez [mailto:lukas.chavez@googlemail.com<mailto:lukas.ch avez@googlemail.com=""><mailto:lukas.chavez@googlemail.com<mailto:lukas.c havez@googlemail.com="">>]
Sent: Friday, December 20, 2013 1:42 PM
To: Vining, Kelly
Cc: Lukas Chavez; bioconductor@r-project.org<mailto:bioconductor@r-pro ject.org=""><mailto:bioconductor@r-project.org<mailto:bioconductor@r-proj ect.org="">>
Subject: Re: [BioC] separate MEDIPS.createSet required for each
replicate? RE: using custom BSgenome with MEDIPS
Hi Kelly,
No, youshould not ignore the error message. I did not fully
understand, if this happens for all of your bam files or only for one
of your bam files? What do you have on your BSgenome variable, I
cannot see this in your code?
Lukas
On Fri, Dec 20, 2013 at 1:00 PM, Vining, Kelly <kelly.vining@oregonsta te.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kelly.vining@oregon state.edu<mailto:kelly.vining@oregonstate.edu="">><mailto:kelly.vining@or egonstate.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kelly.vining @oregonstate.edu<mailto:kelly.vining@oregonstate.edu="">>>> wrote:
A follow-up question about an error message with MEDIPS.createSet
command: when I load my reference genome, the create a MEDIPS set with
one of my .bam files, I get the following error message:
library(BSgenome.Ptrichocarpa.Phytozome.v3)
budFall2_MEDIP = MEDIPS.createSet(file="accepted_hits_F2.bam",
BSgenome=BSgenome, extend=300, uniq=TRUE, shift=0, window_size=1000)
Reading bam alignment accepted_hits_F2.bam
Total number of imported short reads: 46734466
Extending reads...
Creating GRange Object...
Extract unique regions...
Number of unique short reads: 36603869
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 don't know how to interpret this. Is this error generated from
GenomicRanges? Does it indicate that I've done something incorrectly,
or is it safe to ignore it?
Thanks again,
--Kelly V.
________________________________________
From: bioconductor-bounces@r-project.org<mailto:bioconductor- bounces@r-project.org=""><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-="" bounces@r-project.org="">><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-="" bounces@r-project.org=""><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-bounces@r-project.org="">>>
[bioconductor-bounces@r-project.org<mailto:bioconductor- bounces@r-project.org=""><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-="" bounces@r-project.org="">><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-="" bounces@r-project.org=""><mailto:bioconductor- bounces@r-project.org<mailto:bioconductor-bounces@r-project.org="">>>] on
behalf of Vining, Kelly [Kelly.Vining@oregonstate.edu<mailto:kelly.vin ing@oregonstate.edu=""><mailto:kelly.vining@oregonstate.edu<mailto:kelly. vining@oregonstate.edu="">><mailto:kelly.vining@oregonstate.edu<mailto:ke lly.vining@oregonstate.edu=""><mailto:kelly.vining@oregonstate.edu<mailto :kelly.vining@oregonstate.edu="">>>]
Sent: Friday, December 20, 2013 12:34 PM
To: 'Lukas Chavez'
Cc: bioconductor@r-project.org<mailto:bioconductor@r-project.org><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailt o:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">>>
Subject: [BioC] separate MEDIPS.createSet required for each replicate?
RE: using custom BSgenome with MEDIPS
Hi Lukas,
I have a question about the MEDIPS.createSet function with biological
replicates. In the vignette, a MEDIPS.createSet command is issued for
the first replicate of a set, and then other replicates are appended
to make a list of MEDIPS SETs. Does that mean that one "set" equals
one data set from one individual, so that I have to have a separate
MEDIPS.createSet call for each biological replicate? This seems
inefficient, because the BSgenome, uniq, extend etc. parameters must
be repeated for each rep even though all of those parameters are the
same across all replicates.
Can I simply do something like this?
Rep1_MeDIP = system.file("extdata", rep1.bam)
Rep2_MeDIP = system.file("extdata", rep2.bam)
Rep3_MeDIP = system.file("extdata", rep3.bam)
AllReps = MEDIPS.createSet(c(Rep1_MEDIP, Rep2_MEDIP, Rep3_MEDIP,
BSgenome=BSgenome)
Thanks,
--Kelly V.
From: Lukas Chavez [mailto:lukas.chavez.mailings@googlemail.com<mailto :lukas.chavez.mailings@googlemail.com=""><mailto:lukas.chavez.mailings@go oglemail.com<mailto:lukas.chavez.mailings@googlemail.com="">><mailto:luka s.chavez.mailings@googlemail.com<mailto:lukas.chavez.mailings@googlema="" il.com=""><mailto:lukas.chavez.mailings@googlemail.com<mailto:lukas.chave z.mailings@googlemail.com="">>>]
Sent: Wednesday, December 18, 2013 12:42 PM
To: Vining, Kelly
Cc: bioconductor@r-project.org<mailto:bioconductor@r-project.org><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailt o:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">>>
Subject: Re: [BioC] using custom BSgenome with MEDIPS
Hi Kelly,
let us assume that the package name of your successfully created
custom BSgenome for e.g. 'oryza sativa' is 'BSgenome.OryzaSativa' and
that it is installed in your R environment using R CMD INSTALL
<package>.
First load your custom BSgenome object into your R session using
library(BSgenome.OryzaSativa)
Note, BSgenome.OryzaSativa will not appear in the list returned by
available.genomes(). However, in case the library is loaded you should
be able to immediately use it as usual in MEDIPS like
library(MEDIPS)
BSgenome= "BSgenome.OryzaSativa"
MSet = MEDIPS.createSet(file="oryza.bam", BSgenome=BSgenome)
Lukas
On Wed, Dec 18, 2013 at 6:57 AM, Vining, Kelly <kelly.vining@oregonsta te.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kelly.vining@oregon state.edu<mailto:kelly.vining@oregonstate.edu="">><mailto:kelly.vining@or egonstate.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kelly.vining @oregonstate.edu<mailto:kelly.vining@oregonstate.edu="">>><mailto:kelly.v ining@oregonstate.edu<mailto:kelly.vining@oregonstate.edu=""><mailto:kell y.vining@oregonstate.edu<mailto:kelly.vining@oregonstate.edu="">><mailto: kelly.vining@oregonstate.edu<mailto:kelly.vining@oregonstate.edu=""><mail to:kelly.vining@oregonstate.edu<mailto:kelly.vining@oregonstate.edu="">>>
>> wrote:
________________________________
From: Vining, Kelly
Sent: Tuesday, December 17, 2013 10:58 AM
To: 'Lukas Chavez'
Cc: bioconductor@r-project.org<mailto:bioconductor@r-project.org><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailt o:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">>><mail to:bioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailt o:bioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mailt o:bioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailto :bioconductor@r-project.org<mailto:bioconductor@r-project.org="">>>>
Subject: RE: [BioC] Question about MEDIPS extend parameter
Hello BioC group,
have created a BSgenome object for my genome of interest, and it is in
my current working directory. I now want to use it with MEDIPS. How do
I make it appear under
available.genomes()? Is it necessary for my custom genome to appear in
this list for me to be able to use it with MEDIPS? Or is there a
different method for accessing this object?
Thanks much,
--Kelly V.
[[alternative HTML version deleted]]
_______________________________________________
Bioconductor mailing list
Bioconductor@r-project.org<mailto:bioconductor@r-project.org><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailto:bi oconductor@r-project.org<mailto:bioconductor@r-project.org="">>><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailto:bi oconductor@r-project.org<mailto:bioconductor@r-project.org="">><mailto:bi oconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailto:bio conductor@r-project.org<mailto:bioconductor@r-project.org="">>>>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]
_______________________________________________
Bioconductor mailing list
Bioconductor@r-project.org<mailto:bioconductor@r-project.org><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org="">><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org=""><mailto:bi oconductor@r-project.org<mailto:bioconductor@r-project.org="">>>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]
_______________________________________________
Bioconductor mailing list
Bioconductor@r-project.org<mailto:bioconductor@r-project.org><mailto:b ioconductor@r-project.org<mailto:bioconductor@r-project.org="">>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]
_______________________________________________
Bioconductor mailing list
Bioconductor@r-project.org<mailto:bioconductor@r-project.org>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]
_______________________________________________
Bioconductor mailing list
Bioconductor@r-project.org<mailto:bioconductor@r-project.org>
https://stat.ethz.ch/mailman/listinfo/bioconductor
Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
[[alternative HTML version deleted]]
