#!/usr/bin/python3
import subprocess
import time
import os

query = "query.fasta"
db_fasta = "MIDORI2_LONGEST_NUC_GB265_srRNA_RAW1.fasta"
blastn = "./blastn2.10.0"

########################################
#### BLASTN
########################################
print("#### BLASTN ####")

# DBが未作成なら作成
if not os.path.exists(db_fasta + ".nhr"):
    cmd_makedb = f"./makeblastdb2.10.0 -dbtype nucl -in {db_fasta}"
    subprocess.call(cmd_makedb, shell=True)

startTime_B = time.time()
cmd_B = f"{blastn} -num_descriptions 0 -task blastn -word_size 11 -query {query} -db {db_fasta} -num_alignments 5 -evalue 1e-3 -out out_blastn.txt"
print("cmd_B:", cmd_B)
subprocess.call(cmd_B, shell=True)
elapsed_B = round(time.time() - startTime_B, 2)
print(f"#### BLASTN elapsed_time: {elapsed_B} sec\n")

########################################
#### MMseqs2
########################################
print("#### MMseqs2 ####")

# DBが未作成なら作成
if not os.path.exists("mmseqs_db"):
    os.makedirs("mmseqs_db")
if not os.path.exists("mmseqs_db/targetDB"):
    cmd_mmdb = f"mmseqs createdb {db_fasta} mmseqs_db/targetDB"
    subprocess.call(cmd_mmdb, shell=True)
if not os.path.exists("mmseqs_db/queryDB"):
    cmd_mmdb_query = f"mmseqs createdb {query} mmseqs_db/queryDB"
    subprocess.call(cmd_mmdb_query, shell=True)

startTime_M = time.time()
cmd_M = (
    f"mmseqs search mmseqs_db/queryDB mmseqs_db/targetDB mmseqs_db/resultDB mmseqs_db/tmp "
    f"--search-type 3 -e 1e-3 --max-seqs 5 --threads 4"
)
print("cmd_M:", cmd_M)
subprocess.call(cmd_M, shell=True)

cmd_convert = (
    f"mmseqs convertalis mmseqs_db/queryDB mmseqs_db/targetDB mmseqs_db/resultDB "
    f"out_mmseqs.txt --format-mode 0"
)
subprocess.call(cmd_convert, shell=True)
elapsed_M = round(time.time() - startTime_M, 2)
print(f"#### MMseqs2 elapsed_time: {elapsed_M} sec\n")

########################################
#### Summary
########################################
print("==== Summary ====")
print(f"BLASTN:  {elapsed_B} sec")
print(f"MMseqs2: {elapsed_M} sec")
if elapsed_M > 0:
    print(f"Speed ratio (BLASTN/MMseqs2): {round(elapsed_B/elapsed_M, 1)}x")
