Answers.Body: <p>This seems relatively complicated given the structure of a BSGenome object.</p>
  
  <p>The creator of the package answered this question previously on the Bioconductor support forums:
  <a href="https://support.bioconductor.org/p/86665/#86757" rel="nofollow noreferrer">https://support.bioconductor.org/p/86665/#86757</a></p>
  
  <blockquote>
    <p>We don't provide an easy way to inject arbitrary SNPs in an arbitrary
    BSgenome at the moment. However, it should not be too hard to forge
    the BSgenome ... package. First you would need to
    compute the sequences of the mutated chromosomes (you can
    use replaceLetterAt for this), then write them to a 2bit file (put
    them in a DNAStringSet object and call rtracklayer::export on it),
    then use that 2bit file to forge the BSgenome ....
    package (see the BSgenomeForge vignette in the BSgenome package for
    how to do this).
    - Herve Pages</p>
  </blockquote>
  
  <p>Here is another similar question with more info about replaceLetterAt: <a href="https://support.bioconductor.org/p/26199/" rel="nofollow noreferrer">https://support.bioconductor.org/p/26199/</a></p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-11T16:18:45.987Z
Answers.Text: I noticed those. Question: can you think of another algorithmically efficient way to accomplish this? One way would be to convert the GRanges into a BED, then used this to create a new FASTA with the SNPS, and loaded this into BSGenome and then used `getSeq()`. That's very complicated though
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-15T17:00:26.140Z
Answers.Text: I should mention that this is a section in the official documentation---the creators appears to have "not gotten around" to implementing this. See " Injecting known SNPs in the chromosome sequences" here: https://bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/GenomeSearching.pdf
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-11T13:12:29.063Z
Answers.LastActivityDate: 2017-08-11T13:12:29.063Z
Answers.UserDisplayName: story
Body: <p>Let's say I have the genome <code>hg19</code> loaded into R via BSGenome</p>
  
  <pre><code>library("BSgenome")    
  hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)   
  </code></pre>
  
  <p>I then have a list of SNPs loaded as a GRanges object, <code>gr</code></p>
  
  <pre><code>library(GenomicRanges)
  
   &gt; gr
   GRanges object with 212 ranges and 3 metadata columns:                                                                                                                                                                                                                                                                                                                                                
          seqnames               ranges strand |     width     REF    ALT                                                                                                                                                                                                                                                                                                                             
             &lt;Rle&gt;          v &lt;IRanges&gt;  &lt;Rle&gt; | &lt;numeric&gt;                                                                                                                                                                                                                                                                                                                        
      [1]        1 [86099032, 86099032]      * |         1     C      T                                                                                                                                                                                                                                                                                                                  
      [2]        1 [86099033, 86099033]      * |         1     C      A                                                                                                                                                                                                                                                                                                                   
      [3]        1 [86099199, 86099199]      * |         1     G      A  
    ....
  </code></pre>
  
  <p>Is there a straightforward way to inject these SNPs into <code>hg19genome</code>?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-11T10:57:26.930Z
Comments.Text: By inject you mean that the hg19genome should have the alternative base instead of the reference base for those positions?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-11T16:06:15.923Z
Comments.Text: @Llopis Yes, that's what I mean
Comments.UserDisplayName: ShanZhengYang
Tags: <r><bioconductor><reference-genome><snp><biostrings>
Title: Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?
ThreadText: # Given a Genomic Ranges of SNPs, how to inject these SNPs in genome via BSGenome?
  
  (Asked by: username_0 on 2017-08-10)
  
  Let's say I have the genome `hg19` loaded into R via BSGenome
  
  ```
  library("BSgenome")    
  hg19genome = getBSgenome('BSgenome.Hsapiens.UCSC.hg19', masked=FALSE)   
  
  ```
  
  I then have a list of SNPs loaded as a GRanges object, `gr`
  
  ```
  library(GenomicRanges)
  
   > gr
   GRanges object with 212 ranges and 3 metadata columns:                                                                                                                                                                                                                                                                                                                                                
          seqnames               ranges strand |     width     REF    ALT                                                                                                                                                                                                                                                                                                                             
             <Rle>          v <IRanges>  <Rle> | <numeric>                                                                                                                                                                                                                                                                                                                        
      [1]        1 [86099032, 86099032]      * |         1     C      T                                                                                                                                                                                                                                                                                                                  
      [2]        1 [86099033, 86099033]      * |         1     C      A                                                                                                                                                                                                                                                                                                                   
      [3]        1 [86099199, 86099199]      * |         1     G      A  
    ....
  
  ```
  
  Is there a straightforward way to inject these SNPs into `hg19genome`?
  
  **Tags:** r, bioconductor, reference-genome, snp, biostrings
  
  > **Comments:**
  > **username_1**: By inject you mean that the hg19genome should have the alternative base instead of the reference base for those positions?  
  > **username_0**: @Llopis Yes, that's what I mean  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2017-08-11. Score: 3 | Accepted answer
  
  This seems relatively complicated given the structure of a BSGenome object.
  
  The creator of the package answered this question previously on the Bioconductor support forums:
  <https://support.bioconductor.org/p/86665/#86757>
  
  > We don't provide an easy way to inject arbitrary SNPs in an arbitrary
  > BSgenome at the moment. However, it should not be too hard to forge
  > the BSgenome ... package. First you would need to
  > compute the sequences of the mutated chromosomes (you can
  > use replaceLetterAt for this), then write them to a 2bit file (put
  > them in a DNAStringSet object and call rtracklayer::export on it),
  > then use that 2bit file to forge the BSgenome ....
  > package (see the BSgenomeForge vignette in the BSgenome package for
  > how to do this).
  > - Herve Pages
  
  Here is another similar question with more info about replaceLetterAt: <https://support.bioconductor.org/p/26199/>
  
  > **Comments:**
  > **username_0**: I noticed those. Question: can you think of another algorithmically efficient way to accomplish this? One way would be to convert the GRanges into a BED, then used this to create a new FASTA with the SNPS, and loaded this into BSGenome and then used `getSeq()`. That's very complicated though  
  > **username_0**: I should mention that this is a section in the official documentation---the creators appears to have "not gotten around" to implementing this. See " Injecting known SNPs in the chromosome sequences" here: https://bioconductor.org/packages/devel/bioc/vignettes/BSgenome/inst/doc/GenomeSearching.pdf  
  
  ---
  
  