From 3ea3b003267cd2b0f8f64a8cf0ac988fa5104f76 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 10:08:01 +0100 Subject: [PATCH 01/50] update version --- read2tree/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/__init__.py b/read2tree/__init__.py index 8e3e4b2..9063fd7 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -6,7 +6,7 @@ from pkg_resources import resource_string logging.getLogger(__name__).addHandler(logging.NullHandler()) -__version__ = '0.1.5' +__version__ = '0.2.0' __copyright__ = 'read2tree (C) 2017-{:d} David Dylus' \ .format(date.today().year) From f72d0e89b3f47357826bb94ebdb2f714ceeda047 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 13:57:11 +0100 Subject: [PATCH 02/50] add run file --- archive/run_r2t.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 archive/run_r2t.py diff --git a/archive/run_r2t.py b/archive/run_r2t.py new file mode 100644 index 0000000..163378d --- /dev/null +++ b/archive/run_r2t.py @@ -0,0 +1,16 @@ + + +import read2tree +from read2tree.main import main +from read2tree._utils import exe_name + +import sys + +print("start run_r2t 2 ") + +main(sys.argv[1:], exe_name=exe_name(), desc="descr") + + +print("finish run_r2t ") + +a=1 \ No newline at end of file From 2ddbd293d533054f8910b3f63172819259548655 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 13:57:44 +0100 Subject: [PATCH 03/50] [WIP] minimap2 & no temp file --- read2tree/Mapper.py | 119 ++++++++++++++++++++++++-------------------- 1 file changed, 65 insertions(+), 54 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 9348183..0b18655 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -40,6 +40,8 @@ from read2tree.FastxReader import FastxReader +samtools ="/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools " + class Mapper(object): """ Structure for reference @@ -97,43 +99,50 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): output_folder = os.path.join(self.args.output_path, "04_mapping_" + self._species_name) - if len(self._reads) == 2: - ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder.name) - if self.args.threads != None: - ngm_wrapper.options.options['-t'].set_value(self.args.threads) - ngm = ngm_wrapper() - bam_file = ngm['file'] - elif len(self._reads) != 2 and 'short' in self.args.read_type: - ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder.name) - if self.args.threads != None: - ngm_wrapper.options.options['-t'].set_value(self.args.threads) - ngm = ngm_wrapper() - bam_file = ngm['file'] - elif len(self._reads) != 2 and 'long' in self.args.read_type: - ngm_wrapper = NGMLR(ref_file_handle, reads, tmp_output_folder.name) - if self.args.threads != None: - ngm_wrapper.options.options['-t'].set_value(self.args.threads) - if self.args.ngmlr_parameters != None: - par = self.args.ngmlr_parameters.split(',') - ngm_wrapper.options.options['-x'].set_value(str(par[0])) - ngm_wrapper.options \ - .options['--subread-length'].set_value(int(par[1])) - ngm_wrapper.options.options['-R'].set_value(float(par[2])) - ngm = ngm_wrapper() - bam_file = ngm['file'] - self.logger.info('{}: Mapped {} / {} reads to {}'.format(self._species_name, ngm['reads_mapped'], - ngm['total_reads']+ngm['reads_mapped'], os.path.basename(ref_file_handle))) - self._rm_file(ref_file_handle + "-enc.2.ngm", ignore_error=True) - self._rm_file(ref_file_handle + "-ht-13-2.2.ngm", ignore_error=True) - self._rm_file(ref_file_handle + "-ht-13-2.3.ngm", ignore_error=True) - - end = time.time() - self.elapsed_time = end - start - self.logger.info('{}: Mapping to {} references took {}.' - .format(self._species_name, os.path.basename(ref_file_handle), - self.elapsed_time)) + bam_file= tmp_output_folder+"/"+ref_file_handle.split("/")[-1]+".bam" + minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " + + #self._output_shell(minimap2_ex+" -ax sr "+ ref_file_handle+ " -t " + str(self.args.threads) + " "+reads+" | "+samtools_ex+" view -F 4 -bh -S -t" + str(self.args.threads)+ " > " + bam_file) + self._output_shell(minimap2_ex + " -ax sr " + ref_file_handle + " -t " + str( + self.args.threads) + " " + reads + " > " + bam_file) + + # if len(self._reads) == 2: + # ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder) + # if self.args.threads != None: + # ngm_wrapper.options.options['-t'].set_value(self.args.threads) + # ngm = ngm_wrapper() + # bam_file = ngm['file'] + # elif len(self._reads) != 2 and 'short' in self.args.read_type: + # ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder) + # if self.args.threads != None: + # ngm_wrapper.options.options['-t'].set_value(self.args.threads) + # ngm = ngm_wrapper() + # bam_file = ngm['file'] + # elif len(self._reads) != 2 and 'long' in self.args.read_type: + # ngm_wrapper = NGMLR(ref_file_handle, reads, tmp_output_folder) + # if self.args.threads != None: + # ngm_wrapper.options.options['-t'].set_value(self.args.threads) + # if self.args.ngmlr_parameters != None: + # par = self.args.ngmlr_parameters.split(',') + # ngm_wrapper.options.options['-x'].set_value(str(par[0])) + # ngm_wrapper.options \ + # .options['--subread-length'].set_value(int(par[1])) + # ngm_wrapper.options.options['-R'].set_value(float(par[2])) + # ngm = ngm_wrapper() + # bam_file = ngm['file'] # this is a same file ! 'folder/output1/PYGNA_OGs.fa.sam' + # self.logger.info('{}: Mapped {} / {} reads to {}'.format(self._species_name, ngm['reads_mapped'], + # ngm['total_reads']+ngm['reads_mapped'], os.path.basename(ref_file_handle))) + # self._rm_file(ref_file_handle + "-enc.2.ngm", ignore_error=True) # ngmlr deosn't produce these files ! + # self._rm_file(ref_file_handle + "-ht-13-2.2.ngm", ignore_error=True) + # self._rm_file(ref_file_handle + "-ht-13-2.3.ngm", ignore_error=True) + # + # end = time.time() + # self.elapsed_time = end - start + # self.logger.info('{}: Mapping to {} references took {}.' + # .format(self._species_name, os.path.basename(ref_file_handle), + # self.elapsed_time)) - if ngm['reads_mapped'] > 0 and os.path.exists(bam_file) and os.path.getsize(bam_file) > 0: + if os.path.exists(bam_file) and os.path.getsize(bam_file) > 0: # ngm['reads_mapped'] > 0 and shutil.copy(bam_file, os.path.join(output_folder, os.path.basename(bam_file))) return self._post_process_read_mapping(ref_file_handle, bam_file) else: @@ -161,7 +170,7 @@ def _read_mapping_from_folder(self, mapping_name=None, ref_records=None): species+'_OGs.fa') map_reads_species[species] = Reference() self._output_shell( - 'samtools index -@ ' + str(self.args.threads) + ' ' + + samtools+' index -@ ' + str(self.args.threads) + ' ' + file) consensus = self._build_consensus_seq_v2(ref_file, file) records = [] @@ -230,16 +239,17 @@ def _make_tmpdir(self): accelerated when parsing large files from the node drive. ''' - try: - tmp_output_folder = tempfile.TemporaryDirectory( - prefix='ngm', dir=os.environ.get("TMPDIR")) - except NotADirectoryError: - self.logger.info('{}: Environmental variable TMPDIR not set, will use \ - native python tmpdir location.' - .format(self._species_name)) - else: - tmp_output_folder = tempfile.TemporaryDirectory(prefix='ngm_') - self.logger.debug('--- Creating tmp directory on local node ---') + # try: + # tmp_output_folder = tempfile.TemporaryDirectory( + # prefix='ngm', dir=os.environ.get("TMPDIR")) + # except NotADirectoryError: + # self.logger.info('{}: Environmental variable TMPDIR not set, will use \ + # native python tmpdir location.' + # .format(self._species_name)) + # else: + # tmp_output_folder = tempfile.TemporaryDirectory(prefix='ngm_') + # self.logger.debug('--- Creating tmp directory on local node ---') + tmp_output_folder=self.args.output_path return tmp_output_folder def _map_reads_to_references(self, ref): @@ -281,9 +291,10 @@ def _map_reads_to_references(self, ref): # write reference into temporary file ref_file_handle = os.path.join(reference_path, species+'_OGs.fa') - ref_tmp_file_handle = os.path.join(tmp_output_folder.name, - species + '_OGs.fa') - shutil.copy(ref_file_handle, ref_tmp_file_handle) + #ref_tmp_file_handle = os.path.join(tmp_output_folder, #.name, + # species + '_OGs.fa') + #shutil.copy(ref_file_handle, ref_tmp_file_handle) # it doesnt need to be temp. + ref_tmp_file_handle= ref_file_handle # call the WRAPPER here processed_reads = self._call_wrapper(ref_tmp_file_handle, reads, @@ -323,7 +334,7 @@ def _map_reads_to_references(self, ref): # self.progress.set_status('single_map', ref=species) self._rm_file(ref_file_handle+".fai", ignore_error=True) - tmp_output_folder.cleanup() + #tmp_output_folder.cleanup() end = time.time() self.elapsed_time = end - start if len(references) > 1: @@ -541,7 +552,7 @@ def _post_process_read_mapping(self, ref_file, bam_file): bam_file = sam_file.replace(".sam", ".bam") if os.path.exists(sam_file): self._output_shell( - 'samtools view -F 4 -bh -S -@ ' + str(self.args.threads) + + samtools+' view -F 4 -bh -S -@ ' + str(self.args.threads) + ' -o ' + bam_file + " " + sam_file) if self.args.single_mapping: self.logger.debug("{}: ---- Samtools view completed" @@ -549,7 +560,7 @@ def _post_process_read_mapping(self, ref_file, bam_file): if os.path.exists(bam_file): self._output_shell( - 'samtools sort -m 2G -@ ' + str(self.args.threads) + + samtools+' sort -m 2G -@ ' + str(self.args.threads) + # todo 2G might not be enough ' -o ' + outfile_name + "_sorted.bam " + bam_file) if self.args.single_mapping: self.logger.debug("{}: ---- Samtools sort completed" @@ -557,7 +568,7 @@ def _post_process_read_mapping(self, ref_file, bam_file): if os.path.exists(outfile_name + "_sorted.bam"): self._output_shell( - 'samtools index -@ ' + str(self.args.threads) + ' ' + + samtools+' index -@ ' + str(self.args.threads) + ' ' + outfile_name + "_sorted.bam") if self.args.single_mapping: self.logger.debug("{}: ---- Samtools index completed" From 524bc671c0425c43c25341976a7dce58f673a20f Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 14:11:05 +0100 Subject: [PATCH 04/50] [WIP] add short-ont-hifi for minimap2 --- read2tree/Mapper.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 0b18655..96bba37 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -39,7 +39,7 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader - +minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " samtools ="/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools " class Mapper(object): @@ -100,10 +100,17 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): "04_mapping_" + self._species_name) bam_file= tmp_output_folder+"/"+ref_file_handle.split("/")[-1]+".bam" - minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " + #self._output_shell(minimap2_ex+" -ax sr "+ ref_file_handle+ " -t " + str(self.args.threads) + " "+reads+" | "+samtools_ex+" view -F 4 -bh -S -t" + str(self.args.threads)+ " > " + bam_file) - self._output_shell(minimap2_ex + " -ax sr " + ref_file_handle + " -t " + str( + if 'short' in self.args.read_type: + minimap_argm =" -ax sr" + elif 'long-hifi' in self.args.read_type: + minimap_argm = " -ax map-hifi " + elif 'long-ont' in self.args.read_type: + minimap_argm = " -ax map-ont " + + self._output_shell(minimap2_ex + minimap_argm + ref_file_handle + " -t " + str( self.args.threads) + " " + reads + " > " + bam_file) # if len(self._reads) == 2: From 60591a4ccd169a974748c09ccfcd7f69625edb13 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 14:22:48 +0100 Subject: [PATCH 05/50] [WIP] fix --- read2tree/Mapper.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 96bba37..8e24e36 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -39,8 +39,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " -samtools ="/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools " +minimap2_ex= "minimap2" #"/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " +samtools = "samtools" # "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools " class Mapper(object): """ @@ -109,9 +109,12 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): minimap_argm = " -ax map-hifi " elif 'long-ont' in self.args.read_type: minimap_argm = " -ax map-ont " + line_minimp= minimap2_ex +" "+ minimap_argm + ref_file_handle + " -t " + str(self.args.threads) + " " + reads + " > " + bam_file + self._output_shell(line_minimp) + + self.logger.info('mapping with' + line_minimp) + - self._output_shell(minimap2_ex + minimap_argm + ref_file_handle + " -t " + str( - self.args.threads) + " " + reads + " > " + bam_file) # if len(self._reads) == 2: # ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder) From 842f6c373a2f215aaab87c0d06f58f0b701f7058 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:35:22 +0100 Subject: [PATCH 06/50] add comment --- read2tree/stats/Coverage.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/stats/Coverage.py b/read2tree/stats/Coverage.py index e0ac291..ed34723 100644 --- a/read2tree/stats/Coverage.py +++ b/read2tree/stats/Coverage.py @@ -12,7 +12,7 @@ def get_coverage_bam(self, file_name): mybam = pysam.AlignmentFile(file_name, 'rb') for ref in mybam.references: self.coverage[self._get_clean_id(ref)] \ - = self._get_gene_coverage(mybam, ref) + = self._get_gene_coverage(mybam, ref) # {'ASTMX00297_OG1080037': [63.51897184822521, 30.055455801280203], def _get_clean_id(self, id): id = id.split(" ")[0] From 77dd2975d7f317a06fee1905137ccd90b1dfd498 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:37:44 +0100 Subject: [PATCH 07/50] [fix] using threads in minimap2 --- read2tree/Mapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 8e24e36..da6034b 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -109,7 +109,7 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): minimap_argm = " -ax map-hifi " elif 'long-ont' in self.args.read_type: minimap_argm = " -ax map-ont " - line_minimp= minimap2_ex +" "+ minimap_argm + ref_file_handle + " -t " + str(self.args.threads) + " " + reads + " > " + bam_file + line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) + ref_file_handle + " " + reads + " > " + bam_file self._output_shell(line_minimp) self.logger.info('mapping with' + line_minimp) From 7e7d74cd51b4f9db6ab53355c8c149a5b0f899fb Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:42:44 +0100 Subject: [PATCH 08/50] improve logging --- read2tree/parser/OMAOutputParser.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/read2tree/parser/OMAOutputParser.py b/read2tree/parser/OMAOutputParser.py index 90cba48..50488f6 100644 --- a/read2tree/parser/OMAOutputParser.py +++ b/read2tree/parser/OMAOutputParser.py @@ -125,7 +125,10 @@ def _filter_ogs_min_species(self): for x in ignored: print(" - {}".format(x)) if len(names_og) == 0: - raise Exception("could not load any marker genes from {}".format(self.og_fasta_path)) + if "/" in str(self.og_fasta_path): + raise Exception("could not load any marker genes from with path " + str(self.og_fasta_path)) + else: + raise Exception("could not load any marker genes from with path "+str(os.getcwd())+"/"+str(self.og_fasta_path)) return names_og def _estimate_best_number_species(self): From de8d16e5e5a95862094e7019a23c88d989e985bb Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:43:35 +0100 Subject: [PATCH 09/50] add pass --- archive/tests/test_use.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/archive/tests/test_use.py b/archive/tests/test_use.py index 0d87f79..e672129 100644 --- a/archive/tests/test_use.py +++ b/archive/tests/test_use.py @@ -10,10 +10,13 @@ class Use(unittest.TestCase): def test_OGSet(self): + pass def test_write_progress(self): + pass def test_read_progress(self): + pass if __name__ == "__main__": From 385ef09ff34d4e573c25e497150267c487f75ca0 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:52:10 +0100 Subject: [PATCH 10/50] add pass --- tests/test_use.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/tests/test_use.py b/tests/test_use.py index 0d87f79..e672129 100644 --- a/tests/test_use.py +++ b/tests/test_use.py @@ -10,10 +10,13 @@ class Use(unittest.TestCase): def test_OGSet(self): + pass def test_write_progress(self): + pass def test_read_progress(self): + pass if __name__ == "__main__": From aeabe5296d4b9c115c9ad243e41b01eb4d575ed7 Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:56:08 +0100 Subject: [PATCH 11/50] [WIP] update minimap2 options --- read2tree/main.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/read2tree/main.py b/read2tree/main.py index f7e6e33..f10c32b 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -65,12 +65,10 @@ def parse_args(argv, exe_name, desc): arg_parser.add_argument('--read_type', default='short', help='[Default is short] Type of reads to ' - 'use for mapping: short or long. Either ngm for short reads or ' - 'ngmlr for long will be used.') + 'use for mapping: short, long-hifi or long-ont corresponding to sr, map-hifi, or map-ont in minimap2.') arg_parser.add_argument('--threads', type=int, default=1, - help='[Default is 1] Number of threads for the mapping ' - 'using ngm / ngmlr!') + help='[Default is 1] Number of threads for the mapping ') arg_parser.add_argument('--split_reads', action='store_true', help='[Default is off] Splits reads as defined by split_len (200) ' @@ -117,13 +115,13 @@ def parse_args(argv, exe_name, desc): '(number of ACGT basepairs vs length of sequence). ' 'By default, all sequences are selected.') - arg_parser.add_argument('--ngmlr_parameters', default=None, - help='[Default is none] In case this parameters ' - 'need to be changed all 3 values have to be ' - 'changed [x,subread-length,R]. The standard ' - 'is: ont,256,0.25. Possibilities for these ' - 'parameter can be found in the original ' - 'documentation of ngmlr.') + # arg_parser.add_argument('--ngmlr_parameters', default=None, # todo this could be used for minimap2 options + # help='[Default is none] In case this parameters ' + # 'need to be changed all 3 values have to be ' + # 'changed [x,subread-length,R]. The standard ' + # 'is: ont,256,0.25. Possibilities for these ' + # 'parameter can be found in the original ' + # 'documentation of ngmlr.') arg_parser.add_argument('--check_mate_pairing', action='store_true', @@ -239,10 +237,10 @@ def parse_args(argv, exe_name, desc): arg_parser.error( 'Splitting reads does not work for paired end reads.') - if args.read_type == 'short' and args.ngmlr_parameters: - arg_parser.error( - 'Arguments for --ngmlr_parameters only work if --read_type is set ' - 'to "long".') + # if args.read_type == 'short' and args.ngmlr_parameters: + # arg_parser.error( + # 'Arguments for --ngmlr_parameters only work if --read_type is set ' + # 'to "long".') if not args.sample_reads and (args.coverage != 10 or args.genome_len != 2000000): From f602ca06baafeacce2a0896802b039e93b52433e Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 16:58:38 +0100 Subject: [PATCH 12/50] [WIP] cleanup --- read2tree/Mapper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index da6034b..5c146c0 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -39,8 +39,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "minimap2" #"/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2 " -samtools = "samtools" # "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools " +minimap2_ex= "minimap2" +samtools = "samtools" class Mapper(object): """ From d588782ffbeb553dae8d5cc4fa2d55388fb824af Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 17:36:42 +0100 Subject: [PATCH 13/50] [WIP] fix threads, improve logging --- read2tree/Mapper.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 5c146c0..59678ab 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -42,6 +42,7 @@ minimap2_ex= "minimap2" samtools = "samtools" + class Mapper(object): """ Structure for reference @@ -109,7 +110,7 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): minimap_argm = " -ax map-hifi " elif 'long-ont' in self.args.read_type: minimap_argm = " -ax map-ont " - line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) + ref_file_handle + " " + reads + " > " + bam_file + line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads + " > " + bam_file self._output_shell(line_minimp) self.logger.info('mapping with' + line_minimp) @@ -650,15 +651,17 @@ def _output_shell(self, line): :return: """ try: + self.logger.debug("Running " + line) shell_command = subprocess.Popen( line, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True) - except OSError: - return None - except ValueError: + except: + self.logger.debug("Shell command failed to execute by running ") return None (output, err) = shell_command.communicate() + self.logger.debug("Shell output: "+ str(output)) + self.logger.debug("Shell err: " + str(err)) shell_command.wait() if shell_command.returncode != 0: self.logger.debug("Shell command failed to execute") From 95b4776ce34a950ca79a1bc507ad35c8989c008d Mon Sep 17 00:00:00 2001 From: sina Date: Sat, 30 Mar 2024 17:39:28 +0100 Subject: [PATCH 14/50] [WIP] no need to limit memory --- read2tree/Mapper.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 59678ab..86c8759 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -571,7 +571,7 @@ def _post_process_read_mapping(self, ref_file, bam_file): if os.path.exists(bam_file): self._output_shell( - samtools+' sort -m 2G -@ ' + str(self.args.threads) + # todo 2G might not be enough + samtools+' sort -@ ' + str(self.args.threads) + # todo 2G might not be enough ' -o ' + outfile_name + "_sorted.bam " + bam_file) if self.args.single_mapping: self.logger.debug("{}: ---- Samtools sort completed" From 7b569f6cea34265b5b42e87e1a9de6a17c76d9dd Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 12:34:41 +0200 Subject: [PATCH 15/50] [WIP] add todo --- read2tree/Aligner.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/Aligner.py b/read2tree/Aligner.py index f85edf5..3efd5f3 100644 --- a/read2tree/Aligner.py +++ b/read2tree/Aligner.py @@ -178,7 +178,7 @@ def add_mapped_seq(self, ogset_add, species_name=None): self.mapped_aligns[name_og].dna = align_filt.dna end = time.time() self.elapsed_time = end-start - logger.info('{}: Appending {} reconstructed sequences to present Alignments ' + logger.info('{}: Appending {} reconstructed sequences to present Alignments ' # todo when is all cases are 0, needs better handling 'took {}.' .format(self._species_name, num_append_seq, From 9e47b69fe76ccc150aa5b213a4f7ecac56f73718 Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 12:35:37 +0200 Subject: [PATCH 16/50] improve logging in debug --- read2tree/Mapper.py | 4 ++++ read2tree/Progress.py | 26 +++++++++++++++++++++++--- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 86c8759..3608a4b 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -352,6 +352,10 @@ def _map_reads_to_references(self, ref): self.logger.info('{}: Mapping to all references took {}.' .format(self._species_name, self.elapsed_time)) + if not mapped_reads_species : + self.logger.warning('mapped_reads_species is empty') + + return mapped_reads_species def _write_read_query_aling(self, read, og_name_file, write_mode): diff --git a/read2tree/Progress.py b/read2tree/Progress.py index 85c3406..cfadf22 100644 --- a/read2tree/Progress.py +++ b/read2tree/Progress.py @@ -87,7 +87,7 @@ def _extract_line_from_log(self, word, logfile): with open(logfile, "r") as file: bestline = [line.split() for line in file if word in line] if bestline: - return bestline[-1] + return bestline[-1] # the last line with this word is selected return None except FileNotFoundError: print('File {} not accessible'.format(logfile)) @@ -101,6 +101,7 @@ def _get_number_of_OGs(self): ''' log_list = self._extract_line_from_log('Gathering', 'mplog.log') if log_list: + logging.debug(' We are using the info from line #' +" ".join(log_list[:3])+".# So number of OGs is "+str(int(log_list[13])) ) return int(log_list[13]) else: return 0 @@ -113,6 +114,7 @@ def _get_number_of_appeneded_seq_to_OGs(self): ''' log_list = self._extract_line_from_log('Appending', 'mplog.log') if log_list: + logging.debug(' We are using the info from line #' + " ".join(log_list[:3]) + ".# So number of appended sequences to OGs " + str(int(log_list[9]))) return int(log_list[9]) else: return 0 @@ -125,6 +127,7 @@ def _get_number_of_alignments(self): ''' log_list = self._extract_line_from_log('Alignment of', 'mplog.log') if log_list: + logging.debug(' We are using the info from line #' + " ".join(log_list[:3]) + ".# So number of alignments is " + str(int(log_list[10]))) return int(log_list[10]) else: return 0 @@ -135,9 +138,10 @@ def _get_number_of_references(self): 2018-11-23 12:13:53,691 - read2tree.ReferenceSet - INFO - ass: Extracted 6 reference species form 5 ogs took 0.0008709430694580078 :return: Number of reference species ''' - log_list = self._extract_line_from_log('ReferenceSet', 'mplog.log') + log_list = self._extract_line_from_log('ReferenceSet', 'mplog.log') # # the last line with ReferenceSet is selected if log_list: - return int(log_list[9]) + logging.debug("We are using the info from line #" + " ".join(log_list[:3]) + ".# So number of references is " + str(int(log_list[9]))) + return int(log_list[9]) # else: return 0 @@ -159,9 +163,12 @@ def _get_og_set_status(self): if os.path.exists(self._folder_ref_ogs_aa) and os.path.exists(self._folder_ref_ogs_dna): num_ogs_aa = self._count_files(self._folder_ref_ogs_aa, '*fa') num_ogs_dna = self._count_files(self._folder_ref_ogs_dna, '*fa') + logging.debug(' We are counting the number of fa files in folder _ogs_aa and _ogs_daa which are ' + str(num_ogs_aa) + " and "+ str(num_ogs_dna) +" in folder "+str(self._folder_ref_ogs_aa) +" and "+ str(self._folder_ref_ogs_dna)) if (num_ogs_expected-num_ogs_aa) == 0 and (num_ogs_expected-num_ogs_dna) == 0: + logging.debug(' We are counting the number of fa files in folder _ogs_aa and _ogs_daa, which are ' + str(num_ogs_aa) + " and " + str(num_ogs_dna) +" the same as expected"+str(num_ogs_expected) +". So this step is done") return True else: + logging.debug(' We are counting the number of fa files in folder _ogs_aa and _ogs_daa, which are ' + str(num_ogs_aa) + " and " + str(num_ogs_dna) +" but not the same as expected "+str(num_ogs_expected) +". So this step is done") return False else: return False @@ -176,8 +183,11 @@ def _get_append_og_set_status(self): num_ogs_aa = self._count_files(self._folder_append_og_aa, '*fa') num_ogs_dna = self._count_files(self._folder_append_og_dna, '*fa') if (num_ogs_expected-num_ogs_aa) <= 0 and (num_ogs_expected-num_ogs_dna) <= 0: + logging.debug("Number of ogs expected after appending is"+str(num_ogs_expected)+" the same as ogs number in dna and aa " + str(self._folder_ref_ogs_dna) +". So this step is done") return True else: + logging.debug("Number of ogs expected after appending is"+str(num_ogs_expected)+" but the number in dna and aa " + str(self._folder_ref_ogs_dna) +" are " + str(num_ogs_aa) + " and "+ str(num_ogs_dna)+". So this step is not done") + return False else: return False @@ -191,8 +201,10 @@ def _get_reference_status(self): if os.path.exists(self._folder_ref_dna): num_references = self._count_files(self._folder_ref_dna, '*fa') if (num_ref_expected-num_references) == 0: + logging.debug("Number of ref expected is " + str( num_ref_expected) + " and number in " + str( self._folder_ref_ogs_dna) +" is "+ str(num_references)+ ". So this step is done") return True else: + logging.debug("Number of ref expected is " + str( num_ref_expected) + " but number in " + str( self._folder_ref_ogs_dna) +" is "+ str(num_references)+ ". So this step is not done") return False else: return False @@ -207,8 +219,10 @@ def _get_alignment_status(self): num_align_aa = self._count_files(self._folder_align_aa, '*phy') num_align_dna = self._count_files(self._folder_align_dna, '*phy') if (num_aligns_expected-num_align_aa) == 0 and (num_aligns_expected-num_align_dna) == 0: + logging.debug("Number of aligns expected is " + str( num_aligns_expected) + " and number in " + str( self._folder_ref_ogs_dna) +" is "+ str(num_align_dna)+ ", similarly for aa. So this step is done") return True else: + logging.debug("Number of aligns expected is " + str( num_aligns_expected) + " and number in " + str( self._folder_ref_ogs_dna) +" and aa versions are "+ str(num_align_dna)+ " and "+ str(num_align_aa)+ ". So this step is not done") return False else: return False @@ -223,8 +237,10 @@ def _get_append_alignment_status(self): num_align_aa = self._count_files(self._folder_align_append_aa, '*phy') num_align_dna = self._count_files(self._folder_align_append_dna, '*phy') if (num_aligns_expected-num_align_aa) == 0 and (num_aligns_expected-num_align_dna) == 0: + logging.debug("Number of aligns expected is " + str(num_aligns_expected) + " and number in " + str(self._folder_align_append_aa) + " is " + str(num_align_aa) + ", similarly for aa. So this step is done") return True else: + logging.debug("Number of aligns expected is " + str(num_aligns_expected) + " and number in " + str(self._folder_align_append_aa) + " and aa versions are " + str(num_align_aa) + " and " + str(num_align_dna) + ". So this step is not done") return False else: return False @@ -233,6 +249,7 @@ def _get_finished_mapping_folders(self, path): mapping_folders_finished = [] num_expected_mappings = self._get_number_of_references() mapping_folders = [x for x in os.listdir(path) if '04' in x] + logging.debug("Number of mapping expected is " + str(num_expected_mappings) + " we are checking folders in " + str(path)) for folder in mapping_folders: # NOTE: we are calculating the number of completed mappings as the number of existing cov files, # because these are written even if the mapping step did not find any reads to map to a particular reference @@ -258,10 +275,13 @@ def _get_mapping_status(self): if len(mapping_folders) > 0: self.num_completed_mappings = len(mapping_folders) # self.logger.info('{}: Mapping completed!'.format(self._species_name)) + logging.debug("There are some mapping folders " + str(self.num_completed_mappings)) return True else: self.num_completed_mappings = 0 + logging.debug("There is no mapping folder") # self.logger.info('{}: Mapping not completed!'.format(self._species_name)) return False else: + logging.debug("There is no mapping folder") return False From b0120c43c90dea41d36867cc0256487aa6db59b9 Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 12:36:38 +0200 Subject: [PATCH 17/50] add todo --- read2tree/main.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/read2tree/main.py b/read2tree/main.py index f10c32b..66364f0 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -248,7 +248,7 @@ def parse_args(argv, exe_name, desc): 'Arguments --coverage and --genome_len ' 'can only be set if --sample_reads is set.') - progress = Progress(args) + progress = Progress(args) # todo why calling Progress twice? if progress.num_completed_mappings <= 1 and args.merge_all_mappings: arg_parser.error('The number of completed mappings ({}) is too ' 'little to perform a merge.'.format(progress.num_completed_mappings)) @@ -266,10 +266,10 @@ def main(argv, exe_name, desc=''): args = parse_args(argv, exe_name, desc) logger.info('{}: ------- NEW RUN -------'.format(args.species_name)) - x = ', '.join("{!s}={!r}".format(key, val) for (key, val) in vars(args).items()) + x = ', '.join("{!s}={!r}".format(key, val) for (key, val) in vars(args).items()) # todo why calling Progress twice? logger.info('{}: read2tree was run with: {}'.format(args.species_name, x)) - progress = Progress(args) + progress = Progress(args) # todo why calling Progress twice? if not os.path.exists(args.output_path): os.makedirs(args.output_path) logger.info('{}: Progress: ogs_dna {} | ref {} | ref_align {} | mapping {} | append_ogs {} | align {} ' From a131e4675a287a1848b27ada1f88c507211aba5d Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 14:44:06 +0200 Subject: [PATCH 18/50] improve log --- read2tree/Mapper.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 3608a4b..c9e8dae 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -353,8 +353,9 @@ def _map_reads_to_references(self, ref): .format(self._species_name, self.elapsed_time)) if not mapped_reads_species : - self.logger.warning('mapped_reads_species is empty') - + log_line='mapped_reads_species is empty! are the aligner and samtools installed correctly? are there enough marker genes close to the species of interest? Please remove the output folder and mplog file and re-run. ' + self.logger.warning(log_line) + print(log_line) return mapped_reads_species From bc2e0fca1aa7aea92d05da2126cbe8c27f15070b Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 14:44:17 +0200 Subject: [PATCH 19/50] print version --- read2tree/main.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/read2tree/main.py b/read2tree/main.py index 66364f0..a9805e5 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -260,6 +260,10 @@ def main(argv, exe_name, desc=''): ''' Main function. ''' + from . import __version__ as r2t_version + + logger.debug(' ------- Read2Tree version: {} -------'.format(r2t_version)) + print(' ------- Read2Tree version: '+str(r2t_version)+' -------') t1 = timer() # Parse From 503728fa0cc80d347fd83b34915813a48ee091a0 Mon Sep 17 00:00:00 2001 From: sina Date: Sun, 31 Mar 2024 14:48:03 +0200 Subject: [PATCH 20/50] enable logger --- read2tree/main.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/read2tree/main.py b/read2tree/main.py index a9805e5..6c47466 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -30,7 +30,7 @@ COPYRIGHT = '(C) 2017-{:d} David V Dylus'.format(date.today().year) logger = logging.getLogger(__name__) -logger.disabled = True +#logger.disabled = True def parse_args(argv, exe_name, desc): ''' @@ -262,7 +262,7 @@ def main(argv, exe_name, desc=''): ''' from . import __version__ as r2t_version - logger.debug(' ------- Read2Tree version: {} -------'.format(r2t_version)) + logger.info(' ------- Read2Tree version: {} -------'.format(r2t_version)) print(' ------- Read2Tree version: '+str(r2t_version)+' -------') t1 = timer() From 4b9fd4e36ffc054e1cf8387218c9feb806db86c5 Mon Sep 17 00:00:00 2001 From: sina Date: Wed, 17 Apr 2024 14:30:31 -0400 Subject: [PATCH 21/50] fix a bug for paired end reads --- read2tree/Mapper.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index c9e8dae..960359d 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -110,7 +110,12 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): minimap_argm = " -ax map-hifi " elif 'long-ont' in self.args.read_type: minimap_argm = " -ax map-ont " - line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads + " > " + bam_file + if isinstance(reads,list): + reads_str=" ".join(reads) + elif isinstance(reads,str): + reads_str = reads + + line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads_str + " > " + bam_file self._output_shell(line_minimp) self.logger.info('mapping with' + line_minimp) From fcdc8e3afb2eeb86a1d55ace004e8635bda4ac93 Mon Sep 17 00:00:00 2001 From: sina Date: Mon, 29 Apr 2024 14:29:43 -0400 Subject: [PATCH 22/50] fix a bug multi-sample --- read2tree/Mapper.py | 197 +++++++++++++++++++++++--------------------- 1 file changed, 104 insertions(+), 93 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 960359d..f6d9616 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -25,6 +25,7 @@ import functools from Bio import SeqIO, SeqRecord, Seq +import sys try: from Bio.Alphabet import generic_dna except ImportError: @@ -95,13 +96,21 @@ def __init__(self, args, ref_set=None, og_set=None, species_name=None, progress= self._read_mapping_from_folder(mapping_name=self._mapping_name, ref_records=ref_set) self.og_records = self._sort_by_og() + # if not (self.mapped_records or hasattr(self,"og_records")): + # log_line='mapped_records or og_records is empty! are the aligner and samtools installed correctly? are there enough marker genes close to the species of interest? If it is a new run, make sure there is no output folder or a file mplog.log' + # self.logger.error(log_line) + # print(log_line) + # sys.exit(-1) + + + def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): start = time.time() - output_folder = os.path.join(self.args.output_path, - "04_mapping_" + self._species_name) - - bam_file= tmp_output_folder+"/"+ref_file_handle.split("/")[-1]+".bam" + sample_reads = str(self._species_name) # input sequencing read + output_folder = os.path.join(self.args.output_path, "04_mapping_" + sample_reads) + # ref_file_handle 'output/02_ref_dna/PYGNA_OGs.fa' + sam_file= output_folder+"/"+ref_file_handle.split('/')[-1].split('.')[0]+".sam" #self._output_shell(minimap2_ex+" -ax sr "+ ref_file_handle+ " -t " + str(self.args.threads) + " "+reads+" | "+samtools_ex+" view -F 4 -bh -S -t" + str(self.args.threads)+ " > " + bam_file) if 'short' in self.args.read_type: @@ -110,17 +119,16 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): minimap_argm = " -ax map-hifi " elif 'long-ont' in self.args.read_type: minimap_argm = " -ax map-ont " + else: + self.logger.error(" read_type is not valid" +self.args.read_type+" please use one of these: short long-ont long-hifi ") if isinstance(reads,list): reads_str=" ".join(reads) elif isinstance(reads,str): reads_str = reads - line_minimp= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads_str + " > " + bam_file - self._output_shell(line_minimp) - - self.logger.info('mapping with' + line_minimp) - - + line_minimap= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads_str + " > " + sam_file + self._output_shell(line_minimap) + self.logger.info('mapping with ' + line_minimap) # if len(self._reads) == 2: # ngm_wrapper = NGM(ref_file_handle, reads, tmp_output_folder) @@ -152,15 +160,13 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): # self._rm_file(ref_file_handle + "-ht-13-2.2.ngm", ignore_error=True) # self._rm_file(ref_file_handle + "-ht-13-2.3.ngm", ignore_error=True) # - # end = time.time() - # self.elapsed_time = end - start - # self.logger.info('{}: Mapping to {} references took {}.' - # .format(self._species_name, os.path.basename(ref_file_handle), - # self.elapsed_time)) - - if os.path.exists(bam_file) and os.path.getsize(bam_file) > 0: # ngm['reads_mapped'] > 0 and - shutil.copy(bam_file, os.path.join(output_folder, os.path.basename(bam_file))) - return self._post_process_read_mapping(ref_file_handle, bam_file) + end = time.time() + self.elapsed_time = end - start + self.logger.info('{}: Mapping to {} references took {}.' .format(self._species_name, os.path.basename(ref_file_handle),self.elapsed_time)) + + if os.path.exists(sam_file) and os.path.getsize(sam_file) > 2: # ngm['reads_mapped'] > 0 and + #shutil.copy(sam_file, os.path.join(output_folder, os.path.basename(sam_file))) + return self._post_process_read_mapping(ref_file_handle, sam_file) else: open(os.path.join(output_folder, os.path.basename(ref_file_handle).split('.')[0]+'_cov.txt'), 'a').close() return None @@ -182,12 +188,12 @@ def _read_mapping_from_folder(self, mapping_name=None, ref_records=None): if self.args.min_cons_coverage >= 2 and bam_files: for file in tqdm(bam_files, desc='Generating consensus from bam files ', unit=' species'): species = file.split("/")[-1].split("_")[0] - ref_file = os.path.join(self.args.output_path, '02_ref_dna', - species+'_OGs.fa') + ref_file = os.path.join(self.args.output_path, '02_ref_dna',species+'_OGs.fa') map_reads_species[species] = Reference() - self._output_shell( - samtools+' index -@ ' + str(self.args.threads) + ' ' + - file) + line_samtools_index = samtools+' index -@ ' + str(self.args.threads) + ' ' +file + self._output_shell(line_samtools_index) + self.logger.debug(line_samtools_index) + consensus = self._build_consensus_seq_v2(ref_file, file) records = [] @@ -348,7 +354,7 @@ def _map_reads_to_references(self, ref): mapped_reads = [] # self.progress.set_status('single_map', ref=species) - self._rm_file(ref_file_handle+".fai", ignore_error=True) + #self._rm_file(ref_file_handle+".fai", ignore_error=True) #tmp_output_folder.cleanup() end = time.time() @@ -358,7 +364,7 @@ def _map_reads_to_references(self, ref): .format(self._species_name, self.elapsed_time)) if not mapped_reads_species : - log_line='mapped_reads_species is empty! are the aligner and samtools installed correctly? are there enough marker genes close to the species of interest? Please remove the output folder and mplog file and re-run. ' + log_line='mapped_reads_species is empty! are the aligner and samtools installed correctly? are there enough marker genes close to the species of interest? If it is a new run, make sure there is no output folder or a file mplog.log' self.logger.warning(log_line) print(log_line) @@ -446,29 +452,31 @@ def _get_mapping_stats(self, bam_file): for l in x]) return mapped, all_reads - def _bin_reads(self, ref_file, bam_file): - """ - Function that bins reads into their orthologous groups - :param ref_file: Current species reference file - :param bam_file: Mapped bam file - """ - self.logger.debug("{}: --- Binning reads ---".format(self._species_name)) - output_folder = os.path.join(self.args.output_path, "04_read_ogs_" + - self._species_name) - if not os.path.exists(output_folder): - os.makedirs(output_folder) - tmp_folder = os.path.dirname(bam_file) - outfile_name = os.path.join(tmp_folder, - ref_file.split('/')[-1].split('.')[0] + - "_post") - - shutil.copy(bam_file, os.path.join(output_folder, - os.path.basename(outfile_name + - "_sorted.bam"))) - shutil.copy(bam_file + ".bai", - os.path.join(output_folder, - os.path.basename(outfile_name + - "_sorted.bam.bai"))) + # def _bin_reads(self, ref_file, sam_file_base): + # """ + # Function that bins reads into their orthologous groups + # :param ref_file: Current species reference file + # :param bam_file: Mapped bam file + # """ + # bam_file = sam_file_base+ '_sorted.bam' + # + # self.logger.debug("{}: --- Binning reads ---".format(self._species_name)) + # output_folder = os.path.join(self.args.output_path, "04_read_ogs_" + + # self._species_name) + # if not os.path.exists(output_folder): + # os.makedirs(output_folder) + # tmp_folder = os.path.dirname(bam_file) + # outfile_name = os.path.join(tmp_folder, + # ref_file.split('/')[-1].split('.')[0] + + # "_post") + # + # shutil.copy(bam_file, os.path.join(output_folder, + # os.path.basename(outfile_name + + # "_sorted.bam"))) + # shutil.copy(bam_file + ".bai", + # os.path.join(output_folder, + # os.path.basename(outfile_name + + # "_sorted.bam.bai"))) # if os.path.exists(bam_file): # bam = pysam.AlignmentFile(bam_file, "rb") @@ -548,7 +556,7 @@ def _build_consensus_seq_v2(self, ref_file, bam_file): new_records[ref] = ("").join(seq) return new_records - def _post_process_read_mapping(self, ref_file, bam_file): + def _post_process_read_mapping(self, ref_file, sam_file): """ Function that will perform postprocessing of finished read mapping using the pysam functionality @@ -557,50 +565,53 @@ def _post_process_read_mapping(self, ref_file, bam_file): :return: """ # self.logger.info("--- Postprocessing reads to {} ---".format(self._species_name)) - output_folder = os.path.join(self.args.output_path, - "04_mapping_"+self._species_name) - tmp_folder = os.path.dirname(bam_file) - outfile_name = os.path.join(tmp_folder, - ref_file.split('/')[-1].split('.')[0] + - "_post") + output_folder = os.path.join(self.args.output_path, "04_mapping_"+self._species_name) + #tmp_folder = os.path.dirname(bam_file) + #outfile_name = os.path.join(tmp_folder, + # ref_file.split('/')[-1].split('.')[0] + + # "_post") if self.args.single_mapping: - self.logger.debug("{}: --- POSTPROCESSING MAPPING " - "---".format(self._species_name)) + self.logger.debug("{}: --- POSTPROCESSING MAPPING ---".format(self._species_name)) # ngmlr doesn't have the option to write in bam file directly - if 'sam' in bam_file.split(".")[-1]: - sam_file = bam_file - bam_file = sam_file.replace(".sam", ".bam") - if os.path.exists(sam_file): - self._output_shell( - samtools+' view -F 4 -bh -S -@ ' + str(self.args.threads) + - ' -o ' + bam_file + " " + sam_file) + #if 'sam' in bam_file.split(".")[-1]: + # sam_file = bam_file + sam_file_base = sam_file[:-4] + if os.path.exists(sam_file): + line_samtools_view = samtools+' view -F 4 -bh -S -@ ' + str(self.args.threads) +' -o ' + sam_file_base + ".bam " + sam_file + self._output_shell(line_samtools_view) + self.logger.debug(line_samtools_view) + else: + self.logger.debug("Sam file is not generated", sam_file) + if self.args.single_mapping: - self.logger.debug("{}: ---- Samtools view completed" - .format(self._species_name)) + self.logger.error("single mapping is not tested with this version ") + #self.logger.debug("{}: ---- Samtools view completed".format(self._species_name)) + + if os.path.exists(sam_file_base+".bam"): + + line_samtools_sort = samtools+' sort -@ ' + str(self.args.threads) + ' -o ' + sam_file_base + "_sorted.bam " + sam_file_base+".bam" + self._output_shell(line_samtools_sort) + self.logger.debug("running "+ line_samtools_sort) + else: + self.logger.debug("bam file is not generated", sam_file_base+".bam") - if os.path.exists(bam_file): - self._output_shell( - samtools+' sort -@ ' + str(self.args.threads) + # todo 2G might not be enough - ' -o ' + outfile_name + "_sorted.bam " + bam_file) if self.args.single_mapping: - self.logger.debug("{}: ---- Samtools sort completed" - .format(self._species_name)) + self.logger.debug("{}: ---- Samtools sort completed".format(self._species_name)) + + if os.path.exists(sam_file_base + "_sorted.bam"): + line_samtools_index= samtools+' index -@ ' + str(self.args.threads) + ' ' + sam_file_base + "_sorted.bam" + self._output_shell(line_samtools_index) + self.logger.debug(line_samtools_index) - if os.path.exists(outfile_name + "_sorted.bam"): - self._output_shell( - samtools+' index -@ ' + str(self.args.threads) + ' ' + - outfile_name + "_sorted.bam") if self.args.single_mapping: - self.logger.debug("{}: ---- Samtools index completed" - .format(self._species_name)) + self.logger.warning("single_mapping is not tested in this version.") #debug("{}: ---- Samtools index completed".format(self._species_name)) # self._rm_file(bam_file, ignore_error=True) - if self.args.debug: - self._bin_reads(ref_file, outfile_name + '_sorted.bam') + #f self.args.debug: + # self._bin_reads(ref_file, sam_file_base ) - consensus = self._build_consensus_seq_v2(ref_file, outfile_name + - '_sorted.bam') + consensus = self._build_consensus_seq_v2(ref_file, sam_file_base +'_sorted.bam') all_consensus = [] if consensus: @@ -628,7 +639,7 @@ def _post_process_read_mapping(self, ref_file, bam_file): # Get effective coverage of each mapped sequence cov = Coverage(self.args) - cov.get_coverage_bam(outfile_name + "_sorted.bam") + cov.get_coverage_bam(sam_file_base + "_sorted.bam") cov.write_coverage_bam(os.path.join( output_folder, ref_file.split('/')[-1].split('.')[0] + "_cov.txt")) self.all_cov.update(cov.coverage) @@ -643,15 +654,15 @@ def _rm_file(self, *fns, ignore_error=False): if not ignore_error: raise - def _clean_up_tmp_files_single(self, species): - output_folder = os.path.join(self.args.output_path, "04_mapping_" + - self._species_name) - fn_ends = ('_post.bam', '_post_consensus_call.fq', '_post_sorted.bam', - '_post_sorted.bam.bai', '.fa.fai', - '.fa.sam', '.fa-ht-13-2.3.ngm', '.fa-ht-13-2.3.ngm', - '.fa', '.fa-enc.2.ngm') - self._rm_file(*[os.path.join(output_folder, species + fn_end) - for fn_end in fn_ends], ignore_error=True) + # def _clean_up_tmp_files_single(self, species): + # output_folder = os.path.join(self.args.output_path, "04_mapping_" + + # self._species_name) + # fn_ends = ('_post.bam', '_post_consensus_call.fq', '_post_sorted.bam', + # '_post_sorted.bam.bai', '.fa.fai', + # '.fa.sam', '.fa-ht-13-2.3.ngm', '.fa-ht-13-2.3.ngm', + # '.fa', '.fa-enc.2.ngm') + # self._rm_file(*[os.path.join(output_folder, species + fn_end) + # for fn_end in fn_ends], ignore_error=True) def _output_shell(self, line): """ From bc3cd2ad48dc831ec39fd0d12099bf0a9f5fb67d Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 10:31:47 -0400 Subject: [PATCH 23/50] update --- read2tree/__init__.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/read2tree/__init__.py b/read2tree/__init__.py index 9063fd7..1d8fe87 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -6,7 +6,7 @@ from pkg_resources import resource_string logging.getLogger(__name__).addHandler(logging.NullHandler()) -__version__ = '0.2.0' +__version__ = '1.5.0' __copyright__ = 'read2tree (C) 2017-{:d} David Dylus' \ .format(date.today().year) @@ -16,9 +16,9 @@ # config = yaml.load(f.read()) # logging.config.dictConfig(config) -conf = resource_string(__name__, 'logging/log.yaml') +#conf = resource_string(__name__, 'logging/log.yaml') -D = yaml.load(conf, Loader=yaml.FullLoader) -D.setdefault('version', 1) -logging.config.dictConfig(D) +# D = yaml.load(conf, Loader=yaml.FullLoader) +# D.setdefault('version', 1) +# logging.config.dictConfig(D) # del D From 174d04668f84b7de69aa6b7b18ac4f9deda598c9 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 10:33:37 -0400 Subject: [PATCH 24/50] remove progress, define steps --- read2tree/main.py | 432 +++++++++++++++++++++++----------------------- 1 file changed, 213 insertions(+), 219 deletions(-) diff --git a/read2tree/main.py b/read2tree/main.py index 6c47466..551af7a 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -1,4 +1,5 @@ #!/usr/bin/env python + ''' read2tree: from reads to trees using the OMA standalone Output and a set of reads 1) DNA sequences are retrieved from OMA standalone Output based on the predicted @@ -22,15 +23,25 @@ from read2tree.ReferenceSet import ReferenceSet from read2tree.Mapper import Mapper from read2tree.Aligner import Aligner -from read2tree.Progress import Progress +# from read2tree.Progress import Progress from read2tree.TreeInference import TreeInference from read2tree.parser import OMAOutputParser import argparse +import glob + +import sys COPYRIGHT = '(C) 2017-{:d} David V Dylus'.format(date.today().year) +# logger = logging.getLogger(__name__) +logger_level = "DEBUG" # DEBUG INFO # TRACE DEBUG INFO WARN ERROR FATAL +logging.basicConfig(format='%(asctime)s %(levelname)-8s %(message)s', level=logging.INFO, datefmt='%Y-%m-%d %H:%M:%S') logger = logging.getLogger(__name__) -#logger.disabled = True +if logger_level == "INFO": + logger.setLevel(logging.INFO) + + +# logger.disabled = True def parse_args(argv, exe_name, desc): ''' @@ -46,52 +57,51 @@ def parse_args(argv, exe_name, desc): # if not is_standalone: # # If standalone, set in parser. - arg_parser.add_argument('--version', action='version', help='Show programme\'s version number and exit.', version=read2tree.__version__) arg_parser.add_argument('--output_path', default='.', help='[Default is current directory] Path to ' - 'output directory.') + 'output directory.') - arg_parser.add_argument('--standalone_path', default='.', required=True, + arg_parser.add_argument('--standalone_path', default='.', required=True, # todo change name to marker_gene help='[Default is current directory] Path to ' 'oma standalone directory.') arg_parser.add_argument('--reads', nargs='+', default=None, help='[Default is none] Reads to be mapped to reference. If paired ' - 'end add separated by space.') + 'end add separated by space.') arg_parser.add_argument('--read_type', default='short', help='[Default is short] Type of reads to ' - 'use for mapping: short, long-hifi or long-ont corresponding to sr, map-hifi, or map-ont in minimap2.') + 'use for mapping: short, long-hifi or long-ont corresponding to sr, map-hifi, or map-ont in minimap2.') arg_parser.add_argument('--threads', type=int, default=1, help='[Default is 1] Number of threads for the mapping ') - arg_parser.add_argument('--split_reads', action='store_true', - help='[Default is off] Splits reads as defined by split_len (200) ' - 'and split_overlap (0) parameters. ') + # arg_parser.add_argument('--split_reads', action='store_true', + # help='[Default is off] Splits reads as defined by split_len (200) ' + # 'and split_overlap (0) parameters. ') - arg_parser.add_argument('--split_len', type=int, default=200, - help='[Default is 200] Parameter for selection of ' - 'read split length can only be used in combination' - 'with with long read option. ') + # arg_parser.add_argument('--split_len', type=int, default=200, + # help='[Default is 200] Parameter for selection of ' + # 'read split length can only be used in combination' + # 'with with long read option. ') + # + # arg_parser.add_argument('--split_overlap', type=int, default=0, + # help='[Default is 0] Reads are split with an ' + # 'overlap defined by this argument.') - arg_parser.add_argument('--split_overlap', type=int, default=0, - help='[Default is 0] Reads are split with an ' - 'overlap defined by this argument.') - - arg_parser.add_argument('--split_min_read_len', type=int, default=200, - help='[Default is 200] Reads longer than this ' - 'value are cut into smaller values as defined ' - 'by --split_len. ') + # arg_parser.add_argument('--split_min_read_len', type=int, default=200, + # help='[Default is 200] Reads longer than this ' + # 'value are cut into smaller values as defined ' + # 'by --split_len. ') arg_parser.add_argument('--sample_reads', action='store_true', help='[Default is off] Splits reads as defined by split_len (200) ' - 'and split_overlap (0) parameters. ') - + 'and split_overlap (0) parameters. ') + arg_parser.add_argument('--genome_len', type=int, default=2000000, help='[Default is 2000000] Genome size in bp.') @@ -101,19 +111,19 @@ def parse_args(argv, exe_name, desc): arg_parser.add_argument('--min_cons_coverage', type=int, default=1, help='[Default is 1] Minimum number of nucleotides at column.') - arg_parser.add_argument('--dna_reference', default='', + arg_parser.add_argument('--dna_reference', default='', # todo make it mandatory dna_reference no api help='[Default is None] Reference file that contains nucleotide ' - 'sequences (fasta, hdf5). If not given it will use' - 'the RESTapi and retrieve sequences ' - 'from http://omabrowser.org directly. ' - 'NOTE: internet connection required!') + 'sequences (fasta, hdf5). If not given it will use' + 'the RESTapi and retrieve sequences ' + 'from http://omabrowser.org directly. ' + 'NOTE: internet connection required!') arg_parser.add_argument('--sc_threshold', type=float, default=0.25, help='[Default is 0.25; Range 0-1] Parameter for ' - 'selection of sequences from mapping by ' - 'completeness compared to its reference sequence ' - '(number of ACGT basepairs vs length of sequence). ' - 'By default, all sequences are selected.') + 'selection of sequences from mapping by ' + 'completeness compared to its reference sequence ' + '(number of ACGT basepairs vs length of sequence). ' + 'By default, all sequences are selected.') # arg_parser.add_argument('--ngmlr_parameters', default=None, # todo this could be used for minimap2 options # help='[Default is none] In case this parameters ' @@ -123,15 +133,14 @@ def parse_args(argv, exe_name, desc): # 'parameter can be found in the original ' # 'documentation of ngmlr.') + # arg_parser.add_argument('--check_mate_pairing', action='store_true', + # help='Check whether in case of paired end ' + # 'reads we have consistent mate pairing. Setting ' + # 'this option will automatically select the ' + # 'overlapping reads and do not consider single ' + # 'reads.') - arg_parser.add_argument('--check_mate_pairing', action='store_true', - help='Check whether in case of paired end ' - 'reads we have consistent mate pairing. Setting ' - 'this option will automatically select the ' - 'overlapping reads and do not consider single ' - 'reads.') - - arg_parser.add_argument('--debug', action='store_true', + arg_parser.add_argument('--debug', action='store_true', # todo make it active always otherwise change it back help='[Default is false] Changes to debug mode: ' '* bam files are saved!' '* reads are saved by mapping to OG') @@ -148,63 +157,62 @@ def parse_args(argv, exe_name, desc): help='[Default is false] Compute tree, otherwise just ' 'output concatenated alignment!') - arg_parser.add_argument('--merge_all_mappings', action='store_true', - help='[Default is off] In case multiple species were mapped to ' - 'the same reference this allows to merge this ' - 'mappings and build a tree with all included ' - 'species!') + arg_parser.add_argument('--step', default="all", + help='[Default is all] marker map combine ') + + # arg_parser.add_argument('--merge_all_mappings', action='store_true', + # help='[Default is off] In case multiple species were mapped to ' + # 'the same reference this allows to merge this ' + # 'mappings and build a tree with all included ' + # 'species!') # Arguments to generate the reference - arg_parser.add_argument('-r', '--reference', action='store_true', - help='[Default is off] Just generate the reference dataset for ' - 'mapping.') + # arg_parser.add_argument('-r', '--reference', action='store_true', + # help='[Default is off] Just generate the reference dataset for ' + # 'mapping.') arg_parser.add_argument('--min_species', type=int, default=None, help='Min number of species in selected ' - 'orthologous groups. If not selected it will be ' - 'estimated such that around 1000 OGs ' - 'are available.') + 'orthologous groups. If not selected it will be ' + 'estimated such that around 1000 OGs ' + 'are available.') - arg_parser.add_argument('--single_mapping', default=None, - help='[Default is none] Single species file allowing to map in a ' - 'job array.') + # arg_parser.add_argument('--single_mapping', default=None, + # help='[Default is none] Single species file allowing to map in a ' + # 'job array.') # Arguments to map the reads - arg_parser.add_argument('--ref_folder', default=None, - help='[Default is none] Folder containing reference files with ' - 'sequences sorted by species.') - + # arg_parser.add_argument('--ref_folder', default=None, + # help='[Default is none] Folder containing reference files with ' + # 'sequences sorted by species.') arg_parser.add_argument('--remove_species_mapping', default=None, help='[Default is none] Remove species present in ' - 'data set after mapping step completed and only ' - 'do analysis on subset. Input is comma separated ' - 'list without spaces, e.g. XXX,YYY,AAA.') + 'data set after mapping step completed and only ' + 'do analysis on subset. Input is comma separated ' + 'list without spaces, e.g. XXX,YYY,AAA.') arg_parser.add_argument('--remove_species_ogs', default=None, help='[Default is none] Remove species present ' - 'in data set after mapping step completed to ' - 'build OGs. Input is comma separated list ' - 'without spaces, e.g. XXX,YYY,AAA.') + 'in data set after mapping step completed to ' + 'build OGs. Input is comma separated list ' + 'without spaces, e.g. XXX,YYY,AAA.') arg_parser.add_argument('--keep_all_ogs', action='store_true', default=True, help='[Default is on] Keep all orthologs after addition of ' - 'mapped seq, which means also the OGs that ' - 'have no mapped sequence. Otherwise only OGs ' - 'are used that have the mapped sequence for ' - 'alignment and tree inference.') + 'mapped seq, which means also the OGs that ' + 'have no mapped sequence. Otherwise only OGs ' + 'are used that have the mapped sequence for ' + 'alignment and tree inference.') arg_parser.add_argument('--ignore_species', default=None, help='[Default is none] Ignores species part of ' - 'the OMA standalone pipeline. Input is comma ' - 'separated list without spaces, e.g. XXX,YYY,AAA.') + 'the OMA standalone pipeline. Input is comma ' + 'separated list without spaces, e.g. XXX,YYY,AAA.') # Parse the arguments. args = arg_parser.parse_args(argv) - if not os.path.exists(args.output_path): - os.makedirs(args.output_path) - _reads = "" _species_name = "" @@ -220,22 +228,22 @@ def parse_args(argv, exe_name, desc): if args.species_name: _species_name = args.species_name - if args.merge_all_mappings: + if args.step == "combine": # todo why is needed? _species_name = 'merge' args.reads = _reads args.species_name = _species_name - if not args.split_reads and (args.split_len != 200 or - args.split_overlap != 0 or - args.split_min_read_len != 200): - arg_parser.error( - 'Arguments --split_len, --split_overlap and --split_min_read_len' - 'can only be set if --split_reads is set.') + # if not args.split_reads and (args.split_len != 200 or + # args.split_overlap != 0 or + # args.split_min_read_len != 200): + # arg_parser.error( + # 'Arguments --split_len, --split_overlap and --split_min_read_len' + # 'can only be set if --split_reads is set.') - if args.split_reads and len(args.reads) == 2: - arg_parser.error( - 'Splitting reads does not work for paired end reads.') + # if args.split_reads and len(args.reads) == 2: + # arg_parser.error( + # 'Splitting reads does not work for paired end reads.') # if args.read_type == 'short' and args.ngmlr_parameters: # arg_parser.error( @@ -248,14 +256,38 @@ def parse_args(argv, exe_name, desc): 'Arguments --coverage and --genome_len ' 'can only be set if --sample_reads is set.') - progress = Progress(args) # todo why calling Progress twice? - if progress.num_completed_mappings <= 1 and args.merge_all_mappings: - arg_parser.error('The number of completed mappings ({}) is too ' - 'little to perform a merge.'.format(progress.num_completed_mappings)) + # progress = Progress(args) # todo why calling Progress twice? + # if progress.num_completed_mappings <= 1 and args.merge_all_mappings: + # arg_parser.error('The number of completed mappings ({}) is too ' + # 'little to perform a merge.'.format(progress.num_completed_mappings)) return args +def get_finished_mapping_folders2(output_path): + mapping_folders_finished = [] + # num_expected_mappings = self._get_number_of_references() + mapping_folders = [x for x in os.listdir(output_path) if '04' in x] + # logging.debug("Number of mapping expected is " + str(num_expected_mappings) + " we are checking folders in " + str(path)) + for folder in mapping_folders: + # NOTE: we are calculating the number of completed mappings as the number of existing cov files, + # because these are written even if the mapping step did not find any reads to map to a particular reference + computed_cov = [f for f in + glob.glob(os.path.join(output_path, + folder + '/*cov.txt'))] + # it is finished if the number of generated coverage files is the same as the number of references + # if self.args.merge_all_mappings: + # if num_expected_mappings >= len(computed_cov) and len(computed_cov) >= 1: + # mapping_folders_finished.append(folder) + # elif self._species_name in folder: + # if (num_expected_mappings - len(computed_cov)) == 0: + # mapping_folders_finished.append(folder) + if len(computed_cov) >= 1: + mapping_folders_finished.append(folder) # todo needs testing + + return mapping_folders_finished + + def main(argv, exe_name, desc=''): ''' Main function. @@ -263,7 +295,7 @@ def main(argv, exe_name, desc=''): from . import __version__ as r2t_version logger.info(' ------- Read2Tree version: {} -------'.format(r2t_version)) - print(' ------- Read2Tree version: '+str(r2t_version)+' -------') + # print(' ------- Read2Tree version: '+str(r2t_version)+' -------') t1 = timer() # Parse @@ -273,161 +305,123 @@ def main(argv, exe_name, desc=''): x = ', '.join("{!s}={!r}".format(key, val) for (key, val) in vars(args).items()) # todo why calling Progress twice? logger.info('{}: read2tree was run with: {}'.format(args.species_name, x)) - progress = Progress(args) # todo why calling Progress twice? - if not os.path.exists(args.output_path): - os.makedirs(args.output_path) - logger.info('{}: Progress: ogs_dna {} | ref {} | ref_align {} | mapping {} | append_ogs {} | align {} ' - .format(args.species_name, progress.ref_ogs_01, progress.ref_dna_02, - progress.ref_align_03, progress.mapping_04, progress.append_ogs_05, progress.align_06)) + logger.info("Running r2t in mode " + args.step) - # TODO: Check whether all the necessary binaries are available - # TODO: Check all given files and throw error if faulty + if os.path.exists(args.output_path): + if args.step == "all" or args.step == "marker": + logger.error( + "the output folder exist" + args.output_path + ". Since you are running r2t in all mode, you need to specify output folder which will created by r2t.") + sys.exit() + else: + os.makedirs(args.output_path) - if (not progress.ref_ogs_01 and not progress.ref_dna_02 and - not progress.ref_align_03 and not progress.mapping_04 and - not progress.append_ogs_05 and not progress.align_06): + if args.step == "all" or args.step == "map": + if not args.reads: + logger.error("reads are not provided in mode " + str(args.step)) + sys.exit() + if isinstance(args.reads, + list): # in parse_args it is input read is converted to a list, if there are 2 - paired end + for read_file in args.reads: + if not os.path.isfile(read_file): + logger.error("read file doesn't exist" + read_file) + sys.exit() + elif isinstance(args.reads, str): + if not os.path.isfile(args.reads): + logger.error("read file doesn't exist") + sys.exit() + + if args.step == "all": oma_output = OMAOutputParser(args) args.oma_output_path = oma_output.oma_output_path - ogset = OGSet(args, oma_output=oma_output, progress=progress) # Generate the OGs with their DNA sequences - reference = ReferenceSet(args, og_set=ogset.ogs, load=True, progress=progress) - alignments = Aligner(args, ogset.ogs, load=True) - if not args.reference: - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, progress=progress) - alignments.remove_species_from_alignments() - ogset.remove_species_from_ogs() - ogset.add_mapped_seq(mapper) - ogset.write_added_ogs_aa() - ogset.write_added_ogs_dna() - # alignments = Aligner(args, ogset.mapped_ogs, load=True) - alignments.add_mapped_seq(ogset.mapped_ogs) - alignments.write_added_align_aa() - alignments.write_added_align_dna() - concat_alignment = alignments.concat_alignment() - if args.tree: - tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (progress.ref_ogs_01 and not progress.ref_dna_02 and - not progress.ref_align_03 and not progress.mapping_04 and - not progress.append_ogs_05 and not progress.align_06): - ogset = OGSet(args, load=False, progress=progress) - reference = ReferenceSet(args, og_set=ogset.ogs, load=True, progress=progress) # Generate the reference - alignments = Aligner(args, ogset.ogs, load=True) - if not args.reference: # just generate reference - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, progress=progress) - alignments.remove_species_from_alignments() - ogset.remove_species_from_ogs() - ogset.add_mapped_seq(mapper) - ogset.write_added_ogs_aa() - ogset.write_added_ogs_dna() - alignments.add_mapped_seq(ogset.mapped_ogs) - alignments.write_added_align_aa() - alignments.write_added_align_dna() - concat_alignment = alignments.concat_alignment() - if args.tree: - tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (progress.ref_ogs_01 and progress.ref_dna_02 and - not progress.ref_align_03 and not progress.mapping_04 and - not progress.append_ogs_05 and not progress.align_06): - ogset = OGSet(args, load=False, progress=progress) - reference = ReferenceSet(args, load=False, progress=progress) - alignments = Aligner(args, ogset.ogs, load=True) - if not args.reference: # just generate reference - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, progress=progress) - alignments.remove_species_from_alignments() - ogset.remove_species_from_ogs() - ogset.add_mapped_seq(mapper) - ogset.write_added_ogs_aa() - ogset.write_added_ogs_dna() - alignments.add_mapped_seq(ogset.mapped_ogs) - alignments.write_added_align_aa() - alignments.write_added_align_dna() - concat_alignment = alignments.concat_alignment() - if args.tree: - tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (progress.ref_ogs_01 and progress.ref_dna_02 and - progress.ref_align_03 and not progress.mapping_04 and - not progress.append_ogs_05 and not progress.align_06): - if args.single_mapping: - reference = ReferenceSet(args, load=False, progress=progress) - Mapper(args, ref_set=reference.ref) # Run the mapping - else: - ogset = OGSet(args, load=False, progress=progress) - reference = ReferenceSet(args, load=False, progress=progress) - alignments = Aligner(args, load=False) - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, progress=progress) # Run the mapping - alignments.remove_species_from_alignments() - ogset.remove_species_from_ogs() - ogset.add_mapped_seq(mapper) - ogset.write_added_ogs_aa() - ogset.write_added_ogs_dna() - alignments.add_mapped_seq(ogset.mapped_ogs) - alignments.write_added_align_aa() - alignments.write_added_align_dna() - concat_alignment = alignments.concat_alignment() - if args.tree: - tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (args.merge_all_mappings and progress.ref_ogs_01 and progress.ref_dna_02 and - progress.ref_align_03 and progress.mapping_04 and - not progress.append_ogs_05 and not progress.align_06): - ogset = OGSet(args, load=False, progress=progress) - reference = ReferenceSet(args, load=False, progress=progress) - alignments = Aligner(args, load=False) - alignments.remove_species_from_alignments() - ogset.remove_species_from_ogs() - for mapping in progress._get_finished_mapping_folders(args.output_path): - species_name = mapping.split("04_mapping_")[-1] - logger.info('--- Addition of {} to all ogs ' - '---'.format(species_name)) - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, - species_name=species_name, load=False, progress=progress) - ogset.add_mapped_seq(mapper, species_name=species_name) - alignments.add_mapped_seq(ogset.mapped_ogs, species_name=species_name) - ogset.write_added_ogs_aa(folder_name="05_merge_OGs_aa") - ogset.write_added_ogs_dna(folder_name="05_merge_OGs_dna") - alignments.write_added_align_aa() - alignments.write_added_align_dna() - concat_alignment = alignments.concat_alignment() - if args.tree: - tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (not args.merge_all_mappings and progress.ref_ogs_01 and progress.ref_dna_02 and - progress.ref_align_03 and progress.mapping_04 and - (not progress.append_ogs_05 or not progress.align_06)): - ogset = OGSet(args, load=False, progress=progress) - reference = ReferenceSet(args, load=False, progress=progress) - alignments = Aligner(args, load=False) - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, load=False, progress=progress) + ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences + # ogset.ogs['OG1188079'].aa[0] = SeqRecord(seq=Seq('FLGMT ... + # ogset.ogs['OG1188079'].dna[0] = SeqRecord(seq=Seq('TTT + + reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) + # aa and dna of each species + # reference.ref['MNELE'].aa[0] = SeqRecord(seq=Seq('FLGM + + alignments = Aligner(args, ogset.ogs, step=args.step) # multiple sequence alignment of OGs + # alignments.alignments['OG1008242'].aa < instance (5 records of length 573) at 7f9189695c00> + # alignments.alignments['OG1008242'].aa[0] SeqRecord(seq=Seq('---------------MTDFDKL + + mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, step=args.step) # map reads onto OG alignments.remove_species_from_alignments() ogset.remove_species_from_ogs() ogset.add_mapped_seq(mapper) ogset.write_added_ogs_aa() ogset.write_added_ogs_dna() + # alignments = Aligner(args, ogset.mapped_ogs, load=True) alignments.add_mapped_seq(ogset.mapped_ogs) alignments.write_added_align_aa() alignments.write_added_align_dna() concat_alignment = alignments.concat_alignment() if args.tree: tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) - elif (not args.merge_all_mappings and progress.ref_ogs_01 and progress.ref_dna_02 and - progress.ref_align_03 and progress.mapping_04 and progress.append_ogs_05 and progress.align_06): - ogset = OGSet(args, load=False, progress=progress) - alignments = Aligner(args, load=False) - mapper = Mapper(args, og_set=ogset.ogs, load=False, progress=progress) + logger.info(str(tree.tree)) + + logger.info(' ------- Read2Tree finished -*- -------') + # print("done - all") + + if args.step == "marker": + oma_output = OMAOutputParser(args) + args.oma_output_path = oma_output.oma_output_path + ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences + reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) + alignments = Aligner(args, ogset.ogs) + print("done- marker") + logger.info(' ------- Read2Tree step marker finished -*- -------') + + if args.step == "map": + mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, step=args.step) alignments.remove_species_from_alignments() ogset.remove_species_from_ogs() ogset.add_mapped_seq(mapper) ogset.write_added_ogs_aa() ogset.write_added_ogs_dna() + # alignments = Aligner(args, ogset.mapped_ogs, load=True) alignments.add_mapped_seq(ogset.mapped_ogs) alignments.write_added_align_aa() alignments.write_added_align_dna() + + print("done- map") + logger.info(' ------- Read2Tree step map finished -*- -------') + + if args.step == "combine": + ogset = OGSet(args, step=args.step) + reference = ReferenceSet(args, step=args.step) + alignments = Aligner(args, step=args.step) # og set is not input of this, becuase we only read the aligned seq + alignments.remove_species_from_alignments() + ogset.remove_species_from_ogs() + mapping_folders_finished = get_finished_mapping_folders2(args.output_path) + for mapping in mapping_folders_finished: + species_name = mapping.split("04_mapping_")[-1] + logger.info('--- Addition of {} to all ogs ' + '---'.format(species_name)) + mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, + species_name=species_name, step=args.step) + ogset.add_mapped_seq(mapper, species_name=species_name) + alignments.add_mapped_seq(ogset.mapped_ogs, species_name=species_name) + ogset.write_added_ogs_aa(folder_name="05_merge_OGs_aa") + ogset.write_added_ogs_dna(folder_name="05_merge_OGs_dna") + alignments.write_added_align_aa() + alignments.write_added_align_dna() concat_alignment = alignments.concat_alignment() if args.tree: tree = TreeInference(args, concat_alignment=concat_alignment[0]) - print(tree.tree) + logger.info(str(tree.tree)) + + print("done- combine") + logger.info(' ------- Read2Tree step combined finished -*- -------') + + logger.info(' ------- Read2Tree finished -*- -------') + + print("done- main") + # logger.info('{}: Progress: ogs_dna {} | ref {} | ref_align {} | mapping {} | append_ogs {} | align {} ' + # .format(args.species_name, progress.ref_ogs_01, progress.ref_dna_02, + # progress.ref_align_03, progress.mapping_04, progress.append_ogs_05, progress.align_06)) + + # TODO: Check whether all the necessary binaries are available + # TODO: Check all given files and throw error if faulty logger.info('{}: Time taken {}'.format(args.species_name, timer() - t1)) From edd50a2a4fec5bcfab47b6dacc2a5742188e797f Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 11:53:34 -0400 Subject: [PATCH 25/50] [WIP] remove progress, define steps --- {read2tree => archive}/Progress.py | 0 archive/run_r2t.py | 6 +-- read2tree/Aligner.py | 6 +-- read2tree/Mapper.py | 45 +++++++++---------- read2tree/OGSet.py | 18 ++++---- read2tree/Reads.py | 70 +++++++++++++++--------------- read2tree/ReferenceSet.py | 14 +++--- read2tree/__init__.py | 2 +- read2tree/main.py | 28 ++++++------ 9 files changed, 95 insertions(+), 94 deletions(-) rename {read2tree => archive}/Progress.py (100%) diff --git a/read2tree/Progress.py b/archive/Progress.py similarity index 100% rename from read2tree/Progress.py rename to archive/Progress.py diff --git a/archive/run_r2t.py b/archive/run_r2t.py index 163378d..3318bfd 100644 --- a/archive/run_r2t.py +++ b/archive/run_r2t.py @@ -6,11 +6,11 @@ import sys -print("start run_r2t 2 ") - +print("start run_r2t 2 223 3 ") main(sys.argv[1:], exe_name=exe_name(), desc="descr") -print("finish run_r2t ") + +print("finish run_r2t ") a=1 \ No newline at end of file diff --git a/read2tree/Aligner.py b/read2tree/Aligner.py index 3efd5f3..da34cfb 100644 --- a/read2tree/Aligner.py +++ b/read2tree/Aligner.py @@ -27,7 +27,7 @@ class Aligner(object): - def __init__(self, args, og_set=None, load=True): + def __init__(self, args, og_set=None, step=None): self.args = args self.mapped_aligns = {} @@ -45,11 +45,11 @@ def __init__(self, args, og_set=None, load=True): self.alignments = Alignment() #self.placement_dic = {} - if load and og_set is not None: + if step== "all" or step== "1marker":#and og_set is not None: print('--- Alignment of {} OGs ---'.format(len(list(og_set.keys())))) self._og_set = og_set self.alignments = self._align(og_set) - else: + elif step =="3combine" or step =="2map": self.alignments = self._reload_alignments_from_folder() # print(self._get_codon_dict_og(og_set)) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index f6d9616..e5464eb 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -40,8 +40,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "minimap2" -samtools = "samtools" +minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2" +samtools = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin//samtools" class Mapper(object): @@ -49,7 +49,7 @@ class Mapper(object): Structure for reference """ - def __init__(self, args, ref_set=None, og_set=None, species_name=None, progress=None, + def __init__(self, args, ref_set=None, og_set=None, species_name=None, step=None, load=True): self.args = args self.elapsed_time = 0 @@ -72,25 +72,24 @@ def __init__(self, args, ref_set=None, og_set=None, species_name=None, progress= else: self._mapping_name = self._species_name - self.progress = progress + self.step = step self.all_cov = {} self.all_sc = {} self.read_og_set = {} - if load: # compute mapping + if step=="2map" or step=="all": # compute mapping if ref_set != None: - self.mapped_records = \ - self._map_reads_to_references(ref_set) + self.mapped_records = self._map_reads_to_references(ref_set) #if self.progress.get_mapping_status(): # self.progress.set_status('map') if self.mapped_records and og_set != None: self.og_records = self._sort_by_og() - else: # re-load already computed mapping + elif step =="3combine": # re-load already computed mapping if og_set != None and not self.args.merge_all_mappings: self.mapped_records = self._read_mapping_from_folder(ref_records=ref_set) self.og_records = self._sort_by_og() - elif (og_set != None and + elif (og_set != None and # todo to check self.args.merge_all_mappings and species_name != None): self.mapped_records = \ self._read_mapping_from_folder(mapping_name=self._mapping_name, ref_records=ref_set) @@ -298,11 +297,13 @@ def _map_reads_to_references(self, ref): reads = read_container.reads # print(os.path.getsize(reads[0])) - if self.args.single_mapping: - references = [self.args - .single_mapping.split("/")[-1].split("_")[0]] - else: - references = list(ref.keys()) + # if self.args.single_mapping: + # references = [self.args + # .single_mapping.split("/")[-1].split("_")[0]] + # else: + # references = list(ref.keys()) + # + references = list(ref.keys()) # Going through provided references and starting mapping for species in tqdm(references, @@ -570,8 +571,8 @@ def _post_process_read_mapping(self, ref_file, sam_file): #outfile_name = os.path.join(tmp_folder, # ref_file.split('/')[-1].split('.')[0] + # "_post") - if self.args.single_mapping: - self.logger.debug("{}: --- POSTPROCESSING MAPPING ---".format(self._species_name)) + # if self.args.single_mapping: + # self.logger.debug("{}: --- POSTPROCESSING MAPPING ---".format(self._species_name)) # ngmlr doesn't have the option to write in bam file directly #if 'sam' in bam_file.split(".")[-1]: @@ -584,8 +585,8 @@ def _post_process_read_mapping(self, ref_file, sam_file): else: self.logger.debug("Sam file is not generated", sam_file) - if self.args.single_mapping: - self.logger.error("single mapping is not tested with this version ") + # if self.args.single_mapping: + # self.logger.error("single mapping is not tested with this version ") #self.logger.debug("{}: ---- Samtools view completed".format(self._species_name)) if os.path.exists(sam_file_base+".bam"): @@ -596,16 +597,16 @@ def _post_process_read_mapping(self, ref_file, sam_file): else: self.logger.debug("bam file is not generated", sam_file_base+".bam") - if self.args.single_mapping: - self.logger.debug("{}: ---- Samtools sort completed".format(self._species_name)) + # if self.args.single_mapping: + # self.logger.debug("{}: ---- Samtools sort completed".format(self._species_name)) if os.path.exists(sam_file_base + "_sorted.bam"): line_samtools_index= samtools+' index -@ ' + str(self.args.threads) + ' ' + sam_file_base + "_sorted.bam" self._output_shell(line_samtools_index) self.logger.debug(line_samtools_index) - if self.args.single_mapping: - self.logger.warning("single_mapping is not tested in this version.") #debug("{}: ---- Samtools index completed".format(self._species_name)) + # if self.args.single_mapping: + # self.logger.warning("single_mapping is not tested in this version.") #debug("{}: ---- Samtools index completed".format(self._species_name)) # self._rm_file(bam_file, ignore_error=True) #f self.args.debug: diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 370558b..08ad9bc 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -33,7 +33,7 @@ class OGSet(object): - def __init__(self, args, oma_output=None, load=True, progress=None): + def __init__(self, args, oma_output=None, step=None): self.args = args self.logger = logging.getLogger(__name__) @@ -49,7 +49,7 @@ def __init__(self, args, oma_output=None, load=True, progress=None): self._ham_analysis = None - self.progress = progress + self.step = step # self.progress.get_status(species_name=self._species_name) if self.args.remove_species_mapping: @@ -69,9 +69,11 @@ def __init__(self, args, oma_output=None, load=True, progress=None): # self.mapped_ogs = self._reload_ogs_from_folder( # folder_suffix="05_ogs_map_" + self._species_name) # self.ogs = self.mapped_ogs - if not load and self.progress.ref_ogs_01: + #if not load and self.progress.ref_ogs_01: + if step == "2map" or step=="3combine": self.ogs = self._reload_ogs_from_folder() - elif load and oma_output is not None: + + elif step == "all" or step == "1marker" : #load and oma_output is not None: self.min_species = oma_output.min_species self.oma = oma_output self.ogs = oma_output.ogs @@ -176,10 +178,8 @@ def _load_ogs(self): # name = file.split("/")[-1].split(".")[0] ogs[name] = OG() ogs[name].aa = self._get_aa_records(name, records) - output_file_aa = os.path.join(orthologous_groups_aa, - name + ".fa") - output_file_dna = os.path.join(orthologous_groups_dna, - name + ".fa") + output_file_aa = os.path.join(orthologous_groups_aa, name + ".fa") + output_file_dna = os.path.join(orthologous_groups_dna, name + ".fa") if source: try: @@ -451,7 +451,7 @@ def _get_best_record(self, cons_og, og_sc, species_name): best_record_dna.id = species_name return (best_record_aa, best_record_dna) else: - return none + return None def _generate_seq_completeness(self, seqC, mapper, og, best_record_dna): if self.args.remove_species_ogs: diff --git a/read2tree/Reads.py b/read2tree/Reads.py index 3848a0b..22bc664 100644 --- a/read2tree/Reads.py +++ b/read2tree/Reads.py @@ -27,9 +27,9 @@ def __init__(self, args, load=True): self.elapsed_time = 0 self.total_reads = 0 - self.split_len = args.split_len - self.split_overlap = args.split_overlap - self.split_min_read_len = args.split_min_read_len + #self.split_len = args.split_len + #self.split_overlap = args.split_overlap + #self.split_min_read_len = args.split_min_read_len self.logger = logging.getLogger(__name__) @@ -44,34 +44,34 @@ def __init__(self, args, load=True): self._reads = self.args.reads self._species_name = self.args.species_name - if load: - if len(self.args.reads) == 2 and self.args.check_mate_pairing: - mate_pairs = self.check_read_consistency(self._reads) - if mate_pairs: - tmp = self.select_mates_from_reads(self._reads, - mate_pairs) - self._reads = tmp - - if self.args.split_reads: - print('--- Splitting reads from {} ---'.format(self._reads)) - self.logger.info('{}: --- Splitting reads from {} ---' - .format(self._species_name, self._reads)) - # print(memory_usage(self.process_reads)) - # self.split_reads = self._write_to_tmp_file(self\ - # .process_reads()) - self.reads = self.process_reads() - else: - self.reads = self._reads - - if self.args.sample_reads: - print('--- Sampling reads from {} ---'.format(self.reads)) - self.logger.info('{}: --- Sampling reads from {} ---' - .format(self._species_name, self.reads)) - self.reads = self.sample_from_reads(self.reads) - else: - self.reads = self.reads + #if load: + # if len(self.args.reads) == 2 and self.args.check_mate_pairing: + # mate_pairs = self.check_read_consistency(self._reads) + # if mate_pairs: + # tmp = self.select_mates_from_reads(self._reads, + # mate_pairs) + # self._reads = tmp + + # if self.args.split_reads: + # print('--- Splitting reads from {} ---'.format(self._reads)) + # self.logger.info('{}: --- Splitting reads from {} ---' + # .format(self._species_name, self._reads)) + # # print(memory_usage(self.process_reads)) + # # self.split_reads = self._write_to_tmp_file(self\ + # # .process_reads()) + # self.reads = self.process_reads() + # else: + # self.reads = self._reads + self.reads = self._reads + if self.args.sample_reads: + print('--- Sampling reads from {} ---'.format(self.reads)) + self.logger.info('{}: --- Sampling reads from {} ---' + .format(self._species_name, self.reads)) + self.reads = self.sample_from_reads(self.reads) else: - self.reads = self._reads + self.reads = self.reads + #else: + # self.reads = self._reads def process_reads(self): ''' @@ -126,11 +126,11 @@ def process_reads(self): end = time.time() self.elapsed_time = end - start - self.logger.info('{}: Reads larger than {} were split into {} bp long ' - 'fragments with an overlap of {} bp.'.format( - self._species_name, self.split_min_read_len, - self.split_len, - self.split_overlap)) + # self.logger.info('{}: Reads larger than {} were split into {} bp long ' + # 'fragments with an overlap of {} bp.'.format( + # self._species_name, self.split_min_read_len, + # self.split_len, + # self.split_overlap)) self.logger.info('{}: {} reads were split into {} reads.' .format(self._species_name, total_reads, total_new_reads)) diff --git a/read2tree/ReferenceSet.py b/read2tree/ReferenceSet.py index 30882a4..2817421 100644 --- a/read2tree/ReferenceSet.py +++ b/read2tree/ReferenceSet.py @@ -14,15 +14,13 @@ from Bio import SeqIO from Bio.SeqIO.FastaIO import FastaWriter -from read2tree.Progress import Progress - class ReferenceSet(object): ''' Structure for reference ''' - def __init__(self, args, og_set=None, load=True, progress=None): + def __init__(self, args, og_set=None, step=None): """ :param args: list of arguments from command line @@ -30,16 +28,18 @@ def __init__(self, args, og_set=None, load=True, progress=None): :param load: set to True when reference loaded from folder/file of list of arguments """ self.ref = {} - self.load = load + #self.load = load self.args = args - self.progress = progress + self.step = step self.logger = logging.getLogger(__name__) self._species_name = self.args.species_name - if load is False: + #if load is False: + if step == "2map": self.ref = self._load_records_folder() - elif og_set is not None and load is True: + #elif og_set is not None and load is True: + elif step == "all" or step == "1marker": # self.ref = self._generate_reference(og_set) self.write() # self.progress.set_status('ref') diff --git a/read2tree/__init__.py b/read2tree/__init__.py index 1d8fe87..c22e791 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -7,7 +7,7 @@ logging.getLogger(__name__).addHandler(logging.NullHandler()) __version__ = '1.5.0' -__copyright__ = 'read2tree (C) 2017-{:d} David Dylus' \ +__copyright__ = 'read2tree (C) 2017-{:d} David Dylus ' \ .format(date.today().year) # path = './log.yaml' diff --git a/read2tree/main.py b/read2tree/main.py index 551af7a..61b7fb4 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -158,7 +158,7 @@ def parse_args(argv, exe_name, desc): 'output concatenated alignment!') arg_parser.add_argument('--step', default="all", - help='[Default is all] marker map combine ') + help='[Default is all 1marker 2map 3combine ') # arg_parser.add_argument('--merge_all_mappings', action='store_true', # help='[Default is off] In case multiple species were mapped to ' @@ -228,7 +228,7 @@ def parse_args(argv, exe_name, desc): if args.species_name: _species_name = args.species_name - if args.step == "combine": # todo why is needed? + if args.step == "3combine": # todo why is needed? _species_name = 'merge' args.reads = _reads @@ -308,14 +308,14 @@ def main(argv, exe_name, desc=''): logger.info("Running r2t in mode " + args.step) if os.path.exists(args.output_path): - if args.step == "all" or args.step == "marker": + if args.step == "all" or args.step == "1marker": logger.error( - "the output folder exist" + args.output_path + ". Since you are running r2t in all mode, you need to specify output folder which will created by r2t.") + "the output folder exist " + args.output_path + ". Since you are running r2t in all mode, you need to specify output folder which will created by r2t.") sys.exit() else: os.makedirs(args.output_path) - if args.step == "all" or args.step == "map": + if args.step == "all" or args.step == "2map": if not args.reads: logger.error("reads are not provided in mode " + str(args.step)) sys.exit() @@ -363,16 +363,16 @@ def main(argv, exe_name, desc=''): logger.info(' ------- Read2Tree finished -*- -------') # print("done - all") - if args.step == "marker": + if args.step == "1marker": oma_output = OMAOutputParser(args) args.oma_output_path = oma_output.oma_output_path ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) alignments = Aligner(args, ogset.ogs) - print("done- marker") - logger.info(' ------- Read2Tree step marker finished -*- -------') + print("done- 1marker") + logger.info(' ------- Read2Tree step 1marker finished -*- -------') - if args.step == "map": + if args.step == "2map": mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, step=args.step) alignments.remove_species_from_alignments() ogset.remove_species_from_ogs() @@ -384,10 +384,10 @@ def main(argv, exe_name, desc=''): alignments.write_added_align_aa() alignments.write_added_align_dna() - print("done- map") + print("done- 2map") logger.info(' ------- Read2Tree step map finished -*- -------') - if args.step == "combine": + if args.step == "3combine": ogset = OGSet(args, step=args.step) reference = ReferenceSet(args, step=args.step) alignments = Aligner(args, step=args.step) # og set is not input of this, becuase we only read the aligned seq @@ -411,12 +411,12 @@ def main(argv, exe_name, desc=''): tree = TreeInference(args, concat_alignment=concat_alignment[0]) logger.info(str(tree.tree)) - print("done- combine") - logger.info(' ------- Read2Tree step combined finished -*- -------') + print("done- 3combine") + logger.info(' ------- Read2Tree step 3combined finished -*- -------') logger.info(' ------- Read2Tree finished -*- -------') - print("done- main") + print("done- main ") # logger.info('{}: Progress: ogs_dna {} | ref {} | ref_align {} | mapping {} | append_ogs {} | align {} ' # .format(args.species_name, progress.ref_ogs_01, progress.ref_dna_02, # progress.ref_align_03, progress.mapping_04, progress.append_ogs_05, progress.align_06)) From 1e21c40762ddada63c0047e78cdc1a3df939f412 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 12:49:18 -0400 Subject: [PATCH 26/50] [WIP] debug --- read2tree/Mapper.py | 4 ++-- read2tree/main.py | 39 +++++++++++++++++++++++---------------- 2 files changed, 25 insertions(+), 18 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index e5464eb..6eeb563 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -79,11 +79,11 @@ def __init__(self, args, ref_set=None, og_set=None, species_name=None, step=None self.read_og_set = {} if step=="2map" or step=="all": # compute mapping - if ref_set != None: + if ref_set != None: #todo else self.mapped_records = self._map_reads_to_references(ref_set) #if self.progress.get_mapping_status(): # self.progress.set_status('map') - if self.mapped_records and og_set != None: + if self.mapped_records and og_set != None: #todo else self.og_records = self._sort_by_og() elif step =="3combine": # re-load already computed mapping if og_set != None and not self.args.merge_all_mappings: diff --git a/read2tree/main.py b/read2tree/main.py index 61b7fb4..884e5c7 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -300,17 +300,17 @@ def main(argv, exe_name, desc=''): t1 = timer() # Parse args = parse_args(argv, exe_name, desc) - logger.info('{}: ------- NEW RUN -------'.format(args.species_name)) + x = ', '.join("{!s}={!r}".format(key, val) for (key, val) in vars(args).items()) # todo why calling Progress twice? logger.info('{}: read2tree was run with: {}'.format(args.species_name, x)) - logger.info("Running r2t in mode " + args.step) + logger.info("Running read2tree in mode " + args.step) if os.path.exists(args.output_path): if args.step == "all" or args.step == "1marker": logger.error( - "the output folder exist " + args.output_path + ". Since you are running r2t in all mode, you need to specify output folder which will created by r2t.") + "the output folder exist " + args.output_path + ". Since you are running read2tree in all mode, you need to specify output folder which will created by read2tree.") sys.exit() else: os.makedirs(args.output_path) @@ -323,7 +323,7 @@ def main(argv, exe_name, desc=''): list): # in parse_args it is input read is converted to a list, if there are 2 - paired end for read_file in args.reads: if not os.path.isfile(read_file): - logger.error("read file doesn't exist" + read_file) + logger.error("read file doesn't exist " + read_file) sys.exit() elif isinstance(args.reads, str): if not os.path.isfile(args.reads): @@ -331,17 +331,18 @@ def main(argv, exe_name, desc=''): sys.exit() if args.step == "all": + logger.info('{}: ------- NEW RUN -------'.format(args.species_name)) oma_output = OMAOutputParser(args) args.oma_output_path = oma_output.oma_output_path - ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences + ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences. # Write 01_ # ogset.ogs['OG1188079'].aa[0] = SeqRecord(seq=Seq('FLGMT ... # ogset.ogs['OG1188079'].dna[0] = SeqRecord(seq=Seq('TTT - reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) + reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) # write 02_ # aa and dna of each species # reference.ref['MNELE'].aa[0] = SeqRecord(seq=Seq('FLGM - alignments = Aligner(args, ogset.ogs, step=args.step) # multiple sequence alignment of OGs + alignments = Aligner(args, ogset.ogs, step=args.step) # multiple sequence alignment of OGs # write 03_ # alignments.alignments['OG1008242'].aa < instance (5 records of length 573) at 7f9189695c00> # alignments.alignments['OG1008242'].aa[0] SeqRecord(seq=Seq('---------------MTDFDKL @@ -363,16 +364,24 @@ def main(argv, exe_name, desc=''): logger.info(' ------- Read2Tree finished -*- -------') # print("done - all") + if args.step == "1marker": + logger.info('{}: ------- NEW RUN -------'.format(args.species_name)) oma_output = OMAOutputParser(args) args.oma_output_path = oma_output.oma_output_path ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) - alignments = Aligner(args, ogset.ogs) + alignments = Aligner(args, ogset.ogs, step=args.step) print("done- 1marker") logger.info(' ------- Read2Tree step 1marker finished -*- -------') if args.step == "2map": + logger.info('{}: ------- Read2tree RUN step 2map (after running 1marker) -------'.format(args.species_name)) + oma_output = OMAOutputParser(args) + args.oma_output_path = oma_output.oma_output_path + ogset = OGSet(args, oma_output=oma_output, step=args.step) # Generate the OGs with their DNA sequences + reference = ReferenceSet(args, og_set=ogset.ogs, step=args.step) + alignments = Aligner(args, ogset.ogs, step=args.step) # multiple sequence alignment of OGs mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, step=args.step) alignments.remove_species_from_alignments() ogset.remove_species_from_ogs() @@ -380,12 +389,12 @@ def main(argv, exe_name, desc=''): ogset.write_added_ogs_aa() ogset.write_added_ogs_dna() # alignments = Aligner(args, ogset.mapped_ogs, load=True) - alignments.add_mapped_seq(ogset.mapped_ogs) - alignments.write_added_align_aa() - alignments.write_added_align_dna() + # alignments.add_mapped_seq(ogset.mapped_ogs) + # alignments.write_added_align_aa() + # alignments.write_added_align_dna() print("done- 2map") - logger.info(' ------- Read2Tree step map finished -*- -------') + logger.info(' ------- Read2Tree step 2map finished -*- -------') if args.step == "3combine": ogset = OGSet(args, step=args.step) @@ -412,16 +421,14 @@ def main(argv, exe_name, desc=''): logger.info(str(tree.tree)) print("done- 3combine") - logger.info(' ------- Read2Tree step 3combined finished -*- -------') + logger.info(' ------- Read2Tree step 3combine finished -*- -------') logger.info(' ------- Read2Tree finished -*- -------') print("done- main ") - # logger.info('{}: Progress: ogs_dna {} | ref {} | ref_align {} | mapping {} | append_ogs {} | align {} ' - # .format(args.species_name, progress.ref_ogs_01, progress.ref_dna_02, - # progress.ref_align_03, progress.mapping_04, progress.append_ogs_05, progress.align_06)) # TODO: Check whether all the necessary binaries are available # TODO: Check all given files and throw error if faulty logger.info('{}: Time taken {}'.format(args.species_name, timer() - t1)) + # todo remove sam file \ No newline at end of file From c4823c82fd3982dffed0366344ece5e8375fefc4 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 12:57:21 -0400 Subject: [PATCH 27/50] debug --- read2tree/Mapper.py | 12 ++++++------ read2tree/main.py | 1 - 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 6eeb563..90b2af2 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -86,14 +86,14 @@ def __init__(self, args, ref_set=None, og_set=None, species_name=None, step=None if self.mapped_records and og_set != None: #todo else self.og_records = self._sort_by_og() elif step =="3combine": # re-load already computed mapping - if og_set != None and not self.args.merge_all_mappings: + if og_set != None: #and not self.args.merge_all_mappings: self.mapped_records = self._read_mapping_from_folder(ref_records=ref_set) self.og_records = self._sort_by_og() - elif (og_set != None and # todo to check - self.args.merge_all_mappings and species_name != None): - self.mapped_records = \ - self._read_mapping_from_folder(mapping_name=self._mapping_name, ref_records=ref_set) - self.og_records = self._sort_by_og() + # elif (og_set != None and # todo to check + # self.args.merge_all_mappings and species_name != None): + # self.mapped_records = \ + # self._read_mapping_from_folder(mapping_name=self._mapping_name, ref_records=ref_set) + # self.og_records = self._sort_by_og() # if not (self.mapped_records or hasattr(self,"og_records")): # log_line='mapped_records or og_records is empty! are the aligner and samtools installed correctly? are there enough marker genes close to the species of interest? If it is a new run, make sure there is no output folder or a file mplog.log' diff --git a/read2tree/main.py b/read2tree/main.py index 884e5c7..9ecf730 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -388,7 +388,6 @@ def main(argv, exe_name, desc=''): ogset.add_mapped_seq(mapper) ogset.write_added_ogs_aa() ogset.write_added_ogs_dna() - # alignments = Aligner(args, ogset.mapped_ogs, load=True) # alignments.add_mapped_seq(ogset.mapped_ogs) # alignments.write_added_align_aa() # alignments.write_added_align_dna() From ad1ef0e5eed9c925aac13f2e4b4ebac94241f9d2 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 19:05:30 -0400 Subject: [PATCH 28/50] debug --- read2tree/Aligner.py | 22 ++++++++++++---------- read2tree/Mapper.py | 24 +++++++++++------------- read2tree/main.py | 5 +++-- 3 files changed, 26 insertions(+), 25 deletions(-) diff --git a/read2tree/Aligner.py b/read2tree/Aligner.py index da34cfb..dbee3fb 100644 --- a/read2tree/Aligner.py +++ b/read2tree/Aligner.py @@ -46,7 +46,7 @@ def __init__(self, args, og_set=None, step=None): #self.placement_dic = {} if step== "all" or step== "1marker":#and og_set is not None: - print('--- Alignment of {} OGs ---'.format(len(list(og_set.keys())))) + logger.info('--- Alignment of {} OGs ---'.format(len(list(og_set.keys())))) self._og_set = og_set self.alignments = self._align(og_set) elif step =="3combine" or step =="2map": @@ -149,29 +149,31 @@ def add_mapped_seq(self, ogset_add, species_name=None): num_append_seq = 0 if not species_name: species_name = self._species_name - print('--- Add inferred mapped sequence back to alignment ---') + logger.info('--- Add inferred mapped sequence back to alignment ---') # iterate through all existing ogs - for name_og, align in tqdm(self.alignments.items(), - desc='Adding mapped seq to alignments', unit=' alignments'): + for name_og, align in tqdm(self.alignments.items(), desc='Adding mapped seq to alignments', unit=' alignments'): align_filt = align if len(align_filt.aa) >= 2: # get all species that are not mapped from original alignment if name_og in ogset_add.keys(): # find mapped records from appended records in OGSet - map_record_aa = [r for r in ogset_add[name_og].aa if species_name in r.id] - # print(map_records_aa) - map_record_dna = [r for r in ogset_add[name_og].dna if species_name in r.id] + map_record_aa = [r for r in ogset_add[name_og].aa if species_name == r.id] + # print(map_records_aa) # todo this is important which is solved it's a bug in when the sample name is small , probably should be species_name == r.id species_name='3' id='ASTMX01499_OG1003046' + + map_record_dna = [r for r in ogset_add[name_og].dna if species_name == r.id] if map_record_aa and map_record_dna: - ref_species = self._get_species_id(map_record_aa[0]) + ref_species = self._get_species_id(map_record_aa[0]) # why the zeroth? how about the rest? self.mapped_aligns[name_og] = Alignment() - self.mapped_aligns[name_og].aa = self._add_mapseq_align(align_filt.aa, map_record_aa[0], ref_species, species_name) + self.mapped_aligns[name_og].aa = self._add_mapseq_align(align_filt.aa, map_record_aa[0], ref_species, species_name) # why the zeroth? how about the rest? self.mapped_aligns[name_og].dna = self._add_mapseq_align(align_filt.dna, map_record_dna[0], ref_species, species_name) num_append_seq = num_append_seq + 1 - elif self.args.keep_all_ogs: + + elif self.args.keep_all_ogs: # todo both self.args.keep_all_ogs and "map_record_aa and map_record_dna" are true althout it is eilf self.mapped_aligns[name_og] = Alignment() self.mapped_aligns[name_og].aa = align_filt.aa self.mapped_aligns[name_og].dna = align_filt.dna + elif self.args.keep_all_ogs: self.mapped_aligns[name_og] = Alignment() self.mapped_aligns[name_og].aa = align_filt.aa diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 90b2af2..bc916a4 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -86,7 +86,7 @@ def __init__(self, args, ref_set=None, og_set=None, species_name=None, step=None if self.mapped_records and og_set != None: #todo else self.og_records = self._sort_by_og() elif step =="3combine": # re-load already computed mapping - if og_set != None: #and not self.args.merge_all_mappings: + if og_set != None: #and not self.args.merge_all_mappings: # todo else self.mapped_records = self._read_mapping_from_folder(ref_records=ref_set) self.og_records = self._sort_by_og() # elif (og_set != None and # todo to check @@ -177,14 +177,16 @@ def _read_mapping_from_folder(self, mapping_name=None, ref_records=None): :return: dictionary with key og_name and value sequences mapped to each species """ - print('--- Retrieve mapped consensus sequences ---') + map_reads_species = {} if not mapping_name: mapping_name = self._mapping_name - in_folder = os.path.join(self.args.output_path, - "04_mapping_"+mapping_name) + self.logger.debug('--- Retrieve mapped consensus sequences ---'+str(mapping_name)) + in_folder = os.path.join(self.args.output_path, "04_mapping_"+mapping_name) + bam_files = glob.glob(os.path.join(in_folder, "*.bam")) - if self.args.min_cons_coverage >= 2 and bam_files: + if self.args.min_cons_coverage >= 2 and bam_files: # default self.args.min_cons_coverage is 1 + self.logger.debug('Generating consensus from bam files' ) for file in tqdm(bam_files, desc='Generating consensus from bam files ', unit=' species'): species = file.split("/")[-1].split("_")[0] ref_file = os.path.join(self.args.output_path, '02_ref_dna',species+'_OGs.fa') @@ -204,20 +206,16 @@ def _read_mapping_from_folder(self, mapping_name=None, ref_records=None): cov = Coverage(self.args) cov.get_coverage_bam(file) - cov.write_coverage_bam(os.path.join( - in_folder, ref_file.split('/')[-1].split('.')[0] + "_cov.txt")) + cov.write_coverage_bam(os.path.join(in_folder, ref_file.split('/')[-1].split('.')[0] + "_cov.txt")) self.all_cov.update(cov.coverage) seqC = SeqCompleteness(mapped_ref=ref_records[species].dna) seqC.get_seq_completeness(map_reads_species[species].dna) - seqC.write_seq_completeness(os - .path.join(in_folder, - species + "_OGs_sc.txt")) + seqC.write_seq_completeness(os.path.join(in_folder,species + "_OGs_sc.txt")) self.all_sc.update(seqC.seq_completeness) else: - for file in tqdm(glob.glob(os.path.join(in_folder, "*_consensus.fa")), - desc='Loading consensus read mappings ', - unit=' species'): + self.logger.debug('Loading consensus read mappings '+in_folder+ "*_consensus.fa") + for file in tqdm(glob.glob(os.path.join(in_folder, "*_consensus.fa")), desc='Loading consensus read mappings ', unit=' species'): species = file.split("/")[-1].split("_")[0] map_reads_species[species] = Reference() fasta_reader = FastxReader(file) diff --git a/read2tree/main.py b/read2tree/main.py index 9ecf730..6751b31 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -401,13 +401,14 @@ def main(argv, exe_name, desc=''): alignments = Aligner(args, step=args.step) # og set is not input of this, becuase we only read the aligned seq alignments.remove_species_from_alignments() ogset.remove_species_from_ogs() - mapping_folders_finished = get_finished_mapping_folders2(args.output_path) + mapping_folders_finished = get_finished_mapping_folders2(args.output_path) # this is done based on "_all_cov.txt" files for mapping in mapping_folders_finished: species_name = mapping.split("04_mapping_")[-1] logger.info('--- Addition of {} to all ogs ' '---'.format(species_name)) - mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, + mapper = Mapper(args, og_set=ogset.ogs, ref_set=reference.ref, # for those OGs reported in {sample}_all_cov.txt file species_name=species_name, step=args.step) + ogset.add_mapped_seq(mapper, species_name=species_name) alignments.add_mapped_seq(ogset.mapped_ogs, species_name=species_name) ogset.write_added_ogs_aa(folder_name="05_merge_OGs_aa") From 7402d0316de7cb00137af0ea2708ea34463730d1 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 19:10:38 -0400 Subject: [PATCH 29/50] update readme --- README.md | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/README.md b/README.md index 7389035..83dc5ea 100644 --- a/README.md +++ b/README.md @@ -95,7 +95,7 @@ cat marker_genes/*.fna > dna_ref.fa ### output -The output of Read2Tree is the concatenated alignments as a fasta file where each record corresponds to one species. We also provide the option `--tree` for inferring the species tree using IQTREE as defualt. +The output of Read2Tree is the concatenated alignments as a fasta file where each record corresponds to one species. We also provide the option `--tree` for inferring the species tree using IQTREE as default. ### Single species mode @@ -104,12 +104,23 @@ read2tree --tree --standalone_path marker_genes/ --reads read_1.fastq read_2.fas ``` ### Multiple species mode + +#### step1 +``` +read2tree --step 1marker --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --debug # this creates just the reference folder 01 - 03 +``` + +#### step2 +The following could be run in parallel. ``` -read2tree --standalone_path marker_genes/ --output_path output --reference --dna_reference dna_ref.fa # this creates just the reference folder 01 - 03 -read2tree --standalone_path marker_genes/ --output_path output --reads species1_R1.fastq species2_R2.fastq -read2tree --standalone_path marker_genes/ --output_path output --reads species2_R1.fastq species2_R2.fastq -read2tree --standalone_path marker_genes/ --output_path output --reads species3_R1.fastq species3_R2.fastq -read2tree --standalone_path marker_genes/ --output_path output --merge_all_mappings --tree +read2tree --step 2map --standalone_path marker_genes --dna_reference dna_ref.fa --reads species1_R1.fastq species2_R2.fastq --output_path output --debug +read2tree --step 2map --standalone_path marker_genes --dna_reference dna_ref.fa --reads species2_R1.fastq species2_R2.fastq --output_path output --debug +read2tree --step 2map --standalone_path marker_genes --dna_reference dna_ref.fa --reads species3_R1.fastq species3_R2.fastq --output_path output --debug +``` + +#### step3 +``` +read2tree --step 3combine --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --tree --debug ``` ### bootstraping @@ -136,12 +147,12 @@ The goal of this test example is to infer species tree for Mus musculus using it ``` cd tests -read2tree --debug --tree --standalone_path marker_genes/ --reads sample_1.fastq sample_2.fastq --output_path output/ --dna_reference dna_ref.fa +read2tree --tree --standalone_path marker_genes/ --reads read_1.fastq read_2.fastq --output_path output --dna_reference dna_ref.fa ``` #### Run test example using docker - +(to be updated ) ``` docker run --rm -i -v $PWD/tests:/input -v $PWD/tests/:/reads -v $PWD/output:/out -v $PWD/run:/run dessimozlab/read2tree:latest --tree --standalone_path /input/marker_genes --dna_reference /input/cds-marker_genes.fasta.gz --reads /reads/sample_1.fastq --output_path /out ``` From 7e575826e194eef180a9d67bfc953c5b73cbb6bd Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 19:11:14 -0400 Subject: [PATCH 30/50] update readme --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 83dc5ea..a6b0e33 100644 --- a/README.md +++ b/README.md @@ -107,7 +107,7 @@ read2tree --tree --standalone_path marker_genes/ --reads read_1.fastq read_2.fas #### step1 ``` -read2tree --step 1marker --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --debug # this creates just the reference folder 01 - 03 +read2tree --step 1marker --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --debug ``` #### step2 From aa5bb441847b3088ad669f48347e465d37da4559 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 19:12:26 -0400 Subject: [PATCH 31/50] update --- read2tree/Mapper.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index bc916a4..b3fc275 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -40,8 +40,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2" -samtools = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin//samtools" +minimap2_ex= "minimap2" +samtools = "samtools" class Mapper(object): From 5765deb32cc3dcd9f6d3d38d0419ca24416dff9e Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 7 May 2024 20:57:11 -0400 Subject: [PATCH 32/50] update --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index a6b0e33..bca3185 100644 --- a/README.md +++ b/README.md @@ -147,7 +147,7 @@ The goal of this test example is to infer species tree for Mus musculus using it ``` cd tests -read2tree --tree --standalone_path marker_genes/ --reads read_1.fastq read_2.fastq --output_path output --dna_reference dna_ref.fa +read2tree --tree --standalone_path marker_genes/ --reads sample_1.fastq sample_2.fastq --output_path output --dna_reference dna_ref.fa ``` From 8cd78fff0da77bf7dded89bd67196f44956dae03 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Mon, 21 Oct 2024 22:08:23 -0400 Subject: [PATCH 33/50] deactivating dl gene markers via REST --- read2tree/OGSet.py | 102 ++++++++++++++++++++++----------------------- 1 file changed, 51 insertions(+), 51 deletions(-) diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 08ad9bc..13b64af 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -264,57 +264,57 @@ def _get_dna_from_h5(self, record): id=record.id, description="") - def _get_dna_from_REST(self, record): - """ - - :param record: - :param name: - :return: - """ - tmp_id = re.sub(r'\..*', '', record.id.split("_")[0]) - use_id = re.sub(r'\W+', '', tmp_id) - try: - oma_record = requests.get(API_URL + "/protein/" + use_id + "/") - except requests.exceptions.RequestException: - self.logger.debug('DNA not found for {}.'.format(use_id)) - pass - else: - seq = oma_record.json()['cdna'] - rec_id = oma_record.json()['omaid'] - - cleaned_seq = self._clean_DNA_seq(seq) - # else: - # cleaned_seq = dna_record.seq - return SeqRecord.SeqRecord(cleaned_seq, record.id, - description="", name="") - - def _get_dna_from_REST_bulk(self, records, og_name): - """ - - :param record: - :param name: - :return: - """ - record_ids = [r.id for r in records] - dna_records = [] - try: - reply = requests.post('https://omabrowser.org/api/protein/bulk_retrieve/', - json={"ids": record_ids}, - headers={'User-Agent': 'read2tree/'+read2tree_version}) - except requests.exceptions.RequestException as exception_type: - self.logger.warning('DNA not found probably for '+str(record_ids[0])+'. The reason is '+str(exception_type)) - pass - else: - group_members = reply.json() - for memb in group_members: - # print(">{}\n{}\n\n".format(memb['omaid'], memb['cdna'])) - seq = memb['target']['cdna'] - rec_id = memb['target']['omaid']+"_"+og_name - cleaned_seq = self._clean_DNA_seq(seq) - # print(cleaned_seq) - dna_records.append(SeqRecord.SeqRecord(cleaned_seq, id=rec_id, - description="", name="")) - return dna_records + # def _get_dna_from_REST(self, record): + # """ + # + # :param record: + # :param name: + # :return: + # """ + # tmp_id = re.sub(r'\..*', '', record.id.split("_")[0]) + # use_id = re.sub(r'\W+', '', tmp_id) + # try: + # oma_record = requests.get(API_URL + "/protein/" + use_id + "/") + # except requests.exceptions.RequestException: + # self.logger.debug('DNA not found for {}.'.format(use_id)) + # pass + # else: + # seq = oma_record.json()['cdna'] + # rec_id = oma_record.json()['omaid'] + # + # cleaned_seq = self._clean_DNA_seq(seq) + # # else: + # # cleaned_seq = dna_record.seq + # return SeqRecord.SeqRecord(cleaned_seq, record.id, + # description="", name="") + # + # def _get_dna_from_REST_bulk(self, records, og_name): + # """ + # + # :param record: + # :param name: + # :return: + # """ + # record_ids = [r.id for r in records] + # dna_records = [] + # try: + # reply = requests.post('https://omabrowser.org/api/protein/bulk_retrieve/', + # json={"ids": record_ids}, + # headers={'User-Agent': 'read2tree/'+read2tree_version}) + # except requests.exceptions.RequestException as exception_type: + # self.logger.warning('DNA not found probably for '+str(record_ids[0])+'. The reason is '+str(exception_type)) + # pass + # else: + # group_members = reply.json() + # for memb in group_members: + # # print(">{}\n{}\n\n".format(memb['omaid'], memb['cdna'])) + # seq = memb['target']['cdna'] + # rec_id = memb['target']['omaid']+"_"+og_name + # cleaned_seq = self._clean_DNA_seq(seq) + # # print(cleaned_seq) + # dna_records.append(SeqRecord.SeqRecord(cleaned_seq, id=rec_id, + # description="", name="")) + # return dna_records def _get_dna_from_fasta(self, record, db): try: From 912f7481bb50e2dc1966aca7735e48809c48de87 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 09:12:09 -0400 Subject: [PATCH 34/50] deactivating REST cleanup --- archive/run_r2t.py | 2 +- read2tree/OGSet.py | 194 ++++++++++++++++++++++-------------------- read2tree/__init__.py | 3 +- 3 files changed, 107 insertions(+), 92 deletions(-) diff --git a/archive/run_r2t.py b/archive/run_r2t.py index 3318bfd..8ecfc7c 100644 --- a/archive/run_r2t.py +++ b/archive/run_r2t.py @@ -9,8 +9,8 @@ print("start run_r2t 2 223 3 ") main(sys.argv[1:], exe_name=exe_name(), desc="descr") +print("finish run_r2t ") -print("finish run_r2t ") a=1 \ No newline at end of file diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 13b64af..0983bfe 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -8,12 +8,13 @@ import os import re import pyham -import requests +# import requests import logging import random import time import numpy as np import gzip +import sys from tqdm import tqdm from collections import OrderedDict @@ -26,9 +27,9 @@ from . import __version__ as read2tree_version from read2tree.stats.Coverage import Coverage from read2tree.stats.SeqCompleteness import SeqCompleteness -from read2tree.FastxReader import FastxReader +# from read2tree.FastxReader import FastxReader -API_URL = 'http://omabrowser.org/api' +# API_URL = 'http://omabrowser.org/api' class OGSet(object): @@ -86,8 +87,7 @@ def _reload_ogs_from_folder(self, folder_suffix='01_ref_ogs'): folders :return: Dictionary with og name as key and list of SeqRecords """ - print('--- Re-load ogs and find their corresponding DNA seq ' - 'from output folder ---') + print('--- Re-load ogs and find their corresponding DNA seq from output folder ---') ogs = {} ref_ogs_aa = sorted(glob.glob(os.path.join(os.path.join( self.args.output_path, folder_suffix+"_aa"), "*.fa"))) @@ -116,30 +116,37 @@ def _make_output_path(self, prefix): os.makedirs(path) return path - def _load_orthoxml(self): - if self.oma.mode == 'standalone': - og_orthoxml = os.path.join(self.oma_output_path, - 'OrthologousGroups.orthoxml') - tree_str = os.path.join(self.oma_output_path, - 'EstimatedSpeciesTree.nwk') - self._ham_analysis = pyham.Ham(tree_str, og_orthoxml, - use_internal_name=False) + # def _load_orthoxml(self): + # if self.oma.mode == 'standalone': + # og_orthoxml = os.path.join(self.oma_output_path, + # 'OrthologousGroups.orthoxml') + # tree_str = os.path.join(self.oma_output_path, + # 'EstimatedSpeciesTree.nwk') + # self._ham_analysis = pyham.Ham(tree_str, og_orthoxml, + # use_internal_name=False) def _load_dna_db(self): - if '.fa' in self.args.dna_reference or \ - '.fasta' in self.args.dna_reference: - db = {} - self.logger.info('--- Load ogs and find their corresponding ' - 'DNA seq from {} ---'.format(self.args.dna_reference)) - self.logger.info('Loading {} into memory. This might take a ' - 'while . . . '.format(self.args.dna_reference.split("/")[-1])) - - open_ = gzip.open if self.args.dna_reference.endswith('.gz') else open - with open_(self.args.dna_reference, 'rt') as f: - for rec in SeqIO.parse(f, 'fasta'): - db[rec.id.strip()] = str(rec.seq) - source = 'fa' - return db, source + # if not ('.fa' in self.args.dna_reference or \ + # '.fasta' in self.args.dna_reference): + # self.logger.info('--- The input ') + # sys.exit(0) + db = {} + self.logger.info('--- Load ogs and find their corresponding ' + 'DNA seq from {} ---'.format(self.args.dna_reference)) + self.logger.info('Loading {} into memory. This might take a ' + 'while . . . '.format(self.args.dna_reference.split("/")[-1])) + + open_ = gzip.open if self.args.dna_reference.endswith('.gz') else open + with open_(self.args.dna_reference, 'rt') as f: + for rec in SeqIO.parse(f, 'fasta'): + if "_" in rec.id: + self.logger.info( 'The fasta record id should not contain underscore _ . (You cannot use the output gene markers of read2tree in 01_ folders.)') + sys.exit(0) + db[rec.id.strip()] = str(rec.seq) + + source = 'fa' + return db, source + # ---------------- only to be used internally ---------------------- # elif '.h5' in self.args.dna_reference: # print('--- Load ogs and find their corresponding DNA \ @@ -151,11 +158,11 @@ def _load_dna_db(self): # .decode("utf-8") for row in # self._db_id_map.genome_table] # # print(self._db_species_list) - else: - self.logger.info('--- Load ogs and find their corresponding DNA seq using ' - 'the REST api ---') - source = 'REST_api' - return None, source + # else: + # self.logger.info('--- Load ogs and find their corresponding DNA seq using ' + # 'the REST api ---') + # source = 'REST_api' + # return None, source def _load_ogs(self): """ @@ -164,7 +171,7 @@ def _load_ogs(self): :return: Dictionary with og name as key and list of SeqRecords """ db, source = self._load_dna_db() - self._load_orthoxml() + #self._load_orthoxml() start = time.time() ogs = {} @@ -173,44 +180,48 @@ def _load_ogs(self): names_og = self.ogs - for name, records in tqdm(names_og.items(), desc='Loading OGs', - unit=' OGs'): + for name, records in tqdm(names_og.items(), desc='Loading OGs', unit=' OGs'): # name = file.split("/")[-1].split(".")[0] ogs[name] = OG() ogs[name].aa = self._get_aa_records(name, records) output_file_aa = os.path.join(orthologous_groups_aa, name + ".fa") output_file_dna = os.path.join(orthologous_groups_dna, name + ".fa") - - if source: - try: - ogs[name].dna = self._get_dna_records(ogs[name].aa, - db, source, name) - except (ValueError, TypeError): - self.logger.warning('This OG {} did not have any DNA'.format(name)) - pass - else: - all_len_consistent = self._check_dna_aa_length_consistency(name, ogs[name].aa, ogs[name].dna) - if "REST_api" in source and not all_len_consistent: - msg = "The returned DNA sequences from the REST API do not match the protein sequences. " \ - "Most likely this is due to an update of the OMA Browser. Please download the DNA " \ - "sequences from the download page of the OMA Browser that correspond to the release " \ - "of your reference groups." - self.logger.error(msg) - raise Exception(msg) - self._write(output_file_dna, ogs[name].dna) - self._write(output_file_aa, ogs[name].aa) - else: - self.logger.debug('DNA reference was not provided. ' - 'Only amino acid sequences gathered!') + ogs[name].dna = self._get_dna_records(ogs[name].aa, db, source, name) + + all_len_consistent = self._check_dna_aa_consistency(name, ogs[name].aa, ogs[name].dna) + + self._write(output_file_dna, ogs[name].dna) + self._write(output_file_aa, ogs[name].aa) + end = time.time() + self.elapsed_time = end - start + self.logger.info('{}: Gathering of DNA seq for {} OGs took {}.' + .format(self._species_name, len(names_og.keys()), self.elapsed_time)) + + # if source: + # try: + # ogs[name].dna = self._get_dna_records(ogs[name].aa, db, source, name) + # except (ValueError, TypeError): + # self.logger.warning('This OG {} did not have any DNA'.format(name)) + # pass + # else: + # all_len_consistent = self._check_dna_aa_length_consistency(name, ogs[name].aa, ogs[name].dna) + # if "REST_api" in source and not all_len_consistent: + # msg = "The returned DNA sequences from the REST API do not match the protein sequences. " \ + # "Most likely this is due to an update of the OMA Browser. Please download the DNA " \ + # "sequences from the download page of the OMA Browser that correspond to the release " \ + # "of your reference groups." + # self.logger.error(msg) + # raise Exception(msg) + # self._write(output_file_dna, ogs[name].dna) + # self._write(output_file_aa, ogs[name].aa) + # else: + # self.logger.debug('DNA reference was not provided. ' + # 'Only amino acid sequences gathered!') # self.progress.set_status('ogs') - end = time.time() - self.elapsed_time = end-start - self.logger.info('{}: Gathering of DNA seq for {} OGs took {}.' - .format(self._species_name, len(names_og.keys()), self.elapsed_time)) if db: db.clear() - if self._ham_analysis: - self._ham_analysis = None + # if self._ham_analysis: + # self._ham_analysis = None return ogs def _get_aa_records(self, name, records): @@ -318,30 +329,35 @@ def _get_dna_from_h5(self, record): def _get_dna_from_fasta(self, record, db): try: - if record.id.split("_")[0] not in db.keys(): - return self._get_dna_from_REST(record) - else: - dna = db[record.id.split("_")[0]] - except ValueError: - self.logger.debug('DNA not found for {}.'.format(record.id)) - pass - else: - return SeqRecord.SeqRecord(self._clean_DNA_seq(dna), - id=record.id, - description="") + # if record.id.split("_")[0] not in db.keys(): + # return self._get_dna_from_REST(record) + # else: + dna = db[record.id.split("_")[0]] + except: # ValueError + self.logger.debug('DNA not found for {} in the input dna gene marker fasta file. '.format(record.id)) + sys.exit(0) + + return SeqRecord.SeqRecord(self._clean_DNA_seq(dna), id=record.id, description="") + # else: # return SeqRecord.SeqRecord(Seq.Seq(dna.upper()), id=record.id, # description="") - def _check_dna_aa_length_consistency(self, og_name, aa, dna): + def _check_dna_aa_consistency(self, og_name, aa, dna): dna_dic = {r.id.split("_")[0]: r for r in dna} aa_dic = {r.id.split("_")[0]: r for r in aa} all_consistent = True for k, r_dna in dna_dic.items(): r_aa = aa_dic[k] if abs(len(r_dna.seq) - 3*len(r_aa.seq)) > 3: - self.logger.warning('{}: {} has aa-length {} and dna-length {}'.format(self._species_name, og_name+" "+k, 3*len(r_aa.seq), len(r_dna.seq))) + self.logger.error('{}: {} has aa-length {} and dna-length {}'.format(self._species_name, og_name+" "+k, 3*len(r_aa.seq), len(r_dna.seq))) all_consistent = False + if set(dna_dic.keys()) != set(aa_dic.keys()): + self.logger.error(' All aa marker genes should have corresponding dna sequences, not the case for {} '.format(str(set(aa_dic.keys()) - set(dna_dic.keys())))) + all_consistent = False + if not all_consistent: + sys.exit(0) + return all_consistent def _get_dna_records(self, records, db, source, og_name): @@ -351,18 +367,18 @@ def _get_dna_records(self, records, db, source, og_name): :return: """ og_cdna = [] - if 'REST_api' in source: - return self._get_dna_from_REST_bulk(records, og_name) - else: - for i, record in enumerate(records): - if 'h5' in source: - og_cdna.append(self._get_dna_from_h5(record)) - elif 'fa' in source: - og_cdna.append(self._get_dna_from_fasta(record, db)) - # elif 'REST_api' in source: - # og_cdna.append(self._get_dna_from_REST(record)) - - return og_cdna + # if 'REST_api' in source: + # return self._get_dna_from_REST_bulk(records, og_name) + # else: + for i, record in enumerate(records): + # if 'h5' in source: + # og_cdna.append(self._get_dna_from_h5(record)) + #elif 'fa' in source: + og_cdna.append(self._get_dna_from_fasta(record, db)) + # elif 'REST_api' in source: + # og_cdna.append(self._get_dna_from_REST(record)) + + return og_cdna def _clean_DNA_seq(self, record): """ @@ -563,7 +579,6 @@ def add_mapped_seq(self, mapper, species_name=None): def write_added_ogs_aa(self, folder_name=None): """ - :param self: :param folder_name: :return: @@ -582,7 +597,6 @@ def write_added_ogs_aa(self, folder_name=None): def write_added_ogs_dna(self, folder_name=None): """ - :param self: :param folder_name: :return: diff --git a/read2tree/__init__.py b/read2tree/__init__.py index c22e791..7c9df9f 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -6,10 +6,11 @@ from pkg_resources import resource_string logging.getLogger(__name__).addHandler(logging.NullHandler()) -__version__ = '1.5.0' +__version__ = '1.5.2' __copyright__ = 'read2tree (C) 2017-{:d} David Dylus ' \ .format(date.today().year) + # path = './log.yaml' # if os.path.exists(path): # with open(path, 'rt') as f: From 7ba1a41450f92203107668af730ea1202917828c Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 09:12:29 -0400 Subject: [PATCH 35/50] deactivating REST cleanup --- read2tree/OGSet.py | 1 - 1 file changed, 1 deletion(-) diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 0983bfe..1796a74 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -362,7 +362,6 @@ def _check_dna_aa_consistency(self, og_name, aa, dna): def _get_dna_records(self, records, db, source, og_name): """ - :param records: :return: """ From 6651b34c0e29dbeb1bc93096a2d49f4e1f563a72 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 09:12:51 -0400 Subject: [PATCH 36/50] deactivating REST cleanup --- read2tree/OGSet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 1796a74..7043dc4 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -268,7 +268,7 @@ def _get_dna_from_h5(self, record): # if 'X' in seq: cleaned_seq = self._clean_DNA_seq(seq) - # else: + # else: # a # cleaned_seq = seq return SeqRecord.SeqRecord(Seq.Seq(cleaned_seq), From f6d162bac0b72c2a380a8acf462a8fc260ebc36d Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 09:13:10 -0400 Subject: [PATCH 37/50] deactivating REST cleanup --- read2tree/OGSet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 7043dc4..65d4eff 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -268,7 +268,7 @@ def _get_dna_from_h5(self, record): # if 'X' in seq: cleaned_seq = self._clean_DNA_seq(seq) - # else: # a + # else: # # cleaned_seq = seq return SeqRecord.SeqRecord(Seq.Seq(cleaned_seq), From 64e7d07f50cb046fe8378b6219530d05b4ba8773 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 11:09:24 -0400 Subject: [PATCH 38/50] improve logging --- read2tree/Mapper.py | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index b3fc275..cfbb9c6 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -126,6 +126,7 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): reads_str = reads line_minimap= minimap2_ex +" "+ minimap_argm + " -t " + str(self.args.threads) +" "+ ref_file_handle + " " + reads_str + " > " + sam_file + #self._rm_file(ref_file_handle + "-enc.2.ngm", ignore_error=True) self._output_shell(line_minimap) self.logger.info('mapping with ' + line_minimap) @@ -352,6 +353,9 @@ def _map_reads_to_references(self, ref): else: mapped_reads = [] + else: + self.logger.error('No mapped reads.') + sys.exit(0) # self.progress.set_status('single_map', ref=species) #self._rm_file(ref_file_handle+".fai", ignore_error=True) @@ -606,7 +610,7 @@ def _post_process_read_mapping(self, ref_file, sam_file): # if self.args.single_mapping: # self.logger.warning("single_mapping is not tested in this version.") #debug("{}: ---- Samtools index completed".format(self._species_name)) - # self._rm_file(bam_file, ignore_error=True) + #f self.args.debug: # self._bin_reads(ref_file, sam_file_base ) @@ -649,6 +653,7 @@ def _rm_file(self, *fns, ignore_error=False): for fn in fns: try: os.remove(fn) + self.logger.info("we removed {}.".format(fn)) #debug("{}: ---- Samtools index completed".format(self._species_name)) except FileNotFoundError: if not ignore_error: raise @@ -670,23 +675,24 @@ def _output_shell(self, line): :param line: :return: """ - try: - self.logger.debug("Running " + line) - shell_command = subprocess.Popen( - line, stdout=subprocess.PIPE, stderr=subprocess.PIPE, - shell=True) - except: - self.logger.debug("Shell command failed to execute by running ") - return None + #try: + self.logger.debug("Running " + line) + shell_command = subprocess.Popen( line, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True) + + #except: + # self.logger.debug("Shell command failed to execute by running ") + # return None (output, err) = shell_command.communicate() self.logger.debug("Shell output: "+ str(output)) self.logger.debug("Shell err: " + str(err)) + print("Shell output: "+ str(output)) + print("Shell err: " + str(err)) shell_command.wait() if shell_command.returncode != 0: self.logger.debug("Shell command failed to execute") self.logger.debug(line) - return None + sys.exit(1) return output From 74d3872df0add9d714c3324f3afe4b57a3f36489 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 22 Oct 2024 12:53:18 -0400 Subject: [PATCH 39/50] improve logging --- archive/run_r2t.py | 4 ++++ read2tree/Mapper.py | 23 ++++++++++++++--------- read2tree/stats/Coverage.py | 2 +- read2tree/stats/SeqCompleteness.py | 2 +- 4 files changed, 20 insertions(+), 11 deletions(-) diff --git a/archive/run_r2t.py b/archive/run_r2t.py index 8ecfc7c..527f48c 100644 --- a/archive/run_r2t.py +++ b/archive/run_r2t.py @@ -9,6 +9,10 @@ print("start run_r2t 2 223 3 ") main(sys.argv[1:], exe_name=exe_name(), desc="descr") + +# --step 1marker --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --debug +# --step 2map --standalone_path marker_genes --dna_reference dna_ref.fa --reads /work/FAC/FBM/DBC/cdessim2/read2tree/v2_test/t1/reads_20/ERR7350657__2.fastq.gz --output_path output --debug --threads 1 + print("finish run_r2t ") diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index cfbb9c6..380693f 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -40,8 +40,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "minimap2" -samtools = "samtools" +minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2" +samtools = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools" class Mapper(object): @@ -354,8 +354,9 @@ def _map_reads_to_references(self, ref): mapped_reads = [] else: - self.logger.error('No mapped reads.') - sys.exit(0) + self.logger.warning('No mapped reads on {} .'.format(str(species))) + + #sys.exit(0) # for a few species , it is ok not to have mapped reads. # self.progress.set_status('single_map', ref=species) #self._rm_file(ref_file_handle+".fai", ignore_error=True) @@ -607,6 +608,8 @@ def _post_process_read_mapping(self, ref_file, sam_file): self._output_shell(line_samtools_index) self.logger.debug(line_samtools_index) + self._rm_file(sam_file_base + ".bam", ignore_error=True) + self._rm_file(sam_file_base + ".sam", ignore_error=True) # if self.args.single_mapping: # self.logger.warning("single_mapping is not tested in this version.") #debug("{}: ---- Samtools index completed".format(self._species_name)) @@ -653,7 +656,7 @@ def _rm_file(self, *fns, ignore_error=False): for fn in fns: try: os.remove(fn) - self.logger.info("we removed {}.".format(fn)) #debug("{}: ---- Samtools index completed".format(self._species_name)) + # self.logger.info("we removed {}.".format(fn)) #debug("{}: ---- Samtools index completed".format(self._species_name)) except FileNotFoundError: if not ignore_error: raise @@ -684,10 +687,12 @@ def _output_shell(self, line): # return None (output, err) = shell_command.communicate() - self.logger.debug("Shell output: "+ str(output)) - self.logger.debug("Shell err: " + str(err)) - print("Shell output: "+ str(output)) - print("Shell err: " + str(err)) + if output: + self.logger.debug("Shell output: "+ str(output)) + print("Shell output: " + str(output)) + if err: + self.logger.debug("Shell err: " + str(err)) + print("Shell err: " + str(err)) shell_command.wait() if shell_command.returncode != 0: self.logger.debug("Shell command failed to execute") diff --git a/read2tree/stats/Coverage.py b/read2tree/stats/Coverage.py index ed34723..398ffa1 100644 --- a/read2tree/stats/Coverage.py +++ b/read2tree/stats/Coverage.py @@ -32,7 +32,7 @@ def write_coverage_bam(self, file_name): gene_id = key.split("_")[0] coverage = value line = species + "," + og + "," + gene_id + "," + \ - str(coverage[0]) + "," + str(coverage[1]) + "\n" + str(round(coverage[0],2)) + "," + str(round(coverage[1],3)) + "\n" # coverage,std out_text += line with open(file_name, "w") as myfile: diff --git a/read2tree/stats/SeqCompleteness.py b/read2tree/stats/SeqCompleteness.py index 4922f07..ac3c95b 100644 --- a/read2tree/stats/SeqCompleteness.py +++ b/read2tree/stats/SeqCompleteness.py @@ -93,7 +93,7 @@ def write_seq_completeness(self, file_name): gene_id = key.split("_")[0] seq_completeness = value line = species + "," + og + "," + gene_id + "," + \ - str(seq_completeness[0]) + "," + str(seq_completeness[1]) + \ + str(round(seq_completeness[0],5)) + "," + str(round(seq_completeness[1],5)) + \ "," + str(seq_completeness[2]) + "," + \ str(seq_completeness[3]) + "," + \ str(seq_completeness[4]) + "\n" From d3a1b722f79d36e59341651467b187e63600b7d8 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Fri, 8 Nov 2024 10:20:04 -0500 Subject: [PATCH 40/50] improve logging --- read2tree/Mapper.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 380693f..64e9f05 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -40,8 +40,8 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader -minimap2_ex= "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/minimap2" -samtools = "/work/FAC/FBM/DBC/cdessim2/default/smajidi1/software/miniconda/envs/r2t_3.10.8b/bin/samtools" +minimap2_ex= "minimap2" +samtools = "samtools" class Mapper(object): @@ -361,7 +361,7 @@ def _map_reads_to_references(self, ref): #self._rm_file(ref_file_handle+".fai", ignore_error=True) #tmp_output_folder.cleanup() - end = time.time() + end = time.time() self.elapsed_time = end - start if len(references) > 1: self.logger.info('{}: Mapping to all references took {}.' From 43353a750620f1c35705b3b20c8aa96fa0e6e587 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Fri, 8 Nov 2024 10:22:34 -0500 Subject: [PATCH 41/50] update readme with minimap2 --- README.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index bca3185..7610e30 100644 --- a/README.md +++ b/README.md @@ -42,9 +42,10 @@ conda install -c conda-forge biopython numpy Cython ete3 lxml tqdm scipy pyparsi conda install -c bioconda dendropy pysam ``` -Besides, you need softwares including [mafft](http://mafft.cbrc.jp/alignment/software/) (multiple sequence aligner), [iqtree](http://www.iqtree.org/) (phylogenomic inference), [ngmlr](https://github.com/philres/ngmlr), [ngm/nextgenmap](https://github.com/Cibiv/NextGenMap) (long and short read mappers), and [samtools](http://www.htslib.org/download/) which could be installed using conda. +Besides, you need softwares including [mafft](http://mafft.cbrc.jp/alignment/software/) (multiple sequence aligner), [iqtree](http://www.iqtree.org/) (phylogenomic inference), [minimap2](https://github.com/lh3/minimap2) (long and short read mappers), and [samtools](http://www.htslib.org/download/) which could be installed using conda. +For this version, the `--read_type` argument could be either `short`, `long-hifi` or `long-ont`. You could also use `--threads 40` to be used with minimap2. ``` -conda install -c bioconda mafft iqtree ngmlr nextgenmap samtools +conda install -c bioconda mafft iqtree minimap2 samtools ``` Then, you can install the read2tree package after downlaoding the package from this GitHub repo using From 5963f54c531072711c166562b157c2d2d88c5c03 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Fri, 8 Nov 2024 10:24:35 -0500 Subject: [PATCH 42/50] update readme with minimap2 --- README.md | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 7610e30..9dd7740 100644 --- a/README.md +++ b/README.md @@ -68,8 +68,8 @@ conda install -c bioconda read2tree Alternatively, you could also try using [mamba](https://mamba.readthedocs.io/en/latest/). Caution: please read about compatiblity of conda and mamba in one envirnoment. ### 3) Installation using Docker +** not yet ready for this version of r2t (minimap2) The Dockerfile is also available in this repository. There is an example how to run in the [test example](#test-example) section. - A prebuild container can be loaded from dockerhub: ``` docker pull dessimozlab/read2tree:latest @@ -197,33 +197,30 @@ export LANG=en_US.UTF-8 ## Change log - +- version 1.5: + - using minimap2 as the read mapper - version 0.1.5: - fix issue with UnknownSeq being removed in Biopython>1.80 - removing unused modeltester wrappers - - version 0.1.4: - allow reference folders not named marker_genes (#12) - update environment.yml file to contain all dependencies (#16) - documentation improvements - CI/CD pipeline - - version 0.1.3: - improvements of documentation - adding support for docker - - small bugfixes - + - small bugfixes - version 0.1.2: packaging - - version 0.1.0: Adding covid analysis - - version 0.0: Initial work ## Authors -* [David Dylus](https://github.com/dvdylus), (main author) +* [David Dylus](https://github.com/dvdylus) * [Adrian Altenhoff](http://people.inf.ethz.ch/adriaal). +* [Sina Majidian](https://sinamajidian.github.io/) The authors would like to thank Alex Warwick for help how to initiate such a package. From 3664aab109ad53812dfcfb1dccf27194b0ac44af Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Sun, 10 Nov 2024 16:03:10 -0500 Subject: [PATCH 43/50] v1.5.3 --- read2tree/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/__init__.py b/read2tree/__init__.py index 7c9df9f..d2f7770 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -6,7 +6,7 @@ from pkg_resources import resource_string logging.getLogger(__name__).addHandler(logging.NullHandler()) -__version__ = '1.5.2' +__version__ = '1.5.3' __copyright__ = 'read2tree (C) 2017-{:d} David Dylus ' \ .format(date.today().year) From 8514bc1289e1d8b01a1d2d9f2e2e8a799790a959 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Mon, 18 Nov 2024 11:54:42 -0500 Subject: [PATCH 44/50] include reade when only one of paired-ends are mapped --- read2tree/Mapper.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index 64e9f05..a65ae90 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -539,14 +539,15 @@ def _build_consensus_seq_v2(self, ref_file, bam_file): :param bam_file: :return: """ - bam = pysam.AlignmentFile(bam_file) + bam = pysam.AlignmentFile(bam_file) # bam.header.to_dict() references = list(set([read.reference_name for read in bam.fetch()])) records = {rec.id: rec for rec in list(SeqIO.parse(ref_file, "fasta"))} new_records = {} for ref in references: # self.logger.info(read.qual) seq = list('N' * len(records[ref])) - for pileup_column in bam.pileup(ref, 0, 10000000): + for pileup_column in bam.pileup(ref, 0, 10000000, ignore_orphans=False): # assumption max length gene is 10m + # for read in bam.fetch(ref, 1, 1200000): print(read) # TODO: improve the selection of a column by its quality # qualities = [pileupread.alignment.query_alignment_qualities[pileupread.query_position] for pileupread in # pileupcolumn.pileups if not pileupread.is_del and not pileupread.is_refskip] From 7eca9dcfd14681272553781637202ef3172da402 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Mon, 17 Feb 2025 08:48:30 -0500 Subject: [PATCH 45/50] Update README.md --- README.md | 20 -------------------- 1 file changed, 20 deletions(-) diff --git a/README.md b/README.md index 9dd7740..8dbac47 100644 --- a/README.md +++ b/README.md @@ -57,26 +57,6 @@ python setup.py install ``` -### 2) Installation using Conda - - -``` -conda create -n r2t python=3.10.8 -conda install -c bioconda read2tree - -``` -Alternatively, you could also try using [mamba](https://mamba.readthedocs.io/en/latest/). Caution: please read about compatiblity of conda and mamba in one envirnoment. - -### 3) Installation using Docker -** not yet ready for this version of r2t (minimap2) -The Dockerfile is also available in this repository. There is an example how to run in the [test example](#test-example) section. -A prebuild container can be loaded from dockerhub: -``` -docker pull dessimozlab/read2tree:latest -``` - - - ## Run From 7356b64aea4a36df33d4740236969a92b117bb73 Mon Sep 17 00:00:00 2001 From: Sina Majidian Date: Tue, 18 Feb 2025 19:16:26 -0500 Subject: [PATCH 46/50] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 8dbac47..ea0c70e 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,7 @@ conda install -c bioconda mafft iqtree minimap2 samtools Then, you can install the read2tree package after downlaoding the package from this GitHub repo using ``` -git clone https://github.com/DessimozLab/read2tree.git +git clone https://github.com/DessimozLab/read2tree.git -b minimap2 cd read2tree python setup.py install ``` From 2da775e05161f64ba0471acf0f49b486b4a28d19 Mon Sep 17 00:00:00 2001 From: chalco777 Date: Mon, 4 Aug 2025 16:50:12 -0500 Subject: [PATCH 47/50] Edited Mapper to accept any minimap2 command-line option via --read_type. Updated README and help --- README.md | 4 ++-- read2tree/Mapper.py | 14 ++++++-------- read2tree/main.py | 6 +++--- 3 files changed, 11 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index ea0c70e..bed3451 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # read2tree -read2tree is a software tool that allows to obtain alignment matrices for tree inference. For this purpose it makes use of the OMA database and a set of reads. Its strength lies in the fact that it bipasses the several standard steps when obtaining such a matrix in regular analysis. These steps are read filtereing, assembly, gene prediction, gene annotation, all vs all comparison, orthology prediction, alignment and concatination. +read2tree is a software tool that allows to obtain alignment matrices for tree inference. For this purpose it makes use of the OMA database and a set of reads. Its strength lies in the fact that it bipasses the several standard steps when obtaining such a matrix in regular analysis. These steps are read filtereing, assembly, gene prediction, gene annotation, all vs all comparison, orthology prediction, alignment and concatenation. read2tree works in linux with [![Python 3.10.8](https://img.shields.io/badge/python-3.10.8-blue.svg)](https://www.python.org/downloads/release/python-310/) @@ -43,7 +43,7 @@ conda install -c bioconda dendropy pysam ``` Besides, you need softwares including [mafft](http://mafft.cbrc.jp/alignment/software/) (multiple sequence aligner), [iqtree](http://www.iqtree.org/) (phylogenomic inference), [minimap2](https://github.com/lh3/minimap2) (long and short read mappers), and [samtools](http://www.htslib.org/download/) which could be installed using conda. -For this version, the `--read_type` argument could be either `short`, `long-hifi` or `long-ont`. You could also use `--threads 40` to be used with minimap2. +For this version, the `--read_type` argument accepts any minimap2 options string that defines how reads are aligned to the reference. For example, it could be `-ax sr`, `-ax map-hifi` or `-ax map-ont`. You can also pass `--threads 40` to be used with minimap2. ``` conda install -c bioconda mafft iqtree minimap2 samtools ``` diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index a65ae90..82981d4 100644 --- a/read2tree/Mapper.py +++ b/read2tree/Mapper.py @@ -112,14 +112,12 @@ def _call_wrapper(self, ref_file_handle, reads, tmp_output_folder): sam_file= output_folder+"/"+ref_file_handle.split('/')[-1].split('.')[0]+".sam" #self._output_shell(minimap2_ex+" -ax sr "+ ref_file_handle+ " -t " + str(self.args.threads) + " "+reads+" | "+samtools_ex+" view -F 4 -bh -S -t" + str(self.args.threads)+ " > " + bam_file) - if 'short' in self.args.read_type: - minimap_argm =" -ax sr" - elif 'long-hifi' in self.args.read_type: - minimap_argm = " -ax map-hifi " - elif 'long-ont' in self.args.read_type: - minimap_argm = " -ax map-ont " - else: - self.logger.error(" read_type is not valid" +self.args.read_type+" please use one of these: short long-ont long-hifi ") + + rt = (self.args.read_type or "").strip() + if not rt: + self.logger.error("read_type is not set, please use one of these: -ax short, -ax long-ont, -ax long-hifi") + minimap_argm = f" {rt} " + if isinstance(reads,list): reads_str=" ".join(reads) elif isinstance(reads,str): diff --git a/read2tree/main.py b/read2tree/main.py index 6751b31..c6773b4 100644 --- a/read2tree/main.py +++ b/read2tree/main.py @@ -73,9 +73,9 @@ def parse_args(argv, exe_name, desc): help='[Default is none] Reads to be mapped to reference. If paired ' 'end add separated by space.') - arg_parser.add_argument('--read_type', default='short', - help='[Default is short] Type of reads to ' - 'use for mapping: short, long-hifi or long-ont corresponding to sr, map-hifi, or map-ont in minimap2.') + arg_parser.add_argument('--read_type', default='-ax sr', + help='[Default is -ax sr] Minimap2 command-line options for mapping reads to reference. ' + 'Examples: -ax sr , -ax map-hifi , -ax map-pb or -ax map-ont ') arg_parser.add_argument('--threads', type=int, default=1, help='[Default is 1] Number of threads for the mapping ') From f68f626078088255ec838dab33ec847ff3340569 Mon Sep 17 00:00:00 2001 From: chalco777 Date: Wed, 13 Aug 2025 18:20:08 -0500 Subject: [PATCH 48/50] Commented unnecesary import of pyham --- read2tree/OGSet.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 65d4eff..4f62a3c 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -7,7 +7,7 @@ import glob import os import re -import pyham +# import pyham # import requests import logging import random From bba89dee18d857b18c5524813e3572d56eb60a88 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Adri=C3=A1n=20Chalco?= <160282720+chalco777@users.noreply.github.com> Date: Tue, 2 Sep 2025 12:36:21 -0500 Subject: [PATCH 49/50] Update setup.py --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 4c45fa4..ff6f246 100644 --- a/setup.py +++ b/setup.py @@ -12,7 +12,7 @@ # conda install -c bioconda dendropy requirements = ["numpy", "biopython", "ete3", "lxml", "tqdm", "scipy", "pyparsing", "requests", "natsort", "pyyaml", "dendropy", - "pysam", "pyham", "filelock"] + "pysam", "filelock"] with open("README.md", "r", encoding="utf-8") as fh: long_description = fh.read() From 6d935bcb2eba70b70b31dc4bc16c1c5f54f192bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Adri=C3=A1n=20Chalco?= <160282720+chalco777@users.noreply.github.com> Date: Tue, 2 Sep 2025 12:36:56 -0500 Subject: [PATCH 50/50] Update environment.yml --- environment.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/environment.yml b/environment.yml index 8e61bbc..d4eb262 100644 --- a/environment.yml +++ b/environment.yml @@ -22,5 +22,4 @@ dependencies: - nextgenmap - samtools - filelock - - pyham - - pysam \ No newline at end of file + - pysam