How to plot the gap distribution of contigs wrt to reference genome

Question

I have a contigs file that I generated using Minia and then I have used Minimap2 to map these contigs again to the reference genome. Now I want to plot the gap distribution i.e. gap distance of the set of contigs with respect to the positions in the reference genome. For example: something like , if there are m contigs obtained from the Minia, then X axis will show 0 to m-1 gap points (in an ideal case where there is no overlap between the contigs and the contigs are aligned back to back without space).

The file from minimap2 looks like this:

0   17679   0   17679   +   NZ_BCML01000091.1   39426   18736   36415   17679   17679   60  NM:i:0  ms:i:35358  AS:i:35358  nn:i:0  tp:A:P  cm:i:3274   s1:i:17675  s2:i:49 de:f:0  rl:i:75 cg:Z:17679M
1   16011   0   16011   +   NZ_BCML01000021.1   105137  86      16097   16011   16011   60  NM:i:0  ms:i:32022  AS:i:32022  nn:i:0  tp:A:P  cm:i:3008   s1:i:15996  s2:i:0  de:f:0  rl:i:15 cg:Z:16011M
2   2779    0   2779    -   NZ_BCML01000048.1   72397   14902   17681   2779    2779    60  NM:i:0  ms:i:5558   AS:i:5558   nn:i:0  tp:A:P  cm:i:519    s1:i:2751   s2:i:0  de:f:0  rl:i:45 cg:Z:2779M
4   177     0   177     +   NZ_BCML01000150.1   16425   15710   15887   177     177     60  NM:i:0  ms:i:354    AS:i:354    nn:i:0  tp:A:P  cm:i:30     s1:i:160    s2:i:0  de:f:0  rl:i:15 cg:Z:177M    

where the 1st column is the contig ID and 8th and 9th column show the start and end position of the contig on the reference genome.

I have used the samtools for obtaining a length distribution plot, but not understanding how to do a gap distribution plot. I can use awk on linux or python or R. Anything that is going to work. Thank you for all the help.

Answer 1

a better solution

I wasn't too happy with the previous answer that I gave, so I was fiddling around, and I realized that bedtools genomecov will give you this for free with the following workflow:

minimap2 -ax asm10 ref.fasta contigs.fasta | \
samtools sort --output_fmt BAM > contigs_to_ref.bam
# -bga means output in bedgraph format, including all positions
bedtools genomecov -ibam contigs_to_ref.bam -bga > contig_cov.bg

If you look at the resulting file, it will report all regions of identical coverage:

head contig_cov.bg
# prints the following:
# reference chr/contig, start position, end position, coverage
NODE_1_length_442924_cov_598.199    0   8151    0
NODE_1_length_442924_cov_598.199    8151    8392    1
NODE_1_length_442924_cov_598.199    8392    20203   0
NODE_1_length_442924_cov_598.199    20203   20463   1
NODE_1_length_442924_cov_598.199    20463   44953   0
NODE_1_length_442924_cov_598.199    44953   45216   1
NODE_1_length_442924_cov_598.199    45216   52633   0
NODE_1_length_442924_cov_598.199    52633   52899   1
NODE_1_length_442924_cov_598.199    52899   57698   0
NODE_1_length_442924_cov_598.199    57698   58207   1

You can simply filter this file for all regions wherein the fourth column is equal to zero, and take the first 3 columns, and that is a list of gaps in the assembly relative to the reference:

awk '$4 == "0" { print $0 }' contig_cov.bg | cut -f1-3
# prints
NODE_1_length_442924_cov_598.199    0   8151
NODE_1_length_442924_cov_598.199    8392    20203
NODE_1_length_442924_cov_598.199    20463   44953
NODE_1_length_442924_cov_598.199    45216   52633
NODE_1_length_442924_cov_598.199    52899   57698
NODE_1_length_442924_cov_598.199    58207   73398
NODE_1_length_442924_cov_598.199    73733   78601
NODE_1_length_442924_cov_598.199    78959   80147
NODE_1_length_442924_cov_598.199    80419   86967
NODE_1_length_442924_cov_598.199    87223   93872

You could even do this as a one-liner probably though that's getting a little convoluted.

Original answer

It's a bit circuitous, but one way to do this would be to follow this answer's suggestion and generate a pileup for the alignment of your contigs to the reference.

You can get minimap2 to give SAM/BAM output with the -a flag. You then use the SAM/BAM to generate a pileup or a coverage file, using samtools mpileup or samtools depth -a. You can then use a heuristic to find contiguous regions with zero coverage of contigs. You could for example loop through a depth file, here is some python code that gives the basic idea:

#!/usr/bin/env python

# minimap2 -ax asm10 ref.fasta contigs.fasta | \
# samtools sort --output_fmt BAM > contigs_to_ref.bam
# samtools depth contigs_to_ref.bam > depth.txt
depth_file = "depth.txt"

with open(depth_file) as dep_handle:
  gap_sizes = list()
  start = None
  pos = None
  for line in dep_handle:
    # needs also logic to handle multiple chromosomes in reference 
    # but I'm lazy
    fields = line.split()
    if int(fields[2]) == 0:
      pos = int(fields[2])
      if start is None:
        start = int(fields[1])
      
    elif fields[2] >0:
      if start is None:
        continue
      else:
        gap_sizes.append(pos - start)
        start = None

print(gap_sizes)
        

Probably there is a simpler way to do this but I don't know it.