Answers.Body: <p>I took a quick look at the <a href="https://github.com/ablab/quast/blob/master/quast_libs/contigs_analyzer.py" rel="nofollow noreferrer">quast code</a> that I believe is responsible here:</p>
  <pre><code>def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
      indels_info = IndelsInfo()
      genome_mapping = {}
      for chr_name, chr_len in reference_chromosomes.items():
          genome_mapping[chr_name] = [0] * (chr_len + 1)
      with open(used_snps_fpath, 'w') as used_snps_f:
          for chr_name, aligns in ref_aligns.items():
              for align in aligns:
                  ref_pos, ctg_pos = align.s1, align.s2
                  strand_direction = 1 if align.s2 &lt; align.e2 else -1
                  for op in parse_cs_tag(align.cigar):
                      if op.startswith(':'):
                          n_bases = int(op[1:])
                      else:
                          n_bases = len(op) - 1
                      if op.startswith('*'):
                          ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                          if ctg_nucl != 'N' and ref_nucl != 'N':
                              # MISMATCHES INCREMENTED HERE
                              indels_info.mismatches += 1
                              if qconfig.show_snps:
                                  used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ref_pos += 1
                          ctg_pos += 1 * strand_direction
                      elif op.startswith('+'):
                          # INDELS INCREMENTED HERE
                          indels_info.indels_list.append(n_bases)
                          indels_info.insertions += n_bases
                          if qconfig.show_snps and n_bases &lt; qconfig.MAX_INDEL_LENGTH:
                              ref_nucl, ctg_nucl = '.', op[1:].upper()
                              used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ctg_pos += n_bases * strand_direction
                      elif op.startswith('-'):
                          indels_info.indels_list.append(n_bases)
                          indels_info.deletions += n_bases
                          if qconfig.show_snps and n_bases &lt; qconfig.MAX_INDEL_LENGTH:
                              ref_nucl, ctg_nucl = op[1:].upper(), '.'
                              used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ref_pos += n_bases
                      else:
                          ref_pos += n_bases
                          ctg_pos += n_bases * strand_direction
                  if align.s1 &lt; align.e1:
                      for pos in range(align.s1, align.e1 + 1):
                          genome_mapping[align.ref][pos] = 1
                  else:
                      for pos in range(align.s1, len(genome_mapping[align.ref])):
                          genome_mapping[align.ref][pos] = 1
                      for pos in range(1, align.e1 + 1):
                          genome_mapping[align.ref][pos] = 1
              for i in ns_by_chromosomes[align.ref]:
                  genome_mapping[align.ref][i] = 0
  
      covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
      return covered_ref_bases, indels_info
  </code></pre>
  <p>This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.</p>
  <p>This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).</p>
  <p>From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast <code>used_snps</code> file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-01-28T00:48:50.603Z
Answers.LastActivityDate: 2021-01-28T00:48:50.603Z
Answers.UserDisplayName: Maximilian Press
Body: <p>When analyzing Quast results it seems that it doesn't calculate mismatches and indels in a useful way if the &quot;Duplication ratio&quot; is over 1.
  For example, that's what I get for an assembly with two contigs covering almost the full genome
  <a href="https://i.sstatic.net/iGGbL.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/iGGbL.png" alt="enter image description here" /></a></p>
  <p>And that's what I get if I remove one of the contigs:
  <a href="https://i.sstatic.net/L91VH.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/L91VH.png" alt="enter image description here" /></a></p>
  <p>See how the mismatches and indels per 100kbp are almost double in the first case but that does not correspond to the reality.</p>
  <p>Any better way to get these metrics? Would it be safe to divide error rates by the duplication ratio?</p>
  
Tags: <ngs><assembly>
Title: Quast duplication ratio and mismatches percent
ThreadText: # Quast duplication ratio and mismatches percent
  
  (Asked by: username_0 on 2021-01-22)
  
  When analyzing Quast results it seems that it doesn't calculate mismatches and indels in a useful way if the "Duplication ratio" is over 1.
  For example, that's what I get for an assembly with two contigs covering almost the full genome
  [![enter image description here](https://i.sstatic.net/iGGbL.png)](https://i.sstatic.net/iGGbL.png)
  
  And that's what I get if I remove one of the contigs:
  [![enter image description here](https://i.sstatic.net/L91VH.png)](https://i.sstatic.net/L91VH.png)
  
  See how the mismatches and indels per 100kbp are almost double in the first case but that does not correspond to the reality.
  
  Any better way to get these metrics? Would it be safe to divide error rates by the duplication ratio?
  
  **Tags:** ngs, assembly
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-01-28. Score: 1 | Accepted answer
  
  I took a quick look at the [quast code](https://github.com/ablab/quast/blob/master/quast_libs/contigs_analyzer.py) that I believe is responsible here:
  
  ```
  def analyze_coverage(ref_aligns, reference_chromosomes, ns_by_chromosomes, used_snps_fpath):
      indels_info = IndelsInfo()
      genome_mapping = {}
      for chr_name, chr_len in reference_chromosomes.items():
          genome_mapping[chr_name] = [0] * (chr_len + 1)
      with open(used_snps_fpath, 'w') as used_snps_f:
          for chr_name, aligns in ref_aligns.items():
              for align in aligns:
                  ref_pos, ctg_pos = align.s1, align.s2
                  strand_direction = 1 if align.s2 < align.e2 else -1
                  for op in parse_cs_tag(align.cigar):
                      if op.startswith(':'):
                          n_bases = int(op[1:])
                      else:
                          n_bases = len(op) - 1
                      if op.startswith('*'):
                          ref_nucl, ctg_nucl = op[1].upper(), op[2].upper()
                          if ctg_nucl != 'N' and ref_nucl != 'N':
                              # MISMATCHES INCREMENTED HERE
                              indels_info.mismatches += 1
                              if qconfig.show_snps:
                                  used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ref_pos += 1
                          ctg_pos += 1 * strand_direction
                      elif op.startswith('+'):
                          # INDELS INCREMENTED HERE
                          indels_info.indels_list.append(n_bases)
                          indels_info.insertions += n_bases
                          if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                              ref_nucl, ctg_nucl = '.', op[1:].upper()
                              used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ctg_pos += n_bases * strand_direction
                      elif op.startswith('-'):
                          indels_info.indels_list.append(n_bases)
                          indels_info.deletions += n_bases
                          if qconfig.show_snps and n_bases < qconfig.MAX_INDEL_LENGTH:
                              ref_nucl, ctg_nucl = op[1:].upper(), '.'
                              used_snps_f.write('%s\t%s\t%d\t%s\t%s\t%d\n' % (chr_name, align.contig, ref_pos, ref_nucl, ctg_nucl, ctg_pos))
                          ref_pos += n_bases
                      else:
                          ref_pos += n_bases
                          ctg_pos += n_bases * strand_direction
                  if align.s1 < align.e1:
                      for pos in range(align.s1, align.e1 + 1):
                          genome_mapping[align.ref][pos] = 1
                  else:
                      for pos in range(align.s1, len(genome_mapping[align.ref])):
                          genome_mapping[align.ref][pos] = 1
                      for pos in range(1, align.e1 + 1):
                          genome_mapping[align.ref][pos] = 1
              for i in ns_by_chromosomes[align.ref]:
                  genome_mapping[align.ref][i] = 0
  
      covered_ref_bases = sum([sum(genome_mapping[chrom]) for chrom in genome_mapping])
      return covered_ref_bases, indels_info
  
  ```
  
  This indicates to me that it simply iterates over contigs (or rather, called SNPs from contig alignments to a reference) and counts up the number of mismatches/indels that occur in total across all aligned contigs.
  
  This suggests to me that if haplotigs are present, the areas that they cover will be double-counted (or n-counted where there are n haplotigs, assuming that all haplotigs covering a reference position have roughly equal divergence).
  
  From this point of view, I think that it is heuristically reasonable to divide by genome fraction. However, if you want something really rigorous you might consider digging into the quast `used_snps` file that this code iterates over to decide which SNPs/indels/whatever you think are worth keeping.
  
  
  ---
  
  