download viral database

 cp download_db.py down_db.py.backup
 sed -i -e 's/39733/1396/g' download_db.py
 #NOTE the download_db.py cannot run alone, rather than via "vrap.py -u"
 (damian) jhuang@hamburg:~/DATA/Data_Tam_Acinetobacter_baumannii$ ./vrap/vrap.py -u
 mv vrap/database/viral_db vrap/database/viral_db_Bacillus_cereus

run vrap.py

 #--host=/home/jhuang/REFs/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa 
 vrap.py -1 ./results_ATCC19606/trimmed/A10_CraA_HQ_trimmed_P_1.fastq.gz -2 ./results_ATCC19606/trimmed/A10_CraA_HQ_trimmed_P_2.fastq.gz -o A10_vrap_out --virus=/home/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/ATCC19606.fasta -n /media/jhuang/Titisee/GAMOLA2/Blast_db/nt -a /media/jhuang/Titisee/GAMOLA2/Blast_db/nr -t 4
 #vrap/vrap.py -1 220211_M03701_0272_000000000-JWYW9/p12568/BovHepV_Bulgaria_19_S3_R1_001.fastq.gz -2 220211_M03701_0272_000000000-JWYW9/p12568/BovHepV_Bulgaria_19_S3_R2_001.fastq.gz -o BovHepV_Bulgaria_19_on_Hepacivirus --host=/home/jhuang/REFs/Bos_taurus.ARS-UCD1.2.dna.toplevel.fa  -t 3
 #IMPORTANT, it doesn't work, since the update will always downlod the repeated sequences causing ERROR! It needs to delete it as follows!!

 #In summary: --(optional) host --> bowtie2/host (bowtie_database for host-substraction) + blast/db/host (in nt_dbs, custom blastn database before downloaded_db!)
 #            --(optional) virus --> blast/db/virus (in nt_dbs, custom blastn database before downloaded_db!)
 #            --(always) defined in download_db.py (setting by modifying the file) and run default --> database/viral_db/viral_nucleotide (nt_dbs)
 #                                                                                        --> database/viral_db/viral_protein (prot_db)
 #            --(optional) nt+nr
 #Therefore, the results substracted host should not need the blast with host, however it it in nt_dbs!

 #Under (vrap) jhuang@hamm:~/DATA/Data_Tam_Acinetobacter_baumannii$ 
 #DEBUG: {path}/external_tools/Lighter-master/lighter-->/home/jhuang/miniconda3/envs/vrap/bin/lighter
 #DEBUG: {path}/external_tools/SPAdes-3.11.1-Linux/bin/spades.py-->/home/jhuang/miniconda3/envs/vrap/bin/spades.py

 #A6: --host=/home/jhuang/REFs/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa  --virus=/home/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/vrap/database/viral_db_Bacillus_cereus/nucleotide_Bacillus_cereus_0.01.fasta
 vrap/vrap.py -1 ./results_ATCC19606/trimmed/A6_WT_HQ_trimmed_P_1.fastq.gz -2 ./results_ATCC19606/trimmed/A6_WT_HQ_trimmed_P_2.fastq.gz -o A6_vrap_out  -n /mnt/h1/jhuang/blast/nt -a /mnt/h1/jhuang/blast/nr -t 55

 #A10　
 #START checking blast db: no host and with virus --> mapping to 5 dbs: blast/db/virus, viral_nucleotide, viral_protein, nt, nr
 vrap/vrap.py -1 ./results_ATCC19606/trimmed/A10_CraA_HQ_trimmed_P_1.fastq.gz -2 ./results_ATCC19606/trimmed/A10_CraA_HQ_trimmed_P_2.fastq.gz -o A10_vrap_out --virus=/home/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/vrap/database/viral_db_Bacillus_cereus/nucleotide_Bacillus_cereus_0.01.fasta   -n /mnt/h1/jhuang/blast/nt -a /mnt/h1/jhuang/blast/nr -t 55

 #A12
 vrap/vrap.py -1 ./results_AYE/trimmed/A12_AYE_HQ_trimmed_P_1.fastq.gz -2 ./results_AYE/trimmed/A12_AYE_HQ_trimmed_P_2.fastq.gz -o A12_vrap_out  -n /mnt/h1/jhuang/blast/nt -a /mnt/h1/jhuang/blast/nr -t 55

 #A19 
 vrap/vrap.py -1 ./results_ATCC17978/trimmed/A19_17978_HQ_trimmed_P_1.fastq.gz -2 ./results_ATCC17978/trimmed/A19_17978_HQ_trimmed_P_2.fastq.gz -o A19_vrap_out  -n /mnt/h1/jhuang/blast/nt -a /mnt/h1/jhuang/blast/nr -t 55

 #START VrAP
 #START checking blast db
 #START lighter
 #START estimate k-mer size
 #k-mer size: 5964136
 #/home/jhuang/miniconda3/envs/vrap/bin/lighter -r ./results_ATCC19606/trimmed/A6_WT_HQ_trimmed_P_1.fastq.gz -r ./results_ATCC19606/trimmed/A6_WT_HQ_trimmed_P_2.fastq.gz -K 20 5964136 -od /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/lighter -t 55 2>> /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/vrap.log
 #START flash
 #/home/jhuang/Tools/vrap/external_tools/FLASH-1.2.11/flash /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/lighter/A6_WT_HQ_trimmed_P_1.cor.fq.gz /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/lighter/A6_WT_HQ_trimmed_P_2.cor.fq.gz -o flash -d /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/flash -M 1000 >> /mnt/h1/jhuang/DATA/Data_Tam_Acinetobacter_baumannii/A6_vrap_out/vrap.log
 #START spades

vrap.py (code)

 #!/usr/bin/env python3

 from orf import ORF
 from optparse import OptionParser
 import os
 import glob
 from Bio import SeqIO
 from pathlib import Path
 # see blast.py
 # .ncbirc --> BLASTDB=/media/jhuang/Titisee/GAMOLA2/Blast_db/ --> /mnt/Samsung_T5/Blast_db/
 from blast import Blast
 from hmmer import HMMER
 from bowtie import Bowtie
 # download_db.py only for txid2697049 [Severe acute respiratory syndrome coronavirus 2]
 from download_db import Database
 import glob
 import logging
 from math import exp
 from scipy.optimize import brentq

 #old = "_2010"
 old = ""

 class VrAP:
     """ VrAP """

     def __init__ (self):
         """ Class initialiser """
         self.version = "1.0.5"
         option = self.option_parser()
         self.path = os.path.dirname(os.path.realpath(__file__))

         self.paired_end = True
         self.p1  = option.p1
         self.p2  = option.p2
         self.single = ""
         self.cor1 = self.p1
         self.cor2 = self.p2
         self.cpu = option.cpu
         self.contigs = ""
         self.contig_length = option.clength
         self.host = option.host
         self.virus = option.virus
         self.nt = option.nt
         self.nr = option.nr
         self.bt2idx = option.bt2idx
         self.update = option.update
         self.gbt2 = option.gbowtie
         self.pre = option.pre
         self.assembly = option.assembly
         self.noblast = option.noblast

         if option.out!=None:
             self.out = os.path.abspath(option.out)
         else:
             self.out = "./"

         if self.update:
             self.update_database()

         if self.cpu==None:
             self.cpu = 1

         if self.p2==None:
             self.single = self.p1
             self.paired_end = False

         self.create_output_dir()

         open("{}/vrap.log".format(self.out),"w").close()
         self.log("START VrAP")

         with open("{}/vrap.log".format(self.out),"a") as log_opt:
             for opt, value in option.__dict__.items():
                 log_opt.write("{}:\t{}\n".format(opt,value))

         if self.pre:
             self.log("START only read correction and assembly!")

         self.check_input()

     def check_input (self):
         """ check input """

         bt2idx = self.bt2idx
         if self.bt2idx!=None:
             bt2idx = self.bt2idx+".1.bt2"

         input_files = [self.p1,self.p2,self.host,self.virus,bt2idx]

         for f in input_files:
             #print(f)
             if f != None:
                 if not os.path.exists(f):

                     print("{} don't exists\nPlease check input files\nVrAP quit\n".format(f))
                     quit()

         if self.p1 == None:
             self.parser.print_help()
             print("-1/--p1 is mandatory\nPlease check input file\nVrAP quit\n")
             quit()

         #EXPLANATION: in total, mapping to 4 dbs: viral_nucleotide, viral_protein, nt, nr
         dbs = ["{path}/database/viral_db/viral_nucleotide".format(path=self.path),"{path}/database/viral_db/viral_protein".format(path=self.path),self.nr,self.nt]
         status = []
         self.log("START checking blast db")

         for db in dbs:
             x = os.system("{path}/external_tools/blast/blastdbcmd -db {db} -info >> {out}/vrap.log 2>> {out}/vrap.log".format(path=self.path,db=db,out=self.out))
             status.append(x)

         #COMMENTED
         if status[0] != 0 or status[1] != 0:
             self.update_database()

         if self.nr != None and status[2] != 0:
             print("No alias or index file found for {}".format(self.nr))
             quit()

         if self.nt != None and status[3] != 0:
             print("No alias or index file found for {}".format(self.nt))
             quit()

     def log (self,message,err=False):
         """ log bowtie run """
         print(message)
         LOG_FILENAME = "{}/vrap.log".format(self.out)
         FORMAT = '\n[VrAP] %(asctime)-15s %(message)s\n'

         logging.basicConfig(filename=LOG_FILENAME,format=FORMAT,level=logging.DEBUG)
         if err:
             logging.exception(message)
         else:
             logging.info(message)
         #logging.exception()

     def option_parser(self):
         """ Read input options """
         usage = "usage: %prog -1 [FILE] [options]"
         parser = OptionParser(usage=usage,version="%prog {}".format(self.version))

         parser.add_option("-1","--p1", dest="p1",help="first paired end file" , metavar="FILE")
         parser.add_option("-2","--p2", dest="p2",help="second paired end file", metavar="FILE")
         parser.add_option("-t","--threads", dest="cpu",help="number of threads to use", metavar="INT",type="int",default=1)
         parser.add_option("-r","--host", dest="host",help="reference host sequences", metavar="FASTA")
         parser.add_option("-v","--virus", dest="virus",help="reference viral sequences", metavar="FASTA")
         parser.add_option("-n","--nt", dest="nt",help="nucleotide blast database (e.g. nt)", metavar="DATABASE")
         parser.add_option("-a","--nr", dest="nr",help="protein blast database (e.g. nr)", metavar="DATABASE")
         parser.add_option("-o","--outdir", dest="out",help="output directory", metavar="DIR")
         parser.add_option("-b","--bt2idx", dest="bt2idx",help="bowtie2 index", metavar="FILE")
         parser.add_option("-l","--clength", dest="clength",help="minimum contig length", metavar="INT",default=500)
         parser.add_option("-u","--update", dest="update",help="update viral database",  default=False, action="store_true")
         parser.add_option("-g","--gbt2", dest="gbowtie",help="use global PATH of bowtie2",  default=False, action="store_true")
         parser.add_option("-p","--pre", dest="pre",help="read correction only",  default=False, action="store_true")
         parser.add_option("-q","--assembly", dest="assembly",help="assembly only",  default=False, action="store_true")
         parser.add_option("-k","--noblast", dest="noblast",help="skip blast search",  default=False, action="store_true")
         self.parser = parser
         (options, args) = parser.parse_args()

         return options

     def create_output_dir(self):
         """ Create output directory """
         if not os.path.exists(self.out):
             os.makedirs(self.out)

     def EMfit2(self,F0, f1, F1, k):
         epstol = 1e-8
         x = F1 / F0
         e = 0
         e0 = 1
         x0 = 0
         maxit = 1000
         it = 0
         while abs(x - x0) / x + abs(e - e0) > epstol and it < maxit:
             it = it + 1
             #print x, e, abs(x-x0)/x
             x0 = x
             e0 = e
             x1 = -100
             while abs(x - x1) / x > epstol / 2:
                 x1 = x
                 x  = F1 / F0 * (3 * k * (1 - exp(-x * e / (3 * k))) + 1 - exp(-x * (1 - e)))
             def func(t):
                 return f1 / F1 - (t * exp(-x * t / (3 * k)) + (1 - t) * exp(-x * (1 - t)))
             print(func(0))
             print(func(1.2))
             e = brentq(func, 0, 1.1)
         return x, e

     def genome_size (self):
         """ estimate genome size """
         self.log("START estimate k-mer size")

         lighter_path = "{}/lighter".format(self.out)
         if not os.path.exists(lighter_path):
             os.makedirs(lighter_path)

         if self.paired_end:
             input_data = "{} {}".format(os.path.abspath(self.p1),os.path.abspath(self.p2),self.out)
         else:
             input_data = "{}".format(os.path.abspath(self.p1))

         command = "{path}/external_tools/KmerStream-master/KmerStream --tsv -k 20 -t {cpu} {input_data} -o {lighter}/kmer  2>> {out}/vrap.log".format(path=self.path,lighter=lighter_path,cpu=self.cpu,out=self.out,input_data=input_data)

         os.system(command)
         size = 30000

         try:
             kmer = "{lighter}/kmer".format(lighter=lighter_path)
             if os.path.exists(kmer):
                 with open(kmer) as f:
                     head = next(f)
                     value = next(f).split()

                 value = [int(i) for i in value] #list comprehension
                 size = value[2]-value[3]
         except Exception as e:
             vrap.log(e,True)

         self.log("k-mer size: {}".format(size))

         return size

     def lighter(self):
         """ Run ligther for read correction """
         self.log("START lighter")

         if self.paired_end:
             input_data = "-r {} -r {}".format(self.p1,self.p2)
         else:
             input_data = "-r {}".format(self.p1)

         size = self.genome_size()
         if size<30000:
             size = 30000
         #MODIFIED, since lighter in external_tool doesn't work.
         #{path}/external_tools/Lighter-master/lighter
         command = "/home/jhuang/anaconda3/bin/lighter {input} -K 20 {size} -od {out}/lighter -t {cpu} 2>> {out}/vrap.log".format(path=self.path,input=input_data,out=self.out,cpu=self.cpu,size=size)
         print(command)
         os.system(command)

         if self.paired_end:
             b1 = os.path.basename(self.p1).split(".")[0]
             self.cor1 = glob.glob("{}/lighter/{}*cor*".format(self.out,b1))[0]
             b2 = os.path.basename(self.p2).split(".")[0]
             self.cor2 = glob.glob("{}/lighter/{}*cor*".format(self.out,b2))[0]
         else:
             b1 = os.path.basename(self.single).split(".")[0]
             self.single = glob.glob("{}/lighter/{}*cor*".format(self.out,b1))[0]

     def flash (self):
         """ Run flash for paired end reads """
         self.log("START flash")

         if self.paired_end:
             command = "{path}/external_tools/FLASH-1.2.11/flash {cor1} {cor2} -o flash -d {out}/flash -M 1000 >> {out}/vrap.log".format(path=self.path,cor1=self.cor1,cor2=self.cor2,out=self.out)
             print(command)
             os.system(command)

             self.cor1   = "{}/flash/flash.notCombined_1.fastq".format(self.out)
             self.cor2   = "{}/flash/flash.notCombined_2.fastq".format(self.out)
             self.single = "{}/flash/flash.extendedFrags.fastq".format(self.out)

     def bowtie (self):
         """ Map reads to reference sequence """
         bowtie2 = Bowtie(self.gbt2,self.out,self.cpu)

         bowtie_path = "{}/bowtie".format(self.out)
         if not os.path.exists(bowtie_path):
             os.makedirs(bowtie_path)

         if self.bt2idx==None and self.host!=None:
             self.log("START bowtie2-build")
             self.bt2idx = "{}/bowtie/host".format(self.out)
             bowtie2.build(self.host)

         if self.bt2idx != None:
             if self.paired_end:
                 self.log("START bowtie2 pe")
                 bowtie2.pe(self.cor1,self.cor2,self.single,self.bt2idx)
                 self.cor1 = "{}/bowtie/bowtie.un.1.fastq".format(self.out)
                 self.cor2 = "{}/bowtie/bowtie.un.2.fastq".format(self.out)
                 self.single = "{}/bowtie/bowtie.un.fastq".format(self.out)
             else:
                 self.log("START bowtie2 se")
                 bowtie2.se(self.single,self.bt2idx)
                 self.single = "{}/bowtie/bowtie.un.fastq".format(self.out)

     def update_database (self):
         """ Update viral DB"""
         #self.log("START update viral DB")
         #EXPLANATION: maybe it needs to DELETE the directory viral_db.
         #typs = ["nucleotide"]
         #typs = ["protein"]
         typs = ["nucleotide","protein"]
         #path = os.path.dirname(os.path.realpath(__file__))
         db_path = "{}/database/viral_db/".format(self.path)

         if not os.path.exists(db_path):
                 os.makedirs(db_path)

         for typ in typs:
             self.log("START update viral {} DB".format(typ))
             database = Database(typ,db_path,self.cpu)

             database.download_id()
             database.start_download()

             blast = Blast(self.cpu, self.out)
             db_type = "nucl"
             if "protein" in typ:
                 db_type = "prot"

             self.log("\nSTART {} makeblastdb".format(typ))
             blast.makeblastdb(db_path+typ+".fa", "{}viral_{}".format(db_path,typ), db_type)
             #print(db_path+typ+".fa", "{}viral_{}".format(db_path,typ), db_type)
             #/home/jhuang/Tools/vrap/database/viral_db/nucleotide.fa /home/jhuang/Tools/vrap/database/viral_db/viral_nucleotide nucl
             #MANUAL_RUNNING
             #V5 didn't suitable
             #/home/jhuang/DATA/RNAHiSwitch/ncbi-blast-2.11.0+/bin/makeblastdb -in nuleotide.fa -dbtype nucl
             #/home/jhuang/DATA/RNAHiSwitch/ncbi-blast-2.11.0+/bin/makeblastdb -in protein.fa -dbtype prot
             #USING_V4
             #https://www.biostars.org/p/200742/
             #vrap/external_tools/blast/makeblastdb -in nucleotide.fa -dbtype nucl -parse_seqids -out virus_nucleotide  #-taxid_map nucl_gb.accession2taxid 
             #vrap/external_tools/blast/makeblastdb -in protein.fa -dbtype prot -parse_seqids -out virus_protein     #-taxid_map nucl_gb.accession2taxid
             #search_term="txid2697049[Organism:exp] NOT txid2[Organism:exp] NOT txid2759[Organism:exp] NOT txid2157[Organism:exp] {}".format(date)
             #virus_search_term="txid10239[Organism:exp] NOT txid2[Organism:exp] NOT txid2759[Organism:exp] NOT txid2157[Organism:exp] {}".format(date)
             #txid2=Bacteria, txid2759=Eukaryota, txid2157=Archaea, txid10239=Viruses, txid2697049=Severe acute respiratory syndrome coronavirus 2

     def spades (self):
         """ Run spades """
         self.log("START spades")
         paired = ""

         if self.paired_end:
             paired  = "-1 {} -2 {}".format(self.cor1,self.cor2)

         k       = "-k 33,55,77,99,127"
         options = "--cov-cutoff off --only-assembler --careful"
         #options = "--cov-cutoff off --only-assembler"

         single  = self.single
         if len(self.single)>0:
             single  = "--s1 {}".format(self.single)

         cpu     = "-t {}".format(self.cpu)
         output  = "-o {}/spades".format(self.out)
         os.system("{path}/external_tools/SPAdes-3.11.1-Linux/bin/spades.py {paired} {single} {k} {options} {cpu} {output} >> {out}/vrap.log".format(path=self.path, paired=paired, single=single, k=k, options=options, cpu=cpu, output=output,out=self.out))

         #for assembly only
         self.contigs = "{}/spades/contigs.fasta".format(self.out)
         self.exclude = set()
         self.scaffold = {}

     def cap3 (self):
         """ Runs cap3 """
         self.log("START cap3")
         self.contigs = "{}/spades/contigs.fasta".format(self.out)
         cap3_out = "{}/cap3".format(self.out)
         #print(cap3_out)
         if not os.path.exists(cap3_out):
             os.makedirs(cap3_out)

         os.system("{}/external_tools/cap3 '{}' -y 100 >{}/cap3/cap3.out".format(self.path,self.contigs,self.out))

         b = False
         i = 0
         consensi = {}
         consensus = ""
         group = {}
         name = ""
         exclude = set()

         with open("{}/cap3/cap3.out".format(self.out)) as f:
             for line in f:

                 if not b:

                     if "*******************" in line:
                         name = line.split("* ")[1].split(" *")[0]
                         #print(name)
                         group[name] = set()

                     if "NODE" in line:
                         ary = line.split()
                         for n in ary:
                             if "NODE" in n:
                                 tmp = n.replace("+","").replace("-","")
                                 group[name].add(tmp)
                                 exclude.add(tmp)

                 if "DETAILED" in line:
                     b = True

                 if b:

                     if "*******************" in line:
                         name = line.split("* ")[1].split(" *")[0]
                         #print(name)

                     ary = line.split()
                     #print(ary)
                     if len(ary) > 0:

                         if "consensus" in line:
                             seq = ary[1]
                             tmp = consensi.get(name,"")
                             tmp+=seq.replace("-","N")
                             consensi[name] = tmp

         self.scaffold = consensi
         self.exclude = exclude

     def clean_output (self):
         """ Clean spades output and include scaffolds """

         fasta_sequences = SeqIO.parse(open(self.contigs),'fasta')
         self.result = "{}/vrap_contig.fasta".format(self.out)
         output = open(self.result,"w")

         for fasta in fasta_sequences:
             if fasta.name not in self.exclude and len(fasta.seq)>=self.contig_length:
                 output.write(">{}\n{}\n".format(fasta.name,fasta.seq))
         i = 1
         for name,seq in self.scaffold.items():
             if len(seq)>=self.contig_length:
                 output.write(">CAP_{}_length_{}\n{}\n".format(i,len(seq),seq))
                 i+=1
         output.close()

     def calculate_orf_density (self):
         """ calculate orf density """
         self.log("START calculating orf density")
         self.result = "{}/vrap_contig.fasta".format(self.out)
         orf = ORF(self.result,self.out)
         self.orf_density = orf.orf_density()

     def label_as_virus (self):
         """ Add virus tag to name """
         out_name = "{}/blast/custom_viral_seq.fa".format(self.out)
         output = open(out_name,"w")

         with open(self.virus, "r") as handle:
             for record in SeqIO.parse(handle, "fasta"):
                 output.write(">{}\n{}\n".format(record.name+" [virus_user_db]",record.seq))

         output.close()

         return out_name

     def blast (self):
         """ blast """
         self.log("START blast")
         self.result = "{}/vrap_contig.fasta".format(self.out)

         blast_path = "{}/blast".format(self.out)
         blast = Blast(self.cpu, self.out)

         blast.read_fasta(self.result)

         #EXPLANATION: all virus and host-databases are saved in nt_dbs
         #             blast performed only after assembly(self.result=vrap_contig.fasta) against host+virus --> 
         #             since the host-reads are not deleted before the assembly, 
         #             vrap_contigs.fasta should contain contigs of host and virus
         #sort blast/blastn.csv and blast/blastn.fa
         nt_dbs = "{}/database/viral_db{}/viral_nucleotide".format(self.path,old)
         e_val = "1e-4"
         max_target = 1

         if not os.path.exists(blast_path):
             os.makedirs(blast_path)

         if not os.path.exists(blast_path+"/db"):
             os.makedirs(blast_path+"/db")

         #EXPLANATION
         #if the option --virus is set, make the given fasta blast-format and add it to vrap_out/blast/db/virus 
         #NO record added! nt_dbs should keep only 1 record!
         if self.virus != None:
             custom_viral_seq = self.label_as_virus()
             blast.makeblastdb(custom_viral_seq,"{}/blast/db/virus".format(self.out),"nucl")
             nt_dbs += " {}/blast/db/virus".format(self.out)
         #MODIFIED, add the self.host causing error. COMMENTED them!
         #EXPLANATION: if the option --host is set, make the given fasta blast-format and add it to vrap_out/blast/db/host
         #There are in total 3 nt_dbs, 1 prot_db, +[nt][nr], and 1 bowtie database:
         # nt_dbs /home/jhuang/Tools/vrap/database/viral_db/viral_nucleotide, ./blast/db/virus, ./blast/db/host
         # prot_db /home/jhuang/Tools/vrap/database/viral_db/viral_protein
         # bowtie_database: ./bowtie/host for substraction host-reads
         #DEBUG: /home/jhuang/Tools/vrap/external_tools/blast/blastn -num_threads 4 -query /home/jhuang/DATA/Data_Becher_Picornavirus_BovHepV/F21_1400_Liquor_vrap_out/vrap_contig.fasta -db "/home/jhuang/Tools/vrap/database/viral_db/viral_nucleotide /home/jhuang/DATA/Data_Becher_Picornavirus_BovHepV/F21_1400_Liquor_vrap_out/blast/db/host"  -evalue 1e-4 -outfmt "6 qseqid qstart qend sstart send evalue length pident sseqid stitle qcovs qcovhsp sacc slen qlen" -max_target_seqs 1 -out /home/jhuang/DATA/Data_Becher_Picornavirus_BovHepV/F21_1400_Liquor_vrap_out/blast/blastn.csv
         #--> it works!
         #Warning: [blastn] Examining 5 or more matches is recommended
         #grep -v "scrofa" blastn.csv > blastn_.csv;
         #DEBUG: BLAST Database creation error: Error: Duplicate seq_ids are found: GB|MF170920.1
         #if self.host != None:
         #   blast.makeblastdb(self.host,"{}/blast/db/host".format(self.out),"nucl")
         #   nt_dbs += " {}/blast/db/host".format(self.out)

         #EXPLANATION: blastn.csv is empty means blastn finds nonthing!
         ##nucleotide db search
         blast_out = "{}/blast/blastn.csv".format(self.out)
         blast.blastn(self.result, nt_dbs, e_val, max_target, blast_out)
         no_match_fa = "{}/blast/blastn.fa".format(self.out)
         no_results_found = blast.no_match(blast_out,no_match_fa)

         #EXPLANATION: blastx -num_threads 14 -query /home/jhuang/DATA/Data_Becher_Picornavirus_BovHepV/F21_1399_1400_Serum_vrap_out/blast/blastn.fa -db /home/jhuang/Tools/vrap/database/viral_db/viral_protein -evalue 1e-6 -outfmt 6 qseqid qstart qend sstart send evalue le+ 
         ##protein db search
         if no_results_found > 0:
             e_val = "1e-6"
             blast_out = "{}/blast/blastx.csv".format(self.out)
             prot_db = "{}/database/viral_db{}/viral_protein".format(self.path,old)
             blast.blastx(no_match_fa, prot_db, e_val, max_target, blast_out)
             no_match_fa = "{}/blast/blastx.fa".format(self.out)
             no_results_found = blast.no_match(blast_out,no_match_fa)

         #EXPLANATION: since self.nt and self.nr are not given, the following code won't run!
         if self.nt != None and no_results_found > 0:
             e_val = "1e-4"
             blast_out = "{}/blast/blastn_nt.csv".format(self.out)
             blast.blastn(no_match_fa, self.nt, e_val, max_target, blast_out)
             no_match_fa = "{}/blast/blastn_nt.fa".format(self.out)
             no_results_found = blast.no_match(blast_out,no_match_fa)

         if self.nr != None and no_results_found > 0:
             e_val = "1e-6"
             blast_out = "{}/blast/blastx_nr.csv".format(self.out)
             blast.blastx(no_match_fa, self.nr, e_val, max_target, blast_out)
             no_match_fa = "{}/blast/blastx_nr.fa".format(self.out)
             no_results_found = blast.no_match(blast_out,no_match_fa)

     def hmmer (self):
         """ hidden markov modell search """
         self.log("START HMMER")
         hmmer = HMMER(self.cpu, self.out)
         hmmer.run_hmmer()

     def summary (self):
         """ summarize output """

         if len(self.orf_density)>0:
             blast = Blast(self.cpu, self.out)
             self.blast_results = blast.summary(self.out)

             hmmer = HMMER(self.cpu, self.out)
             hmmer_output = hmmer.read_hmmer_output()

         with open("{}/vrap_summary{}.csv".format(self.out,old),"w") as f:
             f.write("name;qleng;orf_d;hmmer_eval;hmm_model;virus?;ident;qcov;tcov;tlength;tid;tname;mean_eval\n")
             for name,orf_d in self.orf_density.items():
                 results = {"qlength":0,"tid":"-","tname":"-","tlength":0,"qcovarage":0,"ident":0,"tcovarage":0,"eval":1}

                 if name in self.blast_results:
                     results = self.blast_results[name]

                 hmmer = 1
                 hmm_model = ""

                 if name in hmmer_output:
                     hmmer,hmm_model = max(hmmer_output[name],key=lambda x:x[0])

                 qleng = int(name.split("_")[3])
                 virus = "host"

                 if "virus" in results["tname"] or "phage" in results["tname"]:
                     virus = "virus"
                 else:
                     if (len(results["tname"])<2):
                         if ((orf_d>0.9 and qleng <=1000) or (orf_d>0.7 and qleng >1000)):
                             virus = "potential virus"
                         else:
                             virus = "-"

                     if hmmer<=1e-5:
                         virus = "virus"

                 tmp ="{0};{1:d};{2:.2f};{3:.2e};{4};{5};{6:.2f};{7:.2f};{8:.2f};{9:d};{10};{11};{12:.2e}\n".format(name,qleng,orf_d,hmmer,hmm_model,virus,results["ident"],results["qcovarage"],results["tcovarage"],results["tlength"],results["tid"],results["tname"].replace(";",","),results["eval"])
                 f.write(tmp)

 vrap = VrAP()

 try:
     if not vrap.assembly:
         vrap.lighter()
         vrap.flash()

     if not vrap.pre:
         if not vrap.assembly:
             vrap.bowtie()

         vrap.spades()

         if not vrap.assembly:
             vrap.cap3()

         vrap.clean_output()

         vrap.calculate_orf_density()

         if len(vrap.orf_density)>0:
             vrap.hmmer()
             if not vrap.noblast:
                 vrap.blast()

         vrap.summary()

     vrap.log("VrAP pipeline finished.\nThank you for using VrAP!")
 except Exception as e:
     vrap.log(e,True)