virannot_blast2tsv: otu.py comparison

comparison otu.py @ 4:bb29ae8708b5 draft default tip

planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/virAnnot commit 7036ce0e06b6dc64332b1a5642fc58928523c5c6

author	iuc
date	Tue, 13 May 2025 11:52:17 +0000
parents	f8ebd1e802d7
children

comparison

equal deleted inserted replaced

-:f8ebd1e802d7
+:bb29ae8708b5
 # Name: virAnnot_otu
 # Author: Marie Lefebvre - INRAE
 # Reuirements: Ete3 toolkit and external apps
-# Aims: Create viral OTUs based on RPS and Blast annotations
+"""Create viral OTUs based on RPS and Blast annotations"""
 import argparse
 import csv
 import logging as log
 import os
 startQ = int(row[5])
 endQ = int(row[6])
 frame = float(row[7])
 description = row[8]
 superkingdom = row[9]
+try:
+pident = row[10]
+except IndexError:
+log.info(rps_file[0])
+log.info(row)
 match = re.search("Viruses", superkingdom)
 # if contig is viral then retrieve sequence
-if match:
+if match and float(pident) >= options.viral_portion:
 options.fasta.sort()
 seq = _retrieve_fasta_seq(options.fasta[i][0], query_id)
 seq_length = len(seq)
 if endQ < seq_length:
 seq = seq[startQ - 1:endQ]
 if "nb" not in collection[cdd_id][query_id]:
 collection[cdd_id][query_id]["nb"] = 0
 if "taxonomy" not in collection[cdd_id][query_id]:
 collection[cdd_id][query_id]["taxonomy"] = "Unknown"
 else:
-log.info("No blast file")
+log.debug("No blast file")
 collection[cdd_id][query_id]["taxonomy"] = "Unknown"
 collection[cdd_id][query_id]["nb"] = 0
+# keep pfamXXX and RdRp 1
-collection[cdd_id]["short_description"] = description.split(",")[0] + description.split(",")[1]  # keep pfamXXX and RdRp 1
+collection[cdd_id]["short_description"] = description.split(",")[0] + description.split(",")[1]
 collection[cdd_id]["full_description"] = description
 i += 1
+if options.merge_rdrp == "yes":
+rdrp_list = ["pfam00680", "pfam02123", "pfam00978", "pfam00998"]
+collection["RdRp_merge"] = {}
+for cdd_id in collection:
+if cdd_id in rdrp_list and cdd_id != "RdRp_merge":
+log.info("Add " + cdd_id + " in merge")
+for query_id in collection[cdd_id]:
+if query_id not in collection["RdRp_merge"]:
+collection["RdRp_merge"][query_id] = {}
+collection["RdRp_merge"][query_id] = collection[cdd_id][query_id]
 return collection
 def _retrieve_fasta_seq(fasta_file, query_id):
 """
 log.info("Align sequences")
 if not os.path.exists(options.output):
 os.mkdir(options.output)
 color_by_sample = {}
 for cdd_id in hits_collection:
-cdd_output = options.output + "/" + hits_collection[cdd_id]["short_description"].replace(" ", "_")
+log.info("align seq for " + cdd_id)
+if cdd_id == "RdRp_merge":
+cdd_output = options.output + "/" + cdd_id
+else:
+cdd_output = options.output + "/" + hits_collection[cdd_id]["short_description"].replace(" ", "_")
 if not os.path.exists(cdd_output):
 os.mkdir(cdd_output)
 if os.path.exists(cdd_output + "/seq_to_align.fasta"):
 os.remove(cdd_output + "/seq_to_align.fasta")
 if os.path.exists(cdd_output + "/seq_nucc.fasta"):
 # create tree plot with colors
 file_matrix = cdd_output + "/identity_matrix.csv"
 log.info("Create tree...")
 _create_tree(tree_file, file_seq_aligned, tree_file + '.png', file_color_config)
-_compute_pairwise_distance(options, file_seq_aligned, file_matrix, cdd_id)
+_compute_pairwise_distance(file_seq_aligned, file_matrix, cdd_id)
 log.info("Retrieve OTUs...")
 # if os.path.exists(file_cluster):
 #     os.remove(file_cluster)
 otu_cmd = os.path.join(options.tool_path, 'seek_otu.R') + ' ' + file_matrix + ' ' + file_cluster + ' ' + str(options.perc)
 log.debug(otu_cmd)
 f = open(file_cluster, "w+")
 f.write('OTU_1,1,' + list(hits_collection[cdd_id].keys())[0] + ',')
 f.close()
-def _compute_pairwise_distance(options, file_seq_aligned, file_matrix, cdd_id):
+def _compute_pairwise_distance(file_seq_aligned, file_matrix, cdd_id):
 """
 Calculate paiwise distance between aligned protein sequences
 from a cdd_id
 """
 log.info("Compute pairwise distance of " + cdd_id)
 file_xlsx = options.output + '/otu_stats.xlsx'  # Create a workbook
 workbook = xlsxwriter.Workbook(file_xlsx)
 log.info("Writing stats to " + file_xlsx)
 for cdd_id in hits_collection:
 otu_collection = {}
-cdd_output = options.output + "/" + hits_collection[cdd_id]["short_description"].replace(" ", "_")
+if cdd_id == "RdRp_merge":
-worksheet = workbook.add_worksheet(hits_collection[cdd_id]["short_description"])  # add a worksheet
+cdd_output = options.output + "/" + cdd_id
+worksheet = workbook.add_worksheet(cdd_id)
+else:
+cdd_output = options.output + "/" + hits_collection[cdd_id]["short_description"].replace(" ", "_")
+worksheet = workbook.add_worksheet(hits_collection[cdd_id]["short_description"])  # add a worksheet
 file_cluster = cdd_output + '/otu_cluster.csv'
 file_fasta_nucc = cdd_output + '/representative_nucc.fasta'
 with open(file_cluster, 'r') as clust:
 otu_reader = csv.reader(clust, delimiter=',')
 samples_list = []
 samples_list.append(sample) if sample not in samples_list else samples_list
 if sample not in otu_collection[row[0]]:
 otu_collection[row[0]][sample] = {}
 otu_collection[row[0]][sample][contig] = {}
 # add read number of the contig and annotation
-if 'nb' in hits_collection[cdd_id][contig]:
+if contig in hits_collection[cdd_id]:
-otu_collection[row[0]][sample][contig]['nb'] = hits_collection[cdd_id][contig]["nb"]
+if 'nb' in hits_collection[cdd_id][contig]:
+otu_collection[row[0]][sample][contig]['nb'] = hits_collection[cdd_id][contig]["nb"]
+else:
+otu_collection[row[0]][sample][contig]['nb'] = 0
+if 'taxonomy' in hits_collection[cdd_id][contig]:
+otu_collection[row[0]][sample][contig]['taxonomy'] = hits_collection[cdd_id][contig]["taxonomy"]
+else:
+otu_collection[row[0]][sample][contig]['taxonomy'] = 'unknown'
 else:
-otu_collection[row[0]][sample][contig]['nb'] = 0
+otu_collection[row[0]][sample][contig] = {'nb': 0, 'taxonomy': 'unknown'}
-if 'taxonomy' in hits_collection[cdd_id][contig]:
-otu_collection[row[0]][sample][contig]['taxonomy'] = hits_collection[cdd_id][contig]["taxonomy"]
-else:
-otu_collection[row[0]][sample][contig]['taxonomy'] = 'unknown'
 else:
 otu_collection[row[0]][sample][contig] = {}
 # add read number of the contig and annotation
-if 'nb' in hits_collection[cdd_id][contig]:
+if contig in hits_collection[cdd_id]:
-otu_collection[row[0]][sample][contig]['nb'] = hits_collection[cdd_id][contig]["nb"]
+if 'nb' in hits_collection[cdd_id][contig]:
+otu_collection[row[0]][sample][contig]['nb'] = hits_collection[cdd_id][contig]["nb"]
+else:
+otu_collection[row[0]][sample][contig]['nb'] = 0
+if 'taxonomy' in hits_collection[cdd_id][contig]:
+otu_collection[row[0]][sample][contig]['taxonomy'] = hits_collection[cdd_id][contig]["taxonomy"]
+else:
+otu_collection[row[0]][sample][contig]['taxonomy'] = 'unknown'
 else:
-otu_collection[row[0]][sample][contig]['nb'] = 0
+otu_collection[row[0]][sample][contig] = {'nb': 0, 'taxonomy': 'unknown'}
-if 'taxonomy' in hits_collection[cdd_id][contig]:
-otu_collection[row[0]][sample][contig]['taxonomy'] = hits_collection[cdd_id][contig]["taxonomy"]
-else:
-otu_collection[row[0]][sample][contig]['taxonomy'] = 'unknown'
 if 'taxonomy' in hits_collection[cdd_id][contig]:
 otu_collection[row[0]]['global_taxonomy'] = hits_collection[cdd_id][contig]["taxonomy"]
 else:
 otu_collection[row[0]]['global_taxonomy'] = 'unknown'
 worksheet.write(row, column + 2, 'contigs_list')
 row = 1
 # column = 0
 with open(file_fasta_nucc, "w+") as f_nucc:
 for otu in otu_collection:
-log.info(otu)
 if isinstance(otu_collection[otu], dict):
 column = 0
 worksheet.write(row, column, otu)
 # prepare table with 0 in each cells
 for sample in otu_collection[otu]:
 with open(map_file_path, "w+") as map_file:
 headers = ['#cdd_id', 'align_files', 'tree_files', 'cluster_files', 'cluster_nb_reads_files', 'pairwise_files', 'description', 'full_description\n']
 map_file.write("\t".join(headers))
 for cdd_id in hits_collection:
-cdd_output = hits_collection[cdd_id]["short_description"].replace(" ", "_")
+if cdd_id == "RdRp_merge":
+cdd_output = "RdRp_merge"
+else:
+cdd_output = hits_collection[cdd_id]["short_description"].replace(" ", "_")
 short_description = cdd_output
 file_seq_aligned = cdd_output + '/seq_aligned.final_tree.fa'
 tree_file = cdd_output + '/tree.dnd.png'
 file_cluster = cdd_output + '/otu_cluster.csv'
 file_matrix = cdd_output + "/identity_matrix.csv"
 log.debug(html_cmd)
 os.system(html_cmd)
 def _set_options():
+"""
+Set parameters
+"""
 parser = argparse.ArgumentParser()
 parser.add_argument('-b', '--blast', help='TAB blast file from blast2ecsv module.', action='append', required=False, dest='blast', nargs='+')
 parser.add_argument('-r', '--rps', help='TAB rpsblast file from rps2ecsv module.', action='append', required=True, dest='rps', nargs='+')
 parser.add_argument('-f', '--fasta', help='FASTA file with contigs', action='append', required=True, dest='fasta', nargs='+')
 parser.add_argument('-p', '--percentage', help='Percentage similarity threshold for OTUs cutoff.', action='store', type=int, default=90, dest='perc')
 parser.add_argument('-vp', '--viral_portion', help='Minimun portion of viral sequences in RPS domain to be included.', action='store', type=float, default=0.3, dest='viral_portion')
 parser.add_argument('-mpl', '--min_protein_length', help='Minimum query protein length.', action='store', type=int, default=100, dest='min_protein_length')
+parser.add_argument('-m', '--merge_rdrp', help='Merge RdRp1, 2, 3 and 4 to create otu on it.', action='store', type=str, default="no", dest='merge_rdrp')
 parser.add_argument('-tp', '--tool_path', help='Path to otu_seek.R', action='store', type=str, default='./', dest='tool_path')
 parser.add_argument('-o', '--out', help='The output directory', action='store', type=str, default='./Rps2tree_OTU', dest='output')
 parser.add_argument('-rgb', '--rgb-conf', help='Color palette for contigs coloration', action='store', type=str, default='rgb.txt', dest='file_rgb')
 parser.add_argument('-v', '--verbosity', help='Verbose level', action='store', type=int, choices=[1, 2, 3, 4], default=1)
 args = parser.parse_args()
 return args
 def _set_log_level(verbosity):
+"""
+Debbug
+"""
 if verbosity == 1:
 log_format = '%(asctime)s %(levelname)-8s %(message)s'
 log.basicConfig(level=log.INFO, format=log_format)
 elif verbosity == 3:
 log_format = '%(filename)s:%(lineno)s - %(asctime)s %(levelname)-8s %(message)s'

Mercurial > repos > iuc > virannot_blast2tsv

comparison otu.py @ 4:bb29ae8708b5 draft default tip