text: user6459: 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
  
  user2108: 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-
  
  user6343: 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
  
  
original_text: 6459: <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>
   2108: <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>
   6343: <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>
  