I uploaded a single FASTA file with multiple gene clusters from different organisms to an online program called antiSMASH, or fungiSMASH in my case. Most clusters have a ketosynthase (KS) gene in them.  antiSMASH identified the putative KS genes and provided me with an output in a text file. Can anyone assist me in extracting (parsing may be the correct term?) the genes of interest (nucleotide sequences) with associated accession number and definition, or just the taxon name from the definition?

I believe all of the genes of interest will say:

    /aSDomain="PKS_KS"

And if this is present then I will want the range of nucleotides indicated adjacent to the heading aSDomain. For example:

    aSDomain        2527..3816

However, sometimes the adjacent numbers will says something like:

    aSDomain        join(1610..1702,1756..2109,2165..3010)

In which case I believe I would want to concatenate each range indicated.

Or:

    aSDomain        complement(join(10640..10648,10717..11439))

In which case I believe I would want to concatenate all ranges and then take the complementary sequence.



I would like to do an alignment and then make a phylogenetic tree based on the extracted KS genes. I believe FASTA format would be a good output to have my KS genes in, but I can convert if necessary.

This is my first time using antiSMASH and I'm new to coding so I apologies for any obvious blunders and I would have preferred to attach a file of my output data but I didn't see that as an option! Thanks in advance for any help!

Here's a link to my output:

https://drive.google.com/open?id=1KWbh3D7jY7u5AytlGCLxC65MR_KKUY7X

If someone has a better way of attaching a large text file (~1.7 million characters), I'm all ears.

Here's a link all of the output that antiSMASH generated (not just text file of all annotated genes):

https://drive.google.com/open?id=1KWbh3D7jY7u5AytlGCLxC65MR_KKUY7X



Hi [mac03pat][1],

This should work:

```
#!/usr/bin/env python

from Bio import SeqIO
from Bio.SeqFeature import SeqFeature

gbk = "antiSMASH_processed_geneclusters.txt"
fa = "antiSMASH_processed_geneclusters.fa"
input_handle = open(gbk, "r")
output_handle = open(fa, "w")

for record in SeqIO.parse(input_handle, "genbank"):
    features = [feature for feature in record.features if feature.type == "aSDomain"]
    for feature in features:
        if feature.qualifiers["aSDomain"][0] == "PKS_KS":
            output_handle.write(">%s %s\n%s\n" % (
                record.id,
                feature.location,
                SeqFeature(feature.location).extract(record.seq)))

output_handle.close()
input_handle.close()
```

  [1]: https://www.biostars.org/u/55218/