Answers.Body: <p>There's a <code>vcf2fq</code> <a href="https://github.com/lh3/samtools/blob/master/bcftools/vcfutils.pl" rel="nofollow noreferrer">sub-program</a> that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote <a href="https://github.com/gringer/bioinfscripts/blob/master/vcf2fq.pl" rel="nofollow noreferrer">a modification</a> to implement INDEL correction as well:</p>
  
  <pre><code>./vcf2fq.pl -f &lt;input.fasta&gt; &lt;all-site.vcf&gt; &gt; &lt;output.fastq&gt;
  </code></pre>
  
  <p>It adds and removes sequence from the reference, as specified by the INDEL information in the VCF file. The code keeps a record of where in the reference sequence the INDELs should go, then inserts (or deletes) the INDELs, updating the reference position as necessary. It should work with things like deletions right next to insertions, but hasn't been extensively tested.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-08T17:46:06.227Z
Answers.Text: Could you explain how your perl script with the INDEL implementation works? i.e. what changes did you make for `vcfutils.pl` to work properly
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-08T19:18:37.573Z
Answers.Text: "It should work with things like deletions right next to insertions, but hasn't been extensively tested." If you do benchmark this, let me know---this would be interesting to see
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-09T16:50:06.840Z
Answers.Text: Another question: have you seen `vcfutils.pl` written in Python/R? I suspect these exist somewhere, as I don't want to re-invent the wheel. Is this worth a separate SO question? It would be good if your perl script was more advertised IMHO
Answers.UserDisplayName: ShanZhengYang
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-09T20:37:50.107Z
Answers.Text: I think an answer to that fits in with this question. The question is about manipulating FASTA files, not about vcfutils.
  
  No, I haven't seen any python/R alternatives. I don't see much of a point in them (certainly not a python version), given that the command-line functionality already exists via vcf2fq.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-11-09T20:40:38.147Z
Answers.Text: You're welcome to have a go at advertising this further. I've brought it up a few times on SeqAnswers, github, and I think reddit; basically every time I see a question asking something like this.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-06T21:23:14.093Z
Answers.LastActivityDate: 2017-11-10T05:49:52.723Z
Answers.LastEditDate: 2017-11-10T05:49:52.723Z
Answers.UserDisplayName: gringer
Answers.Body: <p>GATK has a solution that might work for you:
  <a href="https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_fasta_FastaAlternateReferenceMaker.php" rel="nofollow noreferrer">FastaAlternateReferenceMaker</a>, which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records." </p>
  
  <blockquote>
    <p><strong>Input</strong></p>
    
    <p>The reference, requested intervals, and any number of variant ROD files.</p>
    
    <p><strong>Output</strong></p>
    
    <p>A FASTA file representing the requested intervals.</p>
    
    <p><strong>Usage example</strong></p>
  
  <pre><code>java -jar GenomeAnalysisTK.jar \
     -T FastaAlternateReferenceMaker \
     -R reference.fasta \
     -o output.fasta \
     -L input.intervals \
     -V input.vcf \
  </code></pre>
  </blockquote>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-06T20:38:17.907Z
Answers.LastActivityDate: 2017-11-13T22:48:19.573Z
Answers.LastEditDate: 2017-11-13T22:48:19.573Z
Answers.UserDisplayName: Greg
Answers.Body: <p>Here's an update how this work in 2019. (The answer is quite similar to my update in this <a href="https://bioinformatics.stackexchange.com/a/6818/1310">thread</a>):</p>
  
  <p><code>bgzip</code> and index your <code>vcf</code> file.</p>
  
  <pre><code><span class="math-container">$ bgzip -c input.vcf &gt; input.vcf.gz
  $</span> tabix input.vcf.gz
  </code></pre>
  
  <p>Now use <a href="http://www.htslib.org/doc/bcftools.html#consensus" rel="nofollow noreferrer">bcftools consensus</a>:</p>
  
  <pre><code>$ bcftools consensus -f genome.fa input.vcf.gz &gt; consensus.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2019-01-15T20:10:30.333Z
Answers.LastActivityDate: 2019-01-15T20:10:30.333Z
Answers.UserDisplayName: finswimmer
Answers.Body: <p>You could <a href="http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html" rel="nofollow noreferrer">convert VCF to BED</a> via <code>vcf2bed --snvs</code>, <code>vcf2bed --insertions</code>, and <code>vcf2bed --deletions</code>, and then use <code>samtools faidx</code> by way of a <a href="https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640" rel="nofollow noreferrer">wrapper script</a> to convert BED to FASTA, e.g.:</p>
  
  <pre><code>$ vcf2bed --snvs &lt; variants.vcf | bed2faidxsta.pl &gt; snvs.fa
  $ vcf2bed --insertions &lt; variants.vcf | bed2faidxsta.pl &gt; insertions.fa
  $ vcf2bed --deletions &lt; variants.vcf | bed2faidxsta.pl &gt; deletions.fa
  </code></pre>
  
  <p>You need FASTA files for your reference genome, which have been indexed with <code>samtools faidx</code>, e.g., for <code>hg38</code>:</p>
  
  <pre><code>$ cd /foo/bar/baz
  $ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
  $ for fn in `ls *.fa.gz`; do gunzip $fn; done
  $ for fn in `ls *.fa`; do samtools faidx $fn; done
  </code></pre>
  
  <p>Once you have indexed FASTA files somewhere on your file system, you can pipe BED to the <code>bed2faidxsta.pl</code> script, to get out FASTA sequences.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-08-06T21:47:56.003Z
Answers.LastActivityDate: 2017-08-07T16:41:45.307Z
Answers.LastEditDate: 2017-08-07T16:41:45.307Z
Answers.UserDisplayName: Alex Reynolds
Body: <p>I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.</p>
  
  <p>The problem is, I don't know of an algorithmically sound way to do this.</p>
  
  <ol>
  <li><p>Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)? </p></li>
  <li><p>What about if I wanted to do the same with a bedpe of germline structural variants? </p></li>
  </ol>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-06T19:48:37.557Z
Comments.Text: Could you provide an example of what you're looking for?  I'm having a hard time picturing what you want.
Comments.UserDisplayName: Greg
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-06T20:11:08.953Z
Comments.Text: @Greg The above might be unclear. I have a reference FASTA. I have variants. I would now like to manipulate the reference FASTA to include these variants. The VCF tells us there's a SNP at XX location; I would like to like a FASTA that includes these SNPs. Does this make sense?
Comments.UserDisplayName: ShanZhengYang
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-08-06T20:36:18.050Z
Comments.Text: Oh I understand now.  Here's a link to a similar biostars question that has some answers that look good.  https://www.biostars.org/p/6553/
Comments.UserDisplayName: Greg
Tags: <sam><fasta><vcf><reference-genome><variation>
Title: How to manipulate a reference FASTA or bam to include variants from a VCF?
ThreadText: # How to manipulate a reference FASTA or bam to include variants from a VCF?
  
  (Asked by: username_0 on 2017-08-06)
  
  I have some software which takes fastas as the input. I need to include SNVs and InDels from a VCF into the reference hg38 and then use this.
  
  The problem is, I don't know of an algorithmically sound way to do this.
  
  1. Are there any existing software packages which could do this efficiently? Is it easier to output a FASTA, or a bam (and then convert to a FASTA)?
  2. What about if I wanted to do the same with a bedpe of germline structural variants?
  
  **Tags:** sam, fasta, vcf, reference-genome, variation
  
  > **Comments:**
  > **username_1**: Could you provide an example of what you're looking for?  I'm having a hard time picturing what you want.  
  > **username_0**: @username_1 The above might be unclear. I have a reference FASTA. I have variants. I would now like to manipulate the reference FASTA to include these variants. The VCF tells us there's a SNP at XX location; I would like to like a FASTA that includes these SNPs. Does this make sense?  
  > **username_1**: Oh I understand now.  Here's a link to a similar biostars question that has some answers that look good.  https://www.biostars.org/p/6553/  
  
  ---
  
  ## 4 answers
  
  ### Answer by username_4 on 2017-08-06. Score: 4 | Accepted answer
  
  You could [convert VCF to BED](http://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/vcf2bed.html) via `vcf2bed --snvs`, `vcf2bed --insertions`, and `vcf2bed --deletions`, and then use `samtools faidx` by way of a [wrapper script](https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640) to convert BED to FASTA, e.g.:
  
  ```
  $ vcf2bed --snvs < variants.vcf | bed2faidxsta.pl > snvs.fa
  $ vcf2bed --insertions < variants.vcf | bed2faidxsta.pl > insertions.fa
  $ vcf2bed --deletions < variants.vcf | bed2faidxsta.pl > deletions.fa
  
  ```
  
  You need FASTA files for your reference genome, which have been indexed with `samtools faidx`, e.g., for `hg38`:
  
  ```
  $ cd /foo/bar/baz
  $ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
  $ for fn in `ls *.fa.gz`; do gunzip $fn; done
  $ for fn in `ls *.fa`; do samtools faidx $fn; done
  
  ```
  
  Once you have indexed FASTA files somewhere on your file system, you can pipe BED to the `bed2faidxsta.pl` script, to get out FASTA sequences.
  
  
  ---
  
  ### Answer by username_1 on 2017-08-06. Score: 7
  
  GATK has a solution that might work for you:
  [FastaAlternateReferenceMaker](https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_fasta_FastaAlternateReferenceMaker.php), which : "Given a variant callset, this tool replaces the reference bases at variation sites with the bases supplied in the corresponding callset records."
  
  > **Input**
  >
  > The reference, requested intervals, and any number of variant ROD files.
  >
  > **Output**
  >
  > A FASTA file representing the requested intervals.
  >
  > **Usage example**
  >
  > ```
  > java -jar GenomeAnalysisTK.jar \
  >    -T FastaAlternateReferenceMaker \
  >    -R reference.fasta \
  >    -o output.fasta \
  >    -L input.intervals \
  >    -V input.vcf \
  >
  > ```
  
  
  ---
  
  ### Answer by username_2 on 2017-08-06. Score: 3
  
  There's a `vcf2fq` [sub-program](https://github.com/lh3/samtools/blob/master/bcftools/vcfutils.pl) that was written as part of vcfutils to convert a VCF file into a fastq file given a reference sequence. Unfortunately this doesn't work properly with INDELs (it will just mask them, rather than actually converting them), so I wrote [a modification](https://github.com/username_2/bioinfscripts/blob/master/vcf2fq.pl) to implement INDEL correction as well:
  
  ```
  ./vcf2fq.pl -f <input.fasta> <all-site.vcf> > <output.fastq>
  
  ```
  
  It adds and removes sequence from the reference, as specified by the INDEL information in the VCF file. The code keeps a record of where in the reference sequence the INDELs should go, then inserts (or deletes) the INDELs, updating the reference position as necessary. It should work with things like deletions right next to insertions, but hasn't been extensively tested.
  
  > **Comments:**
  > **username_0**: Could you explain how your perl script with the INDEL implementation works? i.e. what changes did you make for `vcfutils.pl` to work properly  
  > **username_0**: "It should work with things like deletions right next to insertions, but hasn't been extensively tested." If you do benchmark this, let me know---this would be interesting to see  
  > **username_0**: Another question: have you seen `vcfutils.pl` written in Python/R? I suspect these exist somewhere, as I don't want to re-invent the wheel. Is this worth a separate SO question? It would be good if your perl script was more advertised IMHO  
  > **username_2**: I think an answer to that fits in with this question. The question is about manipulating FASTA files, not about vcfutils.  No, I haven't seen any python/R alternatives. I don't see much of a point in them (certainly not a python version), given that the command-line functionality already exists via vcf2fq.  
  > **username_2**: You're welcome to have a go at advertising this further. I've brought it up a few times on SeqAnswers, github, and I think reddit; basically every time I see a question asking something like this.  
  
  ---
  
  ### Answer by username_3 on 2019-01-15. Score: 2
  
  Here's an update how this work in 2019. (The answer is quite similar to my update in this [thread](https://bioinformatics.stackexchange.com/a/6818/1310)):
  
  `bgzip` and index your `vcf` file.
  
  ```
  $ bgzip -c input.vcf > input.vcf.gz
  $ tabix input.vcf.gz
  
  ```
  
  Now use [bcftools consensus](http://www.htslib.org/doc/bcftools.html#consensus):
  
  ```
  $ bcftools consensus -f genome.fa input.vcf.gz > consensus.fa
  
  ```
  
  
  ---
  
  