Answers.Body: <p>ABRicate can handle multiple input files<sup><a href="https://github.com/tseemann/abricate#input" rel="nofollow noreferrer">1</a></sup> and will set the <em>FILE</em> column in it's output accordingly. Assuming all of the FASTA files in this column exist on disk, you can try the following <a href="/questions/tagged/python" class="post-tag" title="show questions tagged &#39;python&#39;" aria-label="show questions tagged &#39;python&#39;" rel="tag" aria-labelledby="tag-python-tooltip-container">python</a> solution. It uses <a href="https://pypi.org/project/pyfaidx/" rel="nofollow noreferrer">pyfaidx</a>, 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 <code>tsv2gff3.py</code>:</p>
  <pre><code>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
  </code></pre>
  <pre><code>def get_argument_parser():
      parser = argparse.ArgumentParser(add_help=False)
      group = parser.add_argument_group('positional arguments')
      group.add_argument(
          &quot;tsv&quot;,
          type=str,
          metavar=&quot;FILE&quot;,
          help=&quot;The ABRicate TSV file&quot;,
      )
      group = parser.add_argument_group('output options')
      group.add_argument(
          &quot;-o&quot;,
          &quot;--output&quot;,
          type=str,
          metavar=&quot;FILE&quot;,
          default='-',
          help=&quot;Write the GFF3 output to FILE (default: stdout)&quot;,
      )
      group.add_argument(
          &quot;-f&quot;,
          &quot;--force&quot;,
          action='store_true',
          help=&quot;Overwrite the output file if it exists&quot;,
      )
      group = parser.add_argument_group('additional options')
      group.add_argument(
          &quot;-h&quot;,
          &quot;--help&quot;,
          action=&quot;help&quot;,
          help=&quot;Show this help message and exit&quot;,
      )
  
      return parser
  </code></pre>
  <pre><code>@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'
  </code></pre>
  <pre><code>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()
  </code></pre>
  <pre><code>        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, &quot;CDS&quot;, start, end, '.', '+', '0', *attrs)
  
              out.write(str(feature))
  </code></pre>
  <pre><code>        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'&gt;{seq_id}\n{seq}\n')
  </code></pre>
  <pre><code>if __name__ == '__main__':
      main()
  </code></pre>
  <p>And <code>pyfaidx</code> installed, the usage is:</p>
  <pre><code>$ 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
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-11T16:09:10.597Z
Answers.Text: Solution provided by @Steve is elaborate, though import of Literal from typing was causing problem.
Answers.UserDisplayName: com
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-04-12T23:37:01.717Z
Answers.Text: @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.
Answers.UserDisplayName: Steve
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-29T06:28:15.943Z
Answers.LastActivityDate: 2023-03-29T06:28:15.943Z
Answers.UserDisplayName: Steve
Body: <p>Since such a neat solution <a href="https://bioinformatics.stackexchange.com/questions/20749/convert-abricate-output-tsv-file-to-gff3-format">(abricate tsv to gff3)</a> was provided by <a href="https://bioinformatics.stackexchange.com/users/10775/steve">Steve</a>, 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.</p>
  <p>I have two files - (1) fasta file with .fna extension, and (2) a data file with .tsv extension generated from <a href="https://github.com/tseemann/abricate" rel="nofollow noreferrer">abricate</a>.</p>
  <p>I need to write a gff3 file (containing sequences) by combining these two files (.fna and .tsv)</p>
  <p><em>A .gff3 file generally has three or four components as below.</em> Details <a href="https://asia.ensembl.org/info/website/upload/gff3.html" rel="nofollow noreferrer">here</a>.</p>
  <ol>
  <li>Declaration of file type and version</li>
  <li>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).</li>
  <li>Sequence features</li>
  <li>Fasta sequence (optional)</li>
  </ol>
  <p>See below a truncated example</p>
  <p>##gff-version 3
  <br>##sequence-region NZ_BHEI01000001 1 13890
  <br>NZ_BHEI01000001 Prodigal    CDS 29  1588    .   +   0   ID=CEMHHHDG_00001;product=hypothetical protein
  <br>NZ_BHEI01000001 Prodigal    CDS 1585    2289    .   +   0   ID=CEMHHHDG_00002;product=Putative beta-lactamase HcpC
  <br>....
  <br>##FASTA
  <br>&gt;NZ_BHEI01000001
  <br>AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
  <br>TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
  <br>TATAGAGGAGTTGGCAAAATACACAGATCAAA..........</p>
  <p>Notice that the sequence names has to match in every instance (declaration of sequence regions, sequence features, and the fasta header at the end)</p>
  <p>The .tsv to .gff3 conversion was elegantly solved by <a href="https://bioinformatics.stackexchange.com/users/10775/steve">Steve</a>.</p>
  <p>Now we need to add the declaration of sequence region(s) in the .gff3 file. The declaration comprises of
  <br>&quot;##sequence-region&quot;[space]the input file name without the file extension .fna[space]sequence start[space]sequence length[\n]</p>
  <p>For the last part, I have used the following awk code on the .fna file to generate the length in desired format.</p>
  <pre><code>awk '/^&gt;/ {if (seqlen); seqlen=0;next; } { seqlen += length($0)}END{print &quot;1 &quot; seqlen}' example.fna
  </code></pre>
  <p>This output needs to be appended in the declaration of sequence region line.</p>
  <p>Now we need to append the fasta sequence in the following format</p>
  <p><br>##FASTA
  <br>&gt;Sequence name as in the parts above in the .gff3 file
  <br>AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
  <br>TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
  <br>TATAGAGGAGTTGGCAAAATACACAGATCAAA..........</p>
  <p><strong>Is there any help in sight, please?</strong></p>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-22T14:17:09.900Z
Comments.Text: Please use the [formatting tools](https://bioinformatics.stackexchange.com/help/formatting) 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T15:52:45.390Z
Comments.Text: 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.
Comments.UserDisplayName: com
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T15:54:45.287Z
Comments.Text: 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 :)
Comments.UserDisplayName: terdon
Tags: <gff3><tsv>
Title: Improving conversion of abricate tsv file to gff3 file
ThreadText: # Improving conversion of abricate tsv file to gff3 file
  
  (Asked by: username_0 on 2023-03-22)
  
  Since such a neat solution [(abricate tsv to gff3)](https://bioinformatics.stackexchange.com/questions/20749/convert-abricate-output-tsv-file-to-gff3-format) was provided by [username_2](https://bioinformatics.stackexchange.com/users/10775/username_2), 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](https://github.com/tseemann/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](https://asia.ensembl.org/info/website/upload/gff3.html).
  
  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 [username_2](https://bioinformatics.stackexchange.com/users/10775/username_2).
  
  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?**
  
  **Tags:** gff3, tsv
  
  > **Comments:**
  > **username_1**: Please use the [formatting tools](https://bioinformatics.stackexchange.com/help/formatting) 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.  
  > **username_0**: Thank you @username_1. I do understand that my formatting was not correct. I was also feeling so. Trying to improve! Please bear with for the time being.  
  > **username_1**: 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 :)  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2023-03-29. Score: 1
  
  ABRicate can handle multiple input files[1](https://github.com/tseemann/abricate#input) 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 [python](/questions/tagged/python "show questions tagged 'python'") solution. It uses [pyfaidx](https://pypi.org/project/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
  
  ```
  
  > **Comments:**
  > **username_0**: Solution provided by @username_2 is elaborate, though import of Literal from typing was causing problem.  
  > **username_2**: @username_0 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.  
  
  ---
  
  