QuasR: how to use an indexed reference genome?
0
Entering edit mode
Paul Shannon • 750
@paul-shannon-5161
Last seen 6.5 years ago
I am new to QuasR, and alos quite new to aligning short reads to reference genomes more generally. I cannot figure out how to use a pre-built indexed reference genome file with QuasR. The examples supplied with the package work nicely. Scaling up to using all of hg19 raises problems for me. I apologize if I am missing the obvious. To illustrate the problem, I call QuasR's qAlign method with just two arguments (quoting from the man page): sampleFile: a text file listing input sequence files and sample names genome: the reference genome for primary alignments, one of: * a string referring to a "BSgenome" package (e.g. ""BSgenome.Hsapiens.UCSC.hg19""), which will be downloaded automatically from Bioconductor if not present * the name of a fasta sequence file containing one or several sequences (chromosomes) to be used as a reference QuasR apparently invokes the bowtie indexing program when supplied either of the two "genome" options: a BSgenome package, or a fasta file. But since indexing takes a long time -- hours, apparently -- I hoped to provide a ready-made index file, and found some described here: http://bowtie-bio.sourceforge.net/tutorial.shtml specifically ftp://ftp.ccb.jhu.edu/pub/data/bowtie_indexes/hg19.ebwt.zip Various attempts to specify this file, or any of its contents (unzipped) to QuasR fail with these messages: Error: The specified genome /Users/pshannon/s/data/public/bowtie/indexes/hg19.1.ebwt does not have the extension of a fasta file (fa,fasta,fna)> Error: The specified genome has to be a file and not a directory: /Users/pshannon/s/data/public/bowtie/indexes I'll be grateful for advice on how to do this properly. Thanks, - Paul > sessionInfo() R version 3.0.0 (2013-04-03) Platform: x86_64-apple-darwin10.8.0 (64-bit) locale: [1] C attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] Rsamtools_1.13.14 BSgenome_1.29.0 Biostrings_2.29.2 QuasR_1.1.4 GenomicRanges_1.13.12 XVector_0.1.0 [7] IRanges_1.19.8 BiocGenerics_0.7.2 Rbowtie_1.1.3 BiocInstaller_1.11.1 loaded via a namespace (and not attached): [1] AnnotationDbi_1.23.11 Biobase_2.21.2 DBI_0.2-7 GenomicFeatures_1.13.8 RCurl_1.95-4.1 [6] RSQLite_0.11.3 ShortRead_1.19.3 XML_3.95-0.2 biomaRt_2.17.0 bitops_1.0-5 [11] compiler_3.0.0 grid_3.0.0 hwriter_1.3 lattice_0.20-15 rtracklayer_1.21.5 [16] stats4_3.0.0 tools_3.0.0 zlibbioc_1.7.0
BSgenome BSgenome genomes QuasR BSgenome BSgenome genomes QuasR • 1.8k views
0
Entering edit mode
Last seen 7 days ago
Switzerland
0
Entering edit mode
0
Entering edit mode
On 05/16/2013 11:41 PM, Michael Stadler wrote: > If more people would prefer to use pre-built indices, and if the BioC > agrees to host such big packages, it would also be conceivable to > provide index packages corresponding to the BSgenome packages for download. We've been experimenting with AnnotationHub (the package, and the cloud resource -- http://bioconductor.org/packages/release/bioc/html/AnnotationHub.html) as a way to provide large(r) data without having to go through the full formalism of creating packages; perhaps the indexes are a good case for this? The user or software would then library(AnnotationHub) hub = AnnotationHub() hub$<tab> ## the preceeding works, now for hypothetical...! hub$bioconductor.QuasR.BSgenome<tab> and when hitting <return> on a complete name the AnnotationHub would do the right thing, downloading and installing the index in an appropriate place (including, e.g., caching in the AnnotationHub cache, if that's appropriate...) Martin -- Computational Biology / Fred Hutchinson Cancer Research Center 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109 Location: Arnold Building M1 B861 Phone: (206) 667-2793
0
Entering edit mode
Dear Martin, Thank you for pointing out AnnotationHub - it seems that this could indeed be a way to host pre-generated alignment indices. We have internally discussed the option of providing pre-built indices to the community and came to the conclusion that with our current resources, we would not be able to create and maintain the various flavours of indices for current and future BSgenome packages. In addition to the "standard" alignment index, there are two more indices for bisulfite alignments per genome. Furthermore, for each additional aligner supported by QuasR in the future, another set of indices will have to be added. Given that QuasR can generate and re-use those indices locally, we do not consider this a priority. Michael and Dimos On 18.05.2013 02:31, Martin Morgan wrote: > On 05/16/2013 11:41 PM, Michael Stadler wrote: >> If more people would prefer to use pre-built indices, and if the BioC >> agrees to host such big packages, it would also be conceivable to >> provide index packages corresponding to the BSgenome packages for >> download. > > We've been experimenting with AnnotationHub (the package, and the cloud > resource -- > http://bioconductor.org/packages/release/bioc/html/AnnotationHub.html) > as a way to provide large(r) data without having to go through the full > formalism of creating packages; perhaps the indexes are a good case for > this? The user or software would then > > library(AnnotationHub) > hub = AnnotationHub() > hub$<tab> > ## the preceeding works, now for hypothetical...! > hub$bioconductor.QuasR.BSgenome<tab> > > and when hitting <return> on a complete name the AnnotationHub would do > the right thing, downloading and installing the index in an appropriate > place (including, e.g., caching in the AnnotationHub cache, if that's > appropriate...) > > Martin > -- -------------------------------------------- Michael Stadler, PhD Head of Computational Biology Friedrich Miescher Institute Basel (Switzerland) Phone : +41 61 697 6492 Fax : +41 61 697 3976 Mail : michael.stadler at fmi.ch
0
Entering edit mode
Last seen 7 days ago
Switzerland
Hi Paul, Please find my answers below: On 17.05.2013 15:01, Paul Shannon wrote: > Hi Michael, > > Thanks for your quick and clarifying response. > > Since it is not possible to use pre-built indices from the bowtie developers, I would be glad to have a small recipe (perhaps featured prominently in the vignette?) which > > 1) explains the need for custom-built indices > 2) provides (perhaps) a standalone QuasR-specific command for creating one It is actually much simpler than you expect: qAlign() creates the index automatically if it does not yet exist. The index is then saved in a default location (as a new R package if your reference is a BSgenome, or else in the same directory containing the fasta reference), and will be automatically re-used when qAlign is called with the same reference. > My somewhat fuzzy grasp of the current approach is that > > 1) QuasR sees the string "BSgenome.Hsapiens.UCSC.hg19" on a call to qAlign > 2) QuasR then spends a few hours building a new package with the proper index Yes, this is described in section 5.3 of the vignette. > 3) and saves this package somewhere (I could not figure out where) This is described in the documentation to qAlign. I agree that it would be better to have this all described more explicitly in a single place, so I added a description to the qAlign documentation (available shortly in the development branch). I hope this makes it all clear. Best, Michael
0
Entering edit mode
0
Entering edit mode
0
Entering edit mode