VirIdAl/database_search.py at main · venus29/VirIdAl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import sys
import subprocess
from collections import namedtuple
import logging

SearchFile = namedtuple("SearchFile", "database report fields fasta evalue")

class DatabaseSearch:
    """blastn/megablast/diamond search class"""
    def __init__(self,
        search_file_nt,
        search_file_nr,
        blastn_mode,
        diamond_mode,
        threads):
        self.search_file_nt = search_file_nt
        self.search_file_nr = search_file_nr
        self.blastn_mode = blastn_mode
        self.diamond_mode = diamond_mode
        self.threads=threads

    def search(self):
        """blastn/megablast/diamond search with search_file parameters"""
        processes = []

        # blastn or megablast search
        logging.info(f"blastn search, database:{self.search_file_nt.database} file:{self.search_file_nt.fasta}")

        blastn_process = subprocess.Popen(f"blastn -db {self.search_file_nt.database} " +
            f" -outfmt '6 delim=\t {self.search_file_nt.fields}' "+
            f" -num_threads {self.threads} "+
            f" -out {self.search_file_nt.report} "+
            f" -query {self.search_file_nt.fasta} "+
            f" -evalue {self.search_file_nt.evalue} "+
            f" -task {self.blastn_mode}", shell=True)

        processes.append(blastn_process)

        # diamond blastx search
        logging.info(f"diamond blastx search, database: {self.search_file_nr.database} file: {self.search_file_nr.fasta}")

        if self.diamond_mode == "sensitive":
            diamond_process = subprocess.Popen(f"diamond blastx --sensitive --db {self.search_file_nr.database}" +
                f" --outfmt 6 {self.search_file_nr.fields}" +
                f" --threads {self.threads}" +
                f" --out {self.search_file_nr.report}" +
                f" --query {self.search_file_nr.fasta}" +
                f" --evalue {self.search_file_nr.evalue}", shell=True)
        else:
            diamond_process = subprocess.Popen(f"diamond blastx --db {self.search_file_nr.database}" +
                f" --outfmt 6 {self.search_file_nr.fields}" +
                f" --threads {self.threads}" +
                f" --out {self.search_file_nr.report}" +
                f" --query {self.search_file_nr.fasta}" +
                f" --evalue {self.search_file_nr.evalue}", shell=True)

        processes.append(diamond_process)

        for process in processes:
            returncode = process.wait()
            if returncode!=0:
                logging.error(f"search failed {process.args}")
                sys.exit()