From nucleotide or proteine sequences to EC number using biopython

Question

[friguiahlem8]

Hi,

if I have a fasta file containing nucleotide sequences or proteines sequences is it possible to get EC number using biopython for example

1.1.1.169
1.1.1.205
1.1.1.25
1.1.1.302
1.1.1.330
1.1.1.34

ps : I'm working on fungus so I need the fungal database from NCBI but I didn't know how to downloaded or construct it from fasta file using makedb command line and the second point is since I have a large fasta query file to pass I believe that I must work woth blast+ locally but also I don't know how to proceed Thanks

[jv]

Can you say a bit more about your nucleotide and protein sequences? Are your nucleotide and protein sequences complete coding sequences of well-defined reference sequences or de novo assembled sequences from experimental data like metatranscriptomics for example? I ask b/c this can inform how to proceed.

[friguiahlem8]

jv Thank you for responding !! my nucleotide/proteine sequences are Illumina raw reads for RNA-Seq they have been deposited in the NCBI Sequence Read Archive . Genome assemblies and annotations have been deposited in the NCBI BioProject database.

I hope you can help me with the appropriate script because I have been trying but no solution !! what I did for instance is : 1/ download executable blast :https://blast.ncbi.nlm.nih.gov/doc/blast-help/downloadblastdata.html 2/ decompress : tar -zxvf ncbi-blast-2.15.0+-x64-linux.tar.gz

3/dowload database :swissprot database fasta file : https://www.uniprot.org/help/downloads 4/decompress : gunzip uniprot_sprot.fasta.gz 5/create the database : /root/blastEC/ncbi-blast-2.15.0+/bin/makeblastdb -in /root/blastEC/database/second/uniprot_sprot.fasta -dbtype prot -parse_seqids

then I runned this script only for blast for instance but I believe it's not correct at all because I don't generate the same thing as when I runned on the web ( I have the same description for all my hits , I don't have the scientific name etc ..) : I also didn't know how to retriefe the EC number :(

from Bio.Blast import NCBIXML import subprocess from io import StringIO

class YourBlastClass: OUT_FMT = "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore"

def __init__(self):
    pass

def blast_sequence_local(self, sequence, database_path, evalue, thread, cov, num_alignments):
    cmd = [
        "/root/blastEC/ncbi-blast-2.15.0+/bin/blastp",
        "-db", database_path,
        "-query", "-",
        "-evalue", str(evalue),
        "-num_threads", str(thread),
        "-qcov_hsp_perc", str(cov),
        "-num_alignments", str(num_alignments),
        "-outfmt", self.OUT_FMT,
        "-show_gis",
        "-task", "blastp-fast",
    ]

    try:
        result = subprocess.run(cmd, input=sequence, text=True, capture_output=True, check=True)
        lines = result.stdout.split('\n')
        best_hits = []

        for line in lines:
            if line:
                fields = line.split('\t')
                qseqid, sseqid, pident, length, mismatch, gapopen, qstart, qend, sstart, send, evalue, bitscore = fields

                # Add more fields as needed
                # For example, you might want to extract additional information from the descriptions

                best_hits.append({
                    "qseqid": qseqid,
                    "sseqid": sseqid,
                    "pident": pident,
                    "length": length,
                    "mismatch": mismatch,
                    "gapopen": gapopen,
                    "qstart": qstart,
                    "qend": qend,
                    "sstart": sstart,
                    "send": send,
                    "evalue": evalue,
                    "bitscore": bitscore,
                    # Add more fields as needed
                    # ...
                })

        if not best_hits:
            print("No significant hits found.")

    except subprocess.CalledProcessError as e:
        print(f"Error in BLAST search. Non-zero exit status: {e.returncode}")
        print("BLAST command:")
        print(" ".join(cmd))
        print("BLAST output:")
        print(e.stdout)
        print("BLAST error:")
        print(e.stderr)
        raise
    except Exception as e:
        print(f"Error in BLAST search: {e}")
        if 'result' in locals():
            print("BLAST output:")
            print(result.stdout)
        raise

    return best_hits

if __name__ == "__main__": blast_obj = YourBlastClass()

fasta_file_path = '/root/blastEC/P_murina_1.fasta'
database_path = '/root/blastEC/database/second/uniprot_sprot.fasta'
evalue_threshold = 1e-4
num_threads = 4
query_coverage_percentage = 60
num_alignments_to_show = 10

with open(fasta_file_path, 'r') as fasta_file:
    sequence = fasta_file.read()

best_hits = blast_obj.blast_sequence_local(sequence, database_path,
                                           evalue=evalue_threshold,
                                           thread=num_threads,
                                           cov=query_coverage_percentage,
                                           num_alignments=num_alignments_to_show)

output_file_path = '/root/blastEC/best_hits_output.txt'
with open(output_file_path, 'w') as output_file:
    if best_hits:
        output_file.write("qseqid\tsseqid\tpident\tlength\tmismatch\tgapopen\tqstart\tqend\tsstart\tsend\tevalue\tbitscore\n")
        for hit in best_hits:
            output_file.write(f"{hit['qseqid']}\t{hit['sseqid']}\t{hit['pident']}\t{hit['length']}\t{hit['mismatch']}\t{hit['gapopen']}\t{hit['qstart']}\t{hit['qend']}\t{hit['sstart']}\t{hit['send']}\t{hit['evalue']}\t{hit['bitscore']}\n")
            # Add more fields as needed
            # ...
    else:
        output_file.write("No significant hits found.\n")

print(f"Best hits have been written to {output_file_path}")

I also tried diamond https://github.com/bbuchfink/diamond/wiki and i generated my output file matches.tsv but I still have the same problem I don't know how to generate the EC number