Annotating/Finding "NNNNNN" scaffold gaps in genome assembly

Question

[Adrian Pelin]

Hello,

After denovo assembly, most algorithms (I am using Allpaths in this case) do scaffolding. This places contigs into scaffolds bridged by NNNNNN sequences.

I wanted to visually validate these bridges of NNNN sequences by mapping mate pairs back to the assemblies and checking the links. This is difficult to do because I have to manually search for these bridges.

Is there any tool, that provided a fasta file can annotate all "NNNNNN" regions into a bed file let's say? I am wondering if this can even be coded with perl simply.

Thanks,

Adrian

[Prakki Rama]

perl one-liner to do this.

[michael]

Hi there and thank you for the script. I have trouble getting it running. It says "

... , line 17 printrecord.id, match.start(), match.end(), sep='\t') ^ SyntaxError: invalid syntax

I did run it with Python27 and Biopython on a Win10 computer. Any suggestions?

[James Ashmore]

The script was written for Python 3. Either install Python 3 or try changing line 17 to :

print record.id, match.start(), match.end()

[michael]

Yes, sorry I did overlook that. However, I made it work with Python3 and changed the output to gff3 . In case someone needs this:

 #!/usr/bin/env python3

 # Import necessary packages
 import argparse
 import re
 from Bio import SeqIO

 # Parse command-line arguments
 parser = argparse.ArgumentParser()
 parser.add_argument("fasta")
 args = parser.parse_args()

 # Open FASTA, search for masked regions, print in GFF3 format
 with open(args.fasta) as handle:
     for record in SeqIO.parse(handle, "fasta"):
         for match in re.finditer('N+', str(record.seq)):
             print record.id, ".", "assembly gap", match.start(), match.end(), ".", ".", ".", ".", sep='\t')

 #use the following at CMD: FILENAME.py FILENAME.fasta >> FILENAME.gff3  here

Answer 1

[James Ashmore]

If you have BioPython installed you could use this small Python script to print such regions in BED3 format.

#!/usr/bin/env python3

# Import necessary packages
import argparse
import re
from Bio import SeqIO

# Parse command-line arguments
parser = argparse.ArgumentParser()
parser.add_argument("fasta")
args = parser.parse_args()

# Open FASTA, search for masked regions, print in BED3 format
with open(args.fasta) as handle:
    for record in SeqIO.parse(handle, "fasta"):
        for match in re.finditer('N+', str(record.seq)):
            print(record.id, match.start(), match.end(), sep='\t')

[Adrian Pelin]

Excellent, this works really well! Since I am not a python expert, what do I need to replace "with open('mm10.fa') as handle:" so that I can simply issue python FindN.py FASTA-file.fasta? Thanks again! I can also use awk '($3-$2) >= 50' x.bed to filter based on how long I want intervals to be.

[James Ashmore]

I've edited the code in my answer so that you can specify the FASTA file on the command-line.

[Adrian Pelin]

Excellent, works! Thanks!