2
$\begingroup$

Since such a neat solution (abricate tsv to gff3) was provided by Steve, here are few other steps that I am looking to add so that the script progress to logical maturity to be usable by many others.

I have two files - (1) fasta file with .fna extension, and (2) a data file with .tsv extension generated from abricate.

I need to write a gff3 file (containing sequences) by combining these two files (.fna and .tsv)

A .gff3 file generally has three or four components as below. Details here.

  1. Declaration of file type and version
  2. Declaration of sequence region(s) starting from Name of the sequence, followed by single space, then 1 separated by single space, ending in numerical value of the sequence length of the fasta sequence (which is the .fna file).
  3. Sequence features
  4. Fasta sequence (optional)

See below a truncated example

##gff-version 3
##sequence-region NZ_BHEI01000001 1 13890
NZ_BHEI01000001 Prodigal CDS 29 1588 . + 0 ID=CEMHHHDG_00001;product=hypothetical protein
NZ_BHEI01000001 Prodigal CDS 1585 2289 . + 0 ID=CEMHHHDG_00002;product=Putative beta-lactamase HcpC
....
##FASTA
>NZ_BHEI01000001
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........

Notice that the sequence names has to match in every instance (declaration of sequence regions, sequence features, and the fasta header at the end)

The .tsv to .gff3 conversion was elegantly solved by Steve.

Now we need to add the declaration of sequence region(s) in the .gff3 file. The declaration comprises of
"##sequence-region"[space]the input file name without the file extension .fna[space]sequence start[space]sequence length[\n]

For the last part, I have used the following awk code on the .fna file to generate the length in desired format.

awk '/^>/ {if (seqlen); seqlen=0;next; } { seqlen += length($0)}END{print "1 " seqlen}' example.fna

This output needs to be appended in the declaration of sequence region line.

Now we need to append the fasta sequence in the following format


##FASTA
>Sequence name as in the parts above in the .gff3 file
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........

Is there any help in sight, please?

$\endgroup$
3
  • 1
    $\begingroup$ Please use the formatting tools as I requested in your previous question. They really help clarify and that way we can know that the data are correct. It is also much easier than laboriously adding <br> tags. $\endgroup$
    – terdon
    Commented Mar 22, 2023 at 14:17
  • $\begingroup$ Thank you @terdon. I do understand that my formatting was not correct. I was also feeling so. Trying to improve! Please bear with for the time being. $\endgroup$
    – com
    Commented Mar 23, 2023 at 15:52
  • $\begingroup$ Oh hey, that wasn't a rebuke! You spent a lot of effort trying to get your post to look good, it's just that the effort was misplaced because we don't make it very easy to find the right help documentation. I was trying to avoid having you put in so much effort, and ensure that the question is clear and easy to read. So, thank you for the clear effort you put in, even if it did happen to be more than was actually necessary :) $\endgroup$
    – terdon
    Commented Mar 23, 2023 at 15:54

1 Answer 1

1
$\begingroup$

ABRicate can handle multiple input files1 and will set the FILE column in it's output accordingly. Assuming all of the FASTA files in this column exist on disk, you can try the following solution. It uses pyfaidx, which allows for fast random access to the sequences in the (indexed) FASTA file(s). If the FASTA indexes don't exist, they will be created. It loads only a minimal amount of the input FASTA file(s) into memory. With the following in a file called tsv2gff3.py:

import argparse
import csv
import contextlib
import re
import sys
import textwrap

from dataclasses import dataclass, astuple
from pathlib import Path
from typing import List, Literal, Tuple

from pyfaidx import Fasta
def get_argument_parser():
    parser = argparse.ArgumentParser(add_help=False)
    group = parser.add_argument_group('positional arguments')
    group.add_argument(
        "tsv",
        type=str,
        metavar="FILE",
        help="The ABRicate TSV file",
    )
    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        metavar="FILE",
        default='-',
        help="Write the GFF3 output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )
    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser
@dataclass
class Feature:

    seq_id: str
    source: str
    type: str
    start: int
    end: int
    score: float
    strand: Literal['-', '+']
    phase: Literal['0', '1', '2']
    attrs: List[Tuple[str, str]]

    @property
    def attributes(self):
        return ';'.join([f'{key}={val}' for key, val in self.attrs])

    def __iter__(self):
        return iter(astuple(self))

    def __str__(self):
        vals = list(map(str, astuple(self)[:-1])) + [self.attributes]
        return '\t'.join(vals) + '\n'
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.output != '-':
            out = stack.enter_context(open(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        tsv = stack.enter_context(open(args.tsv, 'r'))
        first_line = tsv.readline().rstrip()
        fieldnames = first_line.lstrip('#').split('\t')
        reader = csv.DictReader(tsv, dialect='excel-tab', fieldnames=fieldnames)

        fasta_files = dict()
        seq_ids = dict()
        out.write('##gff-version 3\n')
        for idx, row in enumerate(reader, 1):

            filename = Path(row['FILE'])
            if filename not in fasta_files:
                fasta_files[str(filename)] = stack.enter_context(Fasta(filename))

            if filename.suffix in ('.fasta', '.fa', '.fna'):
                seq_id = filename.stem
            else:
                seq_id = filename.name

            fasta = fasta_files[str(filename)]
            seq = fasta[row['SEQUENCE']]

            if not seq_id in seq_ids:
                out.write(f'##sequence-region {seq_id} 1 {len(seq)}\n')
                seq_ids[seq_id] = [row['FILE'], row['SEQUENCE']]

            database = row['DATABASE']
            product = re.sub(r'[\[\]\(\)]', '', row['PRODUCT'])
            start, end = int(row['START']), int(row['END'])
            attrs = [
                ('ID', f'{seq_id}_{idx}'),
                ('Product', product),
            ],

            feature = Feature(seq_id, database, "CDS", start, end, '.', '+', '0', *attrs)

            out.write(str(feature))
        out.write('##FASTA\n')
        for seq_id, (filename, seq_name) in seq_ids.items():
            fasta_rec = fasta_files[filename][seq_name]
            seq = textwrap.fill(str(fasta_rec), width=60)

            out.write(f'>{seq_id}\n{seq}\n')
if __name__ == '__main__':
    main()

And pyfaidx installed, the usage is:

$ python tsv2gff3.py -h
usage: tsv2gff3.py [-o FILE] [-f] [-h] FILE

positional arguments:
  FILE                  The ABRicate TSV file

output options:
  -o FILE, --output FILE
                        Write the GFF3 output to FILE (default: stdout)
  -f, --force           Overwrite the output file if it exists

additional options:
  -h, --help            Show this help message and exit
$\endgroup$
2
  • $\begingroup$ Solution provided by @Steve is elaborate, though import of Literal from typing was causing problem. $\endgroup$
    – com
    Commented Apr 11, 2023 at 16:09
  • $\begingroup$ @com You might not need pyfaidx if memory is not a concern. It just wasn't clear to me how much memory would be required to load all of the sequences into memory. You can remove the type hints if they're causing problems. $\endgroup$
    – Steve
    Commented Apr 12, 2023 at 23:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.