diff --git a/README.md b/README.md index acd7aafb..06ecd1e7 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/) @@ -47,40 +47,21 @@ 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 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 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 ``` -git clone https://github.com/DessimozLab/read2tree.git +git clone https://github.com/DessimozLab/read2tree.git -b minimap2 cd read2tree 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 -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 @@ -100,7 +81,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 @@ -109,12 +90,23 @@ read2tree --tree --standalone_path marker_genes/ --reads read_1.fastq read_2.fas ``` ### Multiple species mode + +#### step1 ``` -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 1marker --standalone_path marker_genes --dna_reference dna_ref.fa --output_path output --debug +``` + +#### step2 +The following could be run in parallel. +``` +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 @@ -141,12 +133,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 sample_1.fastq sample_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 ``` @@ -190,33 +182,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. diff --git a/read2tree/Progress.py b/archive/Progress.py similarity index 75% rename from read2tree/Progress.py rename to archive/Progress.py index 85c3406d..cfadf22d 100644 --- a/read2tree/Progress.py +++ b/archive/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 diff --git a/archive/run_r2t.py b/archive/run_r2t.py new file mode 100644 index 00000000..527f48ce --- /dev/null +++ b/archive/run_r2t.py @@ -0,0 +1,20 @@ + + +import read2tree +from read2tree.main import main +from read2tree._utils import exe_name + +import sys + +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 ") + + + +a=1 \ No newline at end of file diff --git a/archive/tests/test_use.py b/archive/tests/test_use.py index 0d87f799..e6721290 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__": diff --git a/environment.yml b/environment.yml index 8e61bbc1..d4eb262d 100644 --- a/environment.yml +++ b/environment.yml @@ -22,5 +22,4 @@ dependencies: - nextgenmap - samtools - filelock - - pyham - - pysam \ No newline at end of file + - pysam diff --git a/read2tree/Aligner.py b/read2tree/Aligner.py index f85edf5d..dbee3fbb 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: - print('--- Alignment of {} OGs ---'.format(len(list(og_set.keys())))) + if step== "all" or step== "1marker":#and og_set is not None: + logger.info('--- 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)) @@ -149,36 +149,38 @@ 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 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, diff --git a/read2tree/Mapper.py b/read2tree/Mapper.py index ea7b3725..4a47482d 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: @@ -39,13 +40,16 @@ from read2tree.stats.SeqCompleteness import SeqCompleteness from read2tree.FastxReader import FastxReader +minimap2_ex= "minimap2" +samtools = "samtools" + 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 @@ -68,74 +72,99 @@ 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 ref_set != None: - self.mapped_records = \ - self._map_reads_to_references(ref_set) + if step=="2map" or step=="all": # compute mapping + 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() - else: # re-load already computed mapping - if og_set != None and not self.args.merge_all_mappings: + elif step =="3combine": # re-load already computed mapping + 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 - 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' + # 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) - - 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) + 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) + + 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): + 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) + + # 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)) + 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: - 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) + 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 @@ -147,22 +176,24 @@ 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') + 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 = [] @@ -174,20 +205,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) @@ -230,16 +257,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): @@ -266,11 +294,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, @@ -281,9 +311,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, @@ -320,16 +351,25 @@ def _map_reads_to_references(self, ref): else: mapped_reads = [] + else: + 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) + #self._rm_file(ref_file_handle+".fai", ignore_error=True) - tmp_output_folder.cleanup() - end = time.time() + #tmp_output_folder.cleanup() + end = time.time() self.elapsed_time = end - start if len(references) > 1: self.logger.info('{}: Mapping to all references took {}.' .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? 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) + return mapped_reads_species def _write_read_query_aling(self, read, og_name_file, write_mode): @@ -414,29 +454,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") @@ -495,14 +537,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] @@ -516,7 +559,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 @@ -525,50 +568,55 @@ 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") - if self.args.single_mapping: - self.logger.debug("{}: --- POSTPROCESSING MAPPING " - "---".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") + # 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]: - 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 self.args.single_mapping: - self.logger.debug("{}: ---- Samtools view completed" - .format(self._species_name)) - - if os.path.exists(bam_file): - self._output_shell( - 'samtools sort -m 2G -@ ' + str(self.args.threads) + - ' -o ' + outfile_name + "_sorted.bam " + bam_file) - if self.args.single_mapping: - self.logger.debug("{}: ---- Samtools sort completed" - .format(self._species_name)) - - 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._rm_file(bam_file, ignore_error=True) - if self.args.debug: - self._bin_reads(ref_file, outfile_name + '_sorted.bam') - - consensus = self._build_consensus_seq_v2(ref_file, outfile_name + - '_sorted.bam') + #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.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 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) + + 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)) + + + #f self.args.debug: + # self._bin_reads(ref_file, sam_file_base ) + + consensus = self._build_consensus_seq_v2(ref_file, sam_file_base +'_sorted.bam') all_consensus = [] if consensus: @@ -596,7 +644,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) @@ -607,19 +655,20 @@ 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 - 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): """ @@ -628,26 +677,38 @@ 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 + +# 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() 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)) +# if err: +# self.logger.debug("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 diff --git a/read2tree/OGSet.py b/read2tree/OGSet.py index 370558bb..4f62a3c0 100644 --- a/read2tree/OGSet.py +++ b/read2tree/OGSet.py @@ -7,13 +7,14 @@ import glob import os import re -import pyham -import requests +# import pyham +# 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,14 +27,14 @@ 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): - 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 +50,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 +70,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 @@ -84,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"))) @@ -114,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 \ @@ -149,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): """ @@ -162,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 = {} @@ -171,46 +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!') + output_file_aa = os.path.join(orthologous_groups_aa, name + ".fa") + output_file_dna = os.path.join(orthologous_groups_dna, name + ".fa") + 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): @@ -257,112 +268,116 @@ def _get_dna_from_h5(self, record): # if 'X' in seq: cleaned_seq = self._clean_DNA_seq(seq) - # else: + # else: # # cleaned_seq = seq return SeqRecord.SeqRecord(Seq.Seq(cleaned_seq), id=record.id, description="") - def _get_dna_from_REST(self, record): - """ + # 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 - :param record: - :param name: - :return: - """ - tmp_id = re.sub(r'\..*', '', record.id.split("_")[0]) - use_id = re.sub(r'\W+', '', tmp_id) + def _get_dna_from_fasta(self, record, db): 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) + # if record.id.split("_")[0] not in db.keys(): + # return self._get_dna_from_REST(record) # else: - # cleaned_seq = dna_record.seq - return SeqRecord.SeqRecord(cleaned_seq, record.id, - description="", name="") + 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) - def _get_dna_from_REST_bulk(self, records, og_name): - """ + return SeqRecord.SeqRecord(self._clean_DNA_seq(dna), id=record.id, description="") - :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: - 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="") # 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): """ - :param records: :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): """ @@ -451,7 +466,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: @@ -563,7 +578,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 +596,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/Reads.py b/read2tree/Reads.py index 3848a0b2..22bc664d 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 30882a45..28174214 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 8e3e4b20..d2f77707 100644 --- a/read2tree/__init__.py +++ b/read2tree/__init__.py @@ -6,19 +6,20 @@ from pkg_resources import resource_string logging.getLogger(__name__).addHandler(logging.NullHandler()) -__version__ = '0.1.5' -__copyright__ = 'read2tree (C) 2017-{:d} David Dylus' \ +__version__ = '1.5.3' +__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: # 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 diff --git a/read2tree/main.py b/read2tree/main.py index 23b33b01..a14a6fa0 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,54 +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 or long. Either ngm for short reads or ' - 'ngmlr for long will be used.') + 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 ' - 'using ngm / ngmlr!') - - 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_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. ') + 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_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_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.') @@ -103,37 +111,36 @@ 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.') - - 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('--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', + '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 ' + # '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', + # 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', # 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') @@ -150,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 1marker 2map 3combine ') + + # 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 = "" @@ -222,27 +228,27 @@ def parse_args(argv, exe_name, desc): if args.species_name: _species_name = args.species_name - if args.merge_all_mappings: + if args.step == "3combine": # 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( - '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): @@ -250,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) - 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. @@ -265,171 +295,141 @@ 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 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) - 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 read2tree 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 == "1marker": + logger.error( + "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) - 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 == "2map": + 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": + 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, 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. # 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) # 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 # write 03_ + # 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 == "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, 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() ogset.add_mapped_seq(mapper) ogset.write_added_ogs_aa() ogset.write_added_ogs_dna() - alignments.add_mapped_seq(ogset.mapped_ogs) + # alignments.add_mapped_seq(ogset.mapped_ogs) + # alignments.write_added_align_aa() + # alignments.write_added_align_dna() + + print("done- 2map") + logger.info(' ------- Read2Tree step 2map finished -*- -------') + + 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 + alignments.remove_species_from_alignments() + ogset.remove_species_from_ogs() + 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, # 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") + 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- 3combine") + logger.info(' ------- Read2Tree step 3combine finished -*- -------') + + logger.info(' ------- Read2Tree finished -*- -------') + + print("done- main ") + + # 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 diff --git a/read2tree/parser/OMAOutputParser.py b/read2tree/parser/OMAOutputParser.py index 90cba481..50488f6f 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): diff --git a/read2tree/stats/Coverage.py b/read2tree/stats/Coverage.py index e0ac291c..398ffa1a 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] @@ -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 4922f076..ac3c95b4 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" diff --git a/setup.py b/setup.py index 4c45fa4e..ff6f246b 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() diff --git a/tests/test_use.py b/tests/test_use.py index 0d87f799..e6721290 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__":