add local hcp db module

Author: yangfangs

physpetool/phylotree/retrieveprotein.py

@@ -26,6 +26,7 @@
 
 
 """
+import shutil
 import glob
 import ftplib
 import os
@@ -100,7 +101,7 @@ def splist(l, s):
 
 
 
-def retrieveprotein(proindexlist, outpath, matchlist, spelist):
+def retrieveprotein(proindexlist, outpath, matchlist, spelist, local_db):
     """
     Retrieve proteins form Kegg DB
     :param proindexlist: a list contain protein index
@@ -115,33 +116,53 @@ def retrieveprotein(proindexlist, outpath, matchlist, spelist):
     if not os.path.exists(dirname):
         os.makedirs(dirname)
 
-    connect = ftplib.FTP("173.255.208.244")
-    connect.login('anonymous')
-    connect.cwd('/pub/databasehcp')
-    p = 1
-    # deal with none
-    # connect.mkd("tem/")
     fasta = {}
+    p = 1
 
-    for line in spelist:
-        w_file = dirname + "/" + line + ".fa"
-        fw_ = open(w_file, 'ab')
-        retrievename = line + '.fasta'
-        remoteFileName = 'RETR ' + os.path.basename(retrievename)
-        connect.retrbinary(remoteFileName, fw_.write)
-        fw_.write(b'\n')
-        fw_.close()
-        logretrieveprotein.info("Retrieve " + line + " highly conserved proteins completed")
-        # read get sequences
-        with open(w_file, 'r') as f:
-            for line in f:
-                if line != "\n":
-                    tem = line.strip()
-                    if tem[0] == '>':
-                        header = tem[1:]
-                    else:
-                        sequence = tem
-                        fasta[header] = fasta.get(header, '') + sequence
+    # get hcp proteins form ftp server
+    if local_db == "":
+        connect = ftplib.FTP("173.255.208.244")
+        connect.login('anonymous')
+        connect.cwd('/pub/databasehcp')
+
+        for line in spelist:
+            w_file = dirname + "/" + line + ".fa"
+            fw_ = open(w_file, 'ab')
+            retrievename = line + '.fasta'
+            remoteFileName = 'RETR ' + os.path.basename(retrievename)
+            connect.retrbinary(remoteFileName, fw_.write)
+            fw_.write(b'\n')
+            fw_.close()
+            logretrieveprotein.info("Retrieve " + line + " highly conserved proteins completed")
+            # read get sequences
+            with open(w_file, 'r') as f:
+                for line in f:
+                    if line != "\n":
+                        tem = line.strip()
+                        if tem[0] == '>':
+                            header = tem[1:]
+                        else:
+                            sequence = tem
+                            fasta[header] = fasta.get(header, '') + sequence
+    # get protein sequence from local
+    else:
+        for line in spelist:
+            file_name = line +".fasta"
+            file_name_new = line + ".fa"
+            abb_data_path = os.path.join(local_db,file_name)
+            abb_data_path_new = os.path.join(dirname,file_name_new)
+            shutil.copyfile(abb_data_path,abb_data_path_new)
+            logretrieveprotein.info("Retrieve " + line + " highly conserved proteins completed")
+            # read get sequences
+            with open(abb_data_path_new, 'r') as f:
+                for line in f:
+                    if line != "\n":
+                        tem = line.strip()
+                        if tem[0] == '>':
+                            header = tem[1:]
+                        else:
+                            sequence = tem
+                            fasta[header] = fasta.get(header, '') + sequence
 
     for index in proindexlist:
         have_none = False
@@ -184,7 +205,7 @@ def retrieveprotein(proindexlist, outpath, matchlist, spelist):
     return dirname
 
 
-def doretrieve(specieslistfile, outpath):
+def doretrieve(specieslistfile, outpath,local_db):
     '''main to  retrieve protein from kegg db'''
     # spelist = []
     # for line in specieslistfile:
@@ -194,7 +215,7 @@ def doretrieve(specieslistfile, outpath):
     logretrieveprotein.info("Reading organisms's names success!")
     colname = getcolname()
     proindexlist, matchlist = getspecies(spelist, colname)
-    dirpath = retrieveprotein(proindexlist, outpath, matchlist, spelist)
+    dirpath = retrieveprotein(proindexlist, outpath, matchlist, spelist,local_db)
     return dirpath
 
 
@@ -218,4 +239,4 @@ def hcp_name(index):
     #         print("http://rest.kegg.jp/get/" + proid + "/aaseq")
     specieslistfile = ['zma', "ath", "eco"]
     outpath = "/home/yangfang/test/alg2/"
-    doretrieve(specieslistfile, outpath)
+    doretrieve(specieslistfile, outpath,local_db="")

physpetool/physpe/autobuild.py

@@ -55,6 +55,8 @@
 mafftpara = "--auto"
 
 auto_build_log = getLogging('Used time')
+
+
 def start_args(input):
     """
 Argument parse
@@ -92,6 +94,9 @@ def start_args(input):
                                 default=False,
                                 help="The esrna mode is use SSU RNA sequence with extend SSU RNA sequence \
                                 (users provide) to reconstruct phylogenetic tree. ")
+    autobuild_args.add_argument('-db', action='store', dest="db",
+                                default="",
+                                help="The absolute path for local database.")
     advance_args.add_argument('--muscle', action='store_true', dest='muscle',
                               default=True,
                               help="Multiple sequence alignment by muscle. The default multiple sequence \
@@ -134,9 +139,9 @@ def starting(args):
 starting run reconstruct tree
     :param args: arguments
     """
-    print ("Loading organism's names success.....\n")
-    print ("The result are store in:{0}\n".format(args.outdata))
-    print ("Now loading data and constructing phylogenetic tree......")
+    print("Loading organism's names success.....\n")
+    print("The result are store in:{0}\n".format(args.outdata))
+    print("Now loading data and constructing phylogenetic tree......")
 
     in_put = args.spenames
 
@@ -151,7 +156,6 @@ def starting(args):
     else:
         pass
 
-
     # Reconstruct phylogenetic tree by ssu RNA.
     if args.ssurna:
         setlogdir(out_put)
@@ -178,7 +182,8 @@ def starting(args):
                       args.raxml, args.raxml_parameter,
                       args.fasttree, args.fasttree_parameter,
                       args.iqtree, args.iqtree_parameter,
-                      args.thread, args.extenddata)
+                      args.thread, args.extenddata,
+                      args.db)
 
     # Reconstruct phylogenetic tree by extend ssu rna method.
     elif args.essurna:
@@ -197,15 +202,16 @@ def starting(args):
     elif args.HCP:
         setlogdir(out_put)
         starting_hcp(in_put, out_put,
-                 args.muscle, args.muscle_parameter,
-                 args.clustalw, args.clustalw_parameter,
-                 args.mafft, args.mafft_parameter,
-                 args.gblocks, args.gblocks_parameter,
-                 args.trimal, args.trimal_parameter,
-                 args.raxml, args.raxml_parameter,
-                 args.fasttree, args.fasttree_parameter,
-                 args.iqtree, args.iqtree_parameter,
-                 args.thread)
+                     args.muscle, args.muscle_parameter,
+                     args.clustalw, args.clustalw_parameter,
+                     args.mafft, args.mafft_parameter,
+                     args.gblocks, args.gblocks_parameter,
+                     args.trimal, args.trimal_parameter,
+                     args.raxml, args.raxml_parameter,
+                     args.fasttree, args.fasttree_parameter,
+                     args.iqtree, args.iqtree_parameter,
+                     args.thread, args.db)
+
 
 def starting_hcp(in_put, out_put,
                  args_muscle, args_muscle_p,
@@ -215,25 +221,25 @@ def starting_hcp(in_put, out_put,
                  args_trimal, args_trimal_p,
                  args_raxml, args_raxml_p,
                  args_fasttree, args_fasttree_p,
-                 args_iqtree,args_iqtree_p,
-                 args_thread):
+                 args_iqtree, args_iqtree_p,
+                 args_thread, args_db):
     '''reconstruct phylogenetic tree by hcp method'''
     hcp_input, recovery_dic = checkKeggOrganism(in_put)
     start = time.time()
-    out_retrieve = doretrieve(hcp_input, out_put)
+    out_retrieve = doretrieve(hcp_input, out_put,args_db)
     end = time.time()
     auto_build_log.info('Retrieving HCP sequences used time: {} Seconds'.format(end - start))
 
     if not recovery_dic == {}:
-        recovery(out_retrieve,recovery_dic)
+        recovery(out_retrieve, recovery_dic)
     # set default aligned by muscle if not specify clustalw
     start2 = time.time()
     if args_clustalw:
         out_alg = doclustalw_file(out_retrieve, out_put, args_clustalw_p)
     elif args_mafft:
-        out_alg = domafft_file(out_retrieve, out_put, args_mafft_p,args_thread)
+        out_alg = domafft_file(out_retrieve, out_put, args_mafft_p, args_thread)
     elif args_muscle:
-        out_alg = domuscle_file(out_retrieve, out_put, args_muscle_p,args_thread)
+        out_alg = domuscle_file(out_retrieve, out_put, args_muscle_p, args_thread)
 
     out_concat = cocat_path(out_alg)
 
@@ -263,17 +269,16 @@ def starting_srna(in_put, out_put,
                   args_trimal, args_trimal_p,
                   args_raxml, args_raxml_p,
                   args_fasttree, args_fasttree_p,
-                  args_iqtree,args_iqtree_p,
+                  args_iqtree, args_iqtree_p,
                   args_thread):
     '''reconstruct phylogenetic tree by ssu rna method'''
-    ssu_input,recovery_dic = checkSilvaOrganism(in_put)
+    ssu_input, recovery_dic = checkSilvaOrganism(in_put)
     start = time.time()
     out_retrieve = retrieve16srna(ssu_input, out_put)
     end = time.time()
     auto_build_log.info('Retrieving SSU rRNA sequences used time: {} Seconds'.format(end - start))
     if not recovery_dic == []:
-        recovery_silva(out_retrieve,recovery_dic,ssu_input)
-
+        recovery_silva(out_retrieve, recovery_dic, ssu_input)
 
     # set default aligned by muscle if not specify clustalw or mafft
     start2 = time.time()
@@ -306,6 +311,7 @@ def starting_srna(in_put, out_put,
     end2 = time.time()
     auto_build_log.info('Constructing species tree used time: {} Seconds'.format(end2 - start2))
 
+
 def starting_ehcp(in_put, out_put,
                   args_muscle, args_muscle_p,
                   args_clustalw, args_clustalw_p,
@@ -314,17 +320,18 @@ def starting_ehcp(in_put, out_put,
                   args_trimal, args_trimal_p,
                   args_raxml, args_raxml_p,
                   args_fasttree, args_fasttree_p,
-                  args_iqtree,args_iqtree_p,
-                  args_thread, args_extenddata):
+                  args_iqtree, args_iqtree_p,
+                  args_thread, args_extenddata,
+                  args_db):
     '''reconstruct phylogenetic tree by ehcp method'''
     hcp_input, recovery_dic = checkKeggOrganism(in_put)
     start = time.time()
-    out_retrieve = doretrieve(hcp_input, out_put)
+    out_retrieve = doretrieve(hcp_input, out_put,args_db)
     end = time.time()
     auto_build_log.info('Retrieving HCP sequences used time: {} Seconds'.format(end - start))
 
     if not recovery_dic == {}:
-        recovery(out_retrieve,recovery_dic)
+        recovery(out_retrieve, recovery_dic)
     retrieve_pro = os.listdir(out_retrieve)
     for reline in retrieve_pro:
         fw_name = os.path.join(out_retrieve, reline)
@@ -341,7 +348,7 @@ def starting_ehcp(in_put, out_put,
     if args_clustalw:
         out_alg = doclustalw_file(out_retrieve, out_put, args_clustalw_p)
     elif args_mafft:
-        out_alg = domafft_file(out_retrieve, out_put, args_mafft_p,args_thread)
+        out_alg = domafft_file(out_retrieve, out_put, args_mafft_p, args_thread)
     elif args_muscle:
         out_alg = domuscle_file(out_retrieve, out_put, args_muscle_p, args_thread)
 
@@ -363,6 +370,7 @@ def starting_ehcp(in_put, out_put,
     end2 = time.time()
     auto_build_log.info('Constructing species tree used time: {} Seconds'.format(end2 - start2))
 
+
 def starting_esrna(in_put, out_put,
                    args_muscle, args_muscle_p,
                    args_clustalw, args_clustalw_p,
@@ -371,17 +379,17 @@ def starting_esrna(in_put, out_put,
                    args_trimal, args_trimal_p,
                    args_raxml, args_raxml_p,
                    args_fasttree, args_fasttree_p,
-                   args_iqtree,args_iqtree_p,
+                   args_iqtree, args_iqtree_p,
                    args_thread, args_extenddata):
     '''reconstruct phylogenetic tree by ssu rna extend method'''
     extend_check = checkFile(args_extenddata)
-    ssu_input,recovery_dic = checkSilvaOrganism(in_put)
+    ssu_input, recovery_dic = checkSilvaOrganism(in_put)
     start = time.time()
     out_retrieve = retrieve16srna(ssu_input, out_put)
     end = time.time()
     auto_build_log.info('Retrieving SSU rRNA sequences used time: {} Seconds'.format(end - start))
     if not recovery_dic == {}:
-        recovery(out_retrieve,recovery_dic)
+        recovery(out_retrieve, recovery_dic)
     retrieve_srna_path = os.path.join(out_retrieve, 'rna_sequence.fasta')
 
     fw = open(retrieve_srna_path, 'a')
@@ -418,4 +426,4 @@ def starting_esrna(in_put, out_put,
             args_raxml_p = raxmlpara_dna
             doraxml(out_f2p, out_put, args_raxml_p, args_thread)
     end2 = time.time()
-    auto_build_log.info('Constructing species tree used time: {} Seconds'.format(end2 - start2))
+    auto_build_log.info('Constructing species tree used time: {} Seconds'.format(end2 - start2))