Mabs 2.28

calculate_AG.py

@@ -89,7 +89,7 @@
 	s_number_of_busco_orthogroups_to_use = "1000" #сколько ортогрупп BUSCO использовать. Это строка, содержащая или число, или слово "all", если нужно использовать все. Если пользователь укажет больше, чем есть в используемой базе данных BUSCO, то calculate_AG всё равно будет использовать все.
 	s_maximum_allowed_intron_length = "from_BUSCO" #максимальная разрешённая длина интрона. По умолчанию, используется значение из файла dataset.cfg датасета BUSCO. Переменная начинается с "s_", потому что это строка. Ниже будет ещё переменная n_maximum_allowed_intron_length, которая число.
 
-	s_version_of_calculate_AG = "2.27" #версия этой программы. Всегда равна версии Mabs. Поскольку эта программа нужна, в первую очередь, для Mabs, то когда я увеличиваю номер версии Mabs, то увеличивается и номер версии calculate_AG, и наоборот.
+	s_version_of_calculate_AG = "2.28" #версия этой программы. Всегда равна версии Mabs. Поскольку эта программа нужна, в первую очередь, для Mabs, то когда я увеличиваю номер версии Mabs, то увеличивается и номер версии calculate_AG, и наоборот.
 
 	l_errors_in_command_line = [] #список ошибок в командной строке. Если пользователь совершил много ошибок, то calculate_AG напишет про них все, а не только про первую встреченную.
 
@@ -101,7 +101,7 @@
 2) --nanopore_reads        Path to Nanopore reads.
 3) --pacbio_clr_reads        Path to PacBio CLR reads, also known as "old PacBio" reads.
 4) --pacbio_hifi_reads       Path to PacBio HiFi reads, also known as CCS reads.
-5) --download_busco_dataset        Name of a file from http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the most taxonomically narrow dataset for your species. For example, for a human genome assembly, use "--download_busco_dataset primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz". Calculate_AG will download the respective file. This option is mutually exclusive with "--local_busco_dataset".
+5) --download_busco_dataset        Name of a file from http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the most taxonomically narrow dataset for your species. For example, for a human genome assembly use "--download_busco_dataset primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz". Calculate_AG will download the respective file. This option is mutually exclusive with "--local_busco_dataset".
 6) --threads        Number of CPU threads to be used by calculate_AG. The default value is 10.
 7) --use_proovframe         Whether calculate_AG should polish sequences using Proovframe. "true" or "false". The default value is "false".
 8) --output_folder        Output folder for calculate_AG results. The default is "AG_calculation_results".
@@ -388,7 +388,7 @@
 	#если файл с контигами пустой, то сразу останавливаю выполнение calculate_AG, считая что AG=0. Иначе Metaeuk выдаст ошибку (если я правильно помню).
 	n_size_of_the_file_with_contigs = os.stat(s_path_to_the_assembly).st_size
 	if n_size_of_the_file_with_contigs == 0:
-		f_log.write("AG is 0")
+		f_log.write("\nAG is 0. Number of genes in single-copy orthogroups is 0. Number of genes in true multicopy orthogroups is 0. Number of genes in false multicopy orthogroups is 0.\n")
 		f_AG_calculation_results.write("AG is 0")
 		sys.exit()
 
@@ -447,7 +447,7 @@
 
 	#если MetaEuk вообще не выдал результатов, то считаю, что AG=0.
 	if not os.path.exists(s_path_to_the_output_folder + "/MetaEuk_results.fas"):
-		f_log.write("AG is 0. Number of genes in single-copy orthogroups is 0. Number of genes in true multicopy orthogroups is 0. Number of genes in false multicopy orthogroups is 0.\n")
+		f_log.write("\nAG is 0. Number of genes in single-copy orthogroups is 0. Number of genes in true multicopy orthogroups is 0. Number of genes in false multicopy orthogroups is 0.\n")
 		f_AG_calculation_results.write("AG is 0")
 		sys.exit()
 
@@ -889,11 +889,11 @@
 
 			n_AG = n_number_of_single_copy_genes_found_in_the_assembly + n_number_of_true_multicopy_genes
 	else:
-		f_log.write("AG is 0. Number of genes in single-copy orthogroups is 0. Number of genes in true multicopy orthogroups is 0. Number of genes in false multicopy orthogroups is 0.\n")
+		f_log.write("\nAG is 0. Number of genes in single-copy orthogroups is 0. Number of genes in true multicopy orthogroups is 0. Number of genes in false multicopy orthogroups is 0.\n")
 		f_AG_calculation_results.write("AG is 0")
 		sys.exit()
 
-	f_log.write("AG is " + str(n_AG) + ". Number of genes in single-copy orthogroups is " + str(n_number_of_single_copy_genes_found_in_the_assembly) + ". Number of genes in true multicopy orthogroups is " + str(n_number_of_true_multicopy_genes) + ". Number of genes in false multicopy orthogroups is " + str(n_number_of_false_multicopy_genes) +".\n")
+	f_log.write("\nAG is " + str(n_AG) + ". Number of genes in single-copy orthogroups is " + str(n_number_of_single_copy_genes_found_in_the_assembly) + ". Number of genes in true multicopy orthogroups is " + str(n_number_of_true_multicopy_genes) + ". Number of genes in false multicopy orthogroups is " + str(n_number_of_false_multicopy_genes) +".\n")
 	f_AG_calculation_results.write("AG is " + str(n_AG))
 
 	f_log.close()

install.sh

@@ -58,25 +58,15 @@ fi
 #Copying the content of "Additional_src" to "Additional"
 cp -rp ./Additional_src ./Additional
 
-#Bedtools is pre-compiled, I just change permissions.
-chmod 755 ./Additional/Bedtools/bedtools
-
-#DIAMOND is pre-compiled, I just change permissions.
-chmod 755 ./Additional/DIAMOND/diamond
+#Changing permissions of all files to 755. If a user copied the files of Mabs to a Linux machine through a Windows machine, all files will lost their permissions, which will lead to errors.
+chmod --recursive 755 ./Additional/
 
 #Installing HMMER
 cd ./Additional/HMMER
-chmod 755 ./configure
 ./configure
 make
 cd ../..
 
-#MetaEuk is pre-compiled, I just change permissions. The pre-compiled version is for SSE4.1. Actually, there are MetaEuk versions for newer CPUs, but since MetaEuk is not a time-limiting step of Mabs, I don't provide them.
-chmod 755 ./Additional/MetaEuk/metaeuk-linux-sse41
-
-#Minimap2 is pre-compiled, I just change permissions.
-chmod 755 ./Additional/Minimap2/minimap2
-
 #Installing Modified_hifiasm. It is a special version of Hifiasm made for Mabs. Compared to the ordinary Hifiasm, it has an additional option "--only-primary" that forces an assembly to terminate after the GFA file with primary contigs has been made. This preliminary termination allows to save some time if only primary contigs are needed.
 cd ./Additional/Modified_hifiasm
 make
@@ -92,8 +82,5 @@ cd ./Additional/Flye
 make
 cd ../..
 
-#Proovframe is written in Perl, it does not require to be compiled. I just change permissions. The source code of Proovframe was slightly modified by me — mostly for Proovframe to be able to find DIAMOND provided with Mabs.
-chmod 755 ./Additional/Proovframe/bin/*
-
 #Making mabs-hifiasm.py, mabs-flye.py and calculate_AG.py executable
 chmod 755 mabs-hifiasm.py mabs-flye.py calculate_AG.py

mabs-flye.py

@@ -117,7 +117,7 @@
 	s_maximum_allowed_intron_length = "from_BUSCO" #максимальная разрешённая длина интрона. По умолчанию, используется значение из файла dataset.cfg датасета BUSCO.
 	s_additional_flye_parameters = "" #дополнительные параметры Flye.
 
-	s_Mabs_version = "2.27"
+	s_Mabs_version = "2.28"
 
 	l_errors_in_command_line = [] #список ошибок в командной строке. Если пользователь совершил много ошибок, то Mabs-flye напишет про них все, а не только про первую встреченную.
 
@@ -132,7 +132,7 @@
 
 [any of the above files may be in FASTQ or FASTA, gzipped or not]
 
-4) --download_busco_dataset        Name of a file from http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the most taxonomically narrow dataset for your species. For example, for a human genome assembly, use "--download_busco_dataset primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz". Mabs-flye will download the respective file. This option is mutually exclusive with "--local_busco_dataset".
+4) --download_busco_dataset        Name of a file from http://mikeshelk.site/Data/BUSCO_datasets/Latest/ . It should be the most taxonomically narrow dataset for your species. For example, for a human genome assembly use "--download_busco_dataset primates_odb10.2021-02-19.tar.gz" and for a drosophila genome assembly use "--download_busco_dataset diptera_odb10.2020-08-05.tar.gz". Mabs-flye will download the respective file. This option is mutually exclusive with "--local_busco_dataset".
 5) --threads        Number of CPU threads to be used by Mabs-flye. The default value is 10.
 6) --output_folder        Output folder for Mabs-flye results. The default is "Mabs_results".
 7) --number_of_busco_orthogroups        How many BUSCO orthogroups should Mabs-flye use. Should be either a positive integral value or "all" to use all orthogroups. The default value is 1000. 
@@ -768,12 +768,11 @@
 	Теперь нужно определить, какой опцией Flye давать риды (--nano-raw, --nano-corr или другие). Если пользователь хотя бы один файл дал в формате FASTA, то использую --nano-raw, за исключением случая, когда пользователь дал только риды HiFi — тогда использую --pacbio-hifi. Если все риды в формате FASTQ, то я делаю следующее:
 	I) Считаю точность для каждого рида, используя строку с качеством. "Точность" выражается в процентах.
 	II) Считаю медианное значение по значениям из "I)"
-	III) В завимисимости от значения из "II)" выбираю, какой опцией давать риды программе Flye. Задавая эти числа, я ориентировался на описания опций в https://github.com/fenderglass/Flye/blob/flye/docs/USAGE.md
+	III) В завимисимости от значения из "II)" выбираю, какой опцией давать риды программе Flye. Задавая эти числа, я ориентировался на описания опций в https://github.com/fenderglass/Flye/blob/flye/docs/USAGE.md , а также на мой опыт сборки по ридам Нанопора, имеющим высокую (около 99.2%) точность.
 	Соответствие между медианной точностью ридов, и выбранным режимом Flye:
 	(0; 95] - --nano-raw
-	(95; 97] - --nano-hq
-	(97; 99] - --nano-corr
-	(99; 100] - --pacbio-hifi
+	(95; 99.8] - --nano-hq
+	(99.8; 100] - --pacbio-hifi
 	Если среди файлов с ридами, которые пользователь дал программе, хотя бы один в формате FASTA, то Mabs пишет в логи "WARNING: you have provided reads in FASTA, while FASTQ is recommended. Using reads in FASTA may reduce the accuracy of the assembly."
 	
 	Для скорости, медианная точность считается только по ридам, относящимся к генам BUSCO.
@@ -823,11 +822,9 @@
 
 		if n_median_accuracy_of_reads <= 95:
 			s_flye_option_to_provide_reads_with = "--nano-raw"
-		elif (n_median_accuracy_of_reads > 95) and (n_median_accuracy_of_reads <= 97):
+		elif (n_median_accuracy_of_reads > 95) and (n_median_accuracy_of_reads <= 99.8):
 			s_flye_option_to_provide_reads_with = "--nano-hq"
-		elif (n_median_accuracy_of_reads > 97) and (n_median_accuracy_of_reads <= 99):
-			s_flye_option_to_provide_reads_with = "--nano-corr"
-		elif n_median_accuracy_of_reads > 99:
+		elif n_median_accuracy_of_reads > 99.8:
 			s_flye_option_to_provide_reads_with = "--pacbio-hifi"
 		else:
 			o_current_time_and_date = datetime.datetime.now()
@@ -885,6 +882,13 @@
 	n_number_of_the_point_under_analysis = 0 #порядковый номер точки, которую я анализирую. Считается от 1.
 
 	#Тут я описываю функцию, которой на вход даются значения assemble_ovlp_divergence и repeat_graph_ovlp_divergence, а выдаёт функция -AG. С минусом спереди — потому, что scipy.optimize.minimize ищет минимум, а не максимум. Поэтому для максимизации AG нужно минимизировать -AG. Два параметра, дающихся на вход, я даю через список (потому что так нужно scipy.optimize.minimize).
+
+	"""
+	Результаты, которые получаются у этой функции, я буду записывать в словарь словарей. В нём первый ключ это номер рассматриваемой точки, а второй ключ это или строка "assemble_ovlp_divergence", или строка "repeat_graph_ovlp_divergence", или строка "AG". Значением является, соответственно, или то, какой assemble_divergence_relative использовался в этой точке; или то, какой repeat_graph_ovlp_divergence использовался в этой точке; или то, какой AG получился в этой точке.
+	"""
+
+	dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic = {}
+
 	def function_two_Flye_parameters_to_minus_AG(l_two_input_parameters):
 		global n_number_of_the_point_under_analysis #Без этой строки возникает ошибка "local variable 'n_number_of_the_point_under_analysis' referenced before assignment", потому что без этой строки нельзя модифицировать ("n_number_of_the_point_under_analysis += 1") глобальную переменную внутри функции.
 
@@ -894,6 +898,11 @@ def function_two_Flye_parameters_to_minus_AG(l_two_input_parameters):
 
 		n_number_of_the_point_under_analysis += 1
 
+		dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis] = {}
+
+		dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["assemble_ovlp_divergence"] = n_assemble_ovlp_divergence
+		dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["repeat_graph_ovlp_divergence"] = n_repeat_graph_ovlp_divergence
+
 		o_current_time_and_date = datetime.datetime.now()
 		s_current_time_and_date = o_current_time_and_date.strftime("%H:%M:%S %Y-%m-%d")
 		f_log.write(s_current_time_and_date + "\n")
@@ -930,6 +939,8 @@ def function_two_Flye_parameters_to_minus_AG(l_two_input_parameters):
 		f_log.write(s_current_time_and_date + "\n")
 		f_log.write("AG for point " + str(n_number_of_the_point_under_analysis) + " is " + str(n_AG) + "\n\n")
 
+		dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["AG"] = n_AG
+
 		return(-n_AG)
 
 	#Теперь, собственно, делаю сборку, проверяя максимум n_maximum_number_of_points_to_try точек
@@ -944,9 +955,19 @@ def function_two_Flye_parameters_to_minus_AG(l_two_input_parameters):
 	"""
 	o_optimization_results = scipy.optimize.minimize(fun = function_two_Flye_parameters_to_minus_AG, x0 = [n_max_divergence_between_reads_during_disjointig_construction__when_Flye_is_run_with_default_parameters, n_repeat_graph_ovlp_divergence__when_Flye_is_run_with_default_parameters], method = "Nelder-Mead", bounds = ((0, 1), (0, 0.5)), options = {"maxfev" : n_maximum_number_of_points_to_try, "initial_simplex" : [[n_max_divergence_between_reads_during_disjointig_construction__when_Flye_is_run_with_default_parameters, n_repeat_graph_ovlp_divergence__when_Flye_is_run_with_default_parameters], [2 * n_max_divergence_between_reads_during_disjointig_construction__when_Flye_is_run_with_default_parameters, 2 * n_repeat_graph_ovlp_divergence__when_Flye_is_run_with_default_parameters], [n_max_divergence_between_reads_during_disjointig_construction__when_Flye_is_run_with_default_parameters, n_repeat_graph_ovlp_divergence__when_Flye_is_run_with_default_parameters / 2]]}) #Не понял, чем maxfev отличается от maxiter. Но оптимизация останавливается на n_maximum_number_of_points_to_try именно если указать это число как maxfev, а не как maxiter.
 
-	n_optimal_assemble_ovlp_divergence = o_optimization_results.x[0]
-	n_optimal_repeat_graph_ovlp_divergence = o_optimization_results.x[1]
-	n_maximum_AG = - int(o_optimization_results.fun) #конвертирую в int, потому что scipy.optimize.minimize выдаёт это число во float (хоть это всегда и целый float). А в логи я его хочу записать в виде int.
+	#Смотрю, какая из точек дала наиболее высокий AG. Если таких точек было несколько, то беру последнюю. Потому что, возможно, в процессе максимизации AG Mabs-flye приблизился к той комбинации assemble_ovlp_divergence и repeat_graph_ovlp_divergence, которая даёт наилучшую сборку.
+	n_maximum_AG = -100 #-100 это плейсхолдер.
+	n_optimal_assemble_ovlp_divergence = -100 #-100 это плейсхолдер
+	n_optimal_repeat_graph_ovlp_divergence = -100 #-100 это плейсхолдер
+	for n_number_of_the_point_under_analysis in range(1, n_maximum_number_of_points_to_try + 1):
+		n_assemble_ovlp_divergence = dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["assemble_ovlp_divergence"]
+		n_repeat_graph_ovlp_divergence = dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["repeat_graph_ovlp_divergence"]
+		n_AG = dd_point_number__and__point_characteristic_name__to_the_value_of_the_characteristic[n_number_of_the_point_under_analysis]["AG"]
+
+		if n_AG >= n_maximum_AG:
+			n_maximum_AG = n_AG
+			n_optimal_assemble_ovlp_divergence = n_assemble_ovlp_divergence
+			n_optimal_repeat_graph_ovlp_divergence = n_repeat_graph_ovlp_divergence
 
 	o_current_time_and_date = datetime.datetime.now()
 	s_current_time_and_date = o_current_time_and_date.strftime("%H:%M:%S %Y-%m-%d")