<p><a href="https://github.com/USDA-VS/vSNP" rel="nofollow noreferrer">vSNP</a> outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  <code>vsnp_merge_vcf_into_fasta.py</code>.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  <a href="https://github.com/broadinstitute/pilon/wiki/Methods-of-Operation" rel="nofollow noreferrer">Pilon</a> is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".</p>

<pre class="lang-py prettyprint-override"><code>vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta
</code></pre>

<p>Use <code>*zc.vcf</code> output to merge into FASTA</p>

<pre class="lang-py prettyprint-override"><code>vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf
</code></pre>

CC BY-SA 4.0
2020-02-13T12:50:28.870Z
2020-02-13T12:50:28.870Z
Stuber
<p>This is similar to <a href="https://www.biostars.org/p/303439/" rel="nofollow noreferrer">this question</a>, which links to <a href="https://samtools.github.io/bcftools/howtos/consensus-sequence.html" rel="nofollow noreferrer">this solution</a>.</p>

<p>The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this <a href="https://www.htslib.org/doc/samtools-depth.html" rel="nofollow noreferrer">using</a> <code>samtools depth</code> to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:</p>

<pre><code>### Someone could come up with something better if they spent &gt;3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

threshold = 0
fasta_string = ''
for position in reference_fasta:
  if coverage_data[position] &lt;= threshold:
    fasta_string += "N"
  else:
    fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta. 
write_as_fasta(fasta_string)
</code></pre>

<p>I am a little skeptical of this approach because it seems like <strong>if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against</strong>. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.</p>

<p>I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase <code>acgt</code> with your covered regions uppercase <code>ACGT</code>, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.</p>

<p>Hope that helps-</p>

CC BY-SA 4.0
2020-02-05T19:58:07.297Z
2020-02-05T21:09:02.633Z
2020-02-05T21:09:02.633Z
Maximilian Press
<p>Have you googled <em>de novo assembly</em>?</p>

CC BY-SA 4.0
2020-01-31T20:07:01.830Z
I can’t do a de-novo assembly as I have really low coverage. I’m working with really short reads from ancient DNA.
LauraR
CC BY-SA 4.0
2020-01-31T20:36:23.543Z
You can make a consensus fasta from your bam, but that uses the reference.  If your read depth is that low, can you be confident of discrepancies between your reads and the reference?
swbarnes2
CC BY-SA 4.0
2020-02-01T08:22:25.523Z
No, but that’s another problem. Thanks. The thing I want to do can also be useful when you want to reassemble only some targeted genes from medium coverage data I guess.
LauraR
CC BY-SA 4.0
2020-01-31T18:44:38.540Z
2020-01-31T18:44:38.540Z
swbarnes2
<p>I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide. 
Is there a way to do this ? Thanks</p>

<phylogenetics><sam>
Make a fasta out of mapped reads without taking into account the reference
# Make a fasta out of mapped reads without taking into account the reference

(Asked by: username_0 on 2020-01-31)

I've aligned reads onto an haploid reference genome. I'd like to have a consensus sequence of all my aligned reads that doesn't take into account the reference genome I used (i.e. I just want the nucleotides from my reads at each base (taking into account frequency and quality) but if there is no read a one position I rather have N/? than the reference nucleotide.
Is there a way to do this ? Thanks

**Tags:** phylogenetics, sam

---

## 3 answers

### Answer by username_1 on 2020-02-13. Score: 2

vSNP outputs a VCF file that includes positions with no coverage, represented as "N".  Then there is a tool to merge a VCF file into the FASTA reference  vsnp_merge_vcf_into_fasta.py.  It will merge those "N" regions into the FASTA reference. There are options for considering frequency and quality.
 It's limited to just SNPs and those regions along the reference with no coverage.  It does not merge small indels or unmapped reads.  Pilon is a better tool for applying indels, but I'm not aware of it letting regions with no coverage be represented as "N".
vSNP_step1.py -r1 *_R1*fastq.gz -r2 *_R2*fastq.gz -r ref.fasta

Use *zc.vcf output to merge into FASTA
vsnp_merge_vcf_into_fasta.py -f ref.fasta -v *_zc.vcf


---

### Answer by username_2 on 2020-02-05. Score: 1

This is similar to this question, which links to this solution.
The difference here is that you want to go through and replace low coverage (maybe zero-coverage?) positions with Ns rather than using the reference. You can likely accomplish this using samtools depth to measure coverage and then writing a script that will go through each position in your reference, check its coverage, and string-replace with an N. Probably a bash one-liner could do this if someone very devoted to bash one-liners came up with one (I don't care for them, but other people are more creative than I). Very quick and dirty pseudocode might look like this:
### Someone could come up with something better if they spent >3min on it.
### should be easy to get it working if you use biopython
# read_and_preprocess() stands for whatever you need to do to get data in the right shape
coverage_data = read_and_preprocess('samtools_coverage.txt')
# reference_with_variants.fasta comes from the solution linked above
reference_fasta = read_and_preprocess('reference_with_variants.fasta')

threshold = 0
fasta_string = ''
for position in reference_fasta:
  if coverage_data[position] <= threshold:
    fasta_string += "N"
  else:
    fasta_string += reference_fasta[position]

# nonexistent fn that writes a fasta. 
write_as_fasta(fasta_string)

I am a little skeptical of this approach because it seems like if you do not trust the reference sequence to stand in for zero-coverage areas, then you probably shouldn't trust it to align to and call variants against. Read mapping can be a rather noisy process, especially when you are using a more distantly related reference. I think this is related to the @swbarnes2 comment.
I would suggest instead indicating the regions of low/zero coverage by putting them in lowercase acgt with your covered regions uppercase ACGT, which preserves the coverage information that you care about but still allows you access to the reference information. It's not a perfect solution, but I think it makes the most out of the marginal data you're working with, and it is trivial to change the above example code to do it.
Hope that helps-


---

### Answer by username_3 on 2020-01-31. Score: 0

Have you googled *de novo assembly*?

> **Comments:**
> **username_0**: I can’t do a de-novo assembly as I have really low coverage. I’m working with really short reads from ancient DNA.  
> **username_3**: You can make a consensus fasta from your bam, but that uses the reference.  If your read depth is that low, can you be confident of discrepancies between your reads and the reference?  
> **username_0**: No, but that’s another problem. Thanks. The thing I want to do can also be useful when you want to reassemble only some targeted genes from medium coverage data I guess.  

---

