This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Automatically annotating a feature using a genbank file

I’m wondering if there’s a tool out there that I can use for this use case.

I have a linear DNA cassette that I transformed into my organism. This linear cassette is fully sequenced and annotated and has a GenBank file.

The transformed organism will often integrate this cassette with concatamers and partial containers. For example, if I transformed AATTGG into the genome, the integration might be AATTGGAATTGGAAT.

I have a contig assembled that has this integration in it. I’m wondering there’s a tool that can take the FASTA file of the contig, the GENBANK file of the initial cassette, and return a GENBANK file of the contig.

For example. In the initial construct AATTGG has the annotation “promoter” for AAT and “cds” for TGG, I would like the resulting genbank file for AATTGGAATTGGAA to have an annotation something like “promoter.1” “cds.1” “promoter.2” “cds.2” “promoter.3 - partial” or something like that.

annotation genbank wgs
• 10,628 views • link
updated by Kevin Blighe Kevin Blighe • written by Mark • Mark

Don't forget to follow up on your threads. If an answer was helpful, you should upvote it; if the answer resolved your question, you should mark it as accepted. You can accept more than one answer if they all work. If an answer was not really helpful or did not work, provide detailed feedback so others know not to use that answer.

Upvote|Bookmark|Accept

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

1 answer

To annotate your contig with repeated/partial features from the cassette GenBank, align the cassette sequence to the contig using minimap2 (v2.28, current as of 2025). This handles concatamers well.

Extract the cassette sequence from GenBank using BioPython:

from Bio import SeqIO
cassette = SeqIO.read("cassette.gb", "genbank")
SeqIO.write(cassette, "cassette.fa", "fasta")

Align with minimap2:

minimap2 -a -x map-ont contig.fa cassette.fa > alignments.sam  # Adjust preset if not ONT

Convert to BAM and sort:

samtools view -b alignments.sam | samtools sort -o alignments.bam
samtools index alignments.bam

Query alignments with samtools or bedtools to get positions. For each hit, transfer features from the original GenBank, appending suffixes like ".1", ".2", and "-partial" if clipped (check CIGAR for soft-clips > threshold, e.g., 10bp).

Use BioPython to build the new GenBank:

from Bio import SeqIO
from Bio.SeqFeature import SeqFeature, FeatureLocation

contig = SeqIO.read("contig.fa", "fasta")
# Parse alignments, add features with modified qualifiers
# e.g., feature.qualifiers['note'] = "promoter.1"
SeqIO.write(contig, "annotated_contig.gb", "genbank")

This is customizable; handle partials by prorating feature lengths. Test on small data. For automation, wrap in Snakemake.

If alignments are complex, try GMAP (v2024) instead of minimap2 for spliced-like mapping, though not typical for DNA cassettes.

Kevin

• 0 views • link
written by Kevin Blighe Kevin Blighe

Log in to answer this question.

More posts like this