From d6b43c4ad26419cffb89ec3eeb2ae9d35cd01cbe Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Sun, 2 Jun 2024 23:57:10 +0200 Subject: [PATCH 01/10] Add coverage and identity constraints for GMAP --- pepti_map/aligning/gmap_wrapper.py | 34 +++++++++++++++++++++++++++++- pepti_map/main.py | 33 +++++++++++++++++++++++++++-- 2 files changed, 64 insertions(+), 3 deletions(-) diff --git a/pepti_map/aligning/gmap_wrapper.py b/pepti_map/aligning/gmap_wrapper.py index 9166bce..07b19f3 100644 --- a/pepti_map/aligning/gmap_wrapper.py +++ b/pepti_map/aligning/gmap_wrapper.py @@ -8,7 +8,7 @@ class GmapWrapper: - def __init__(self): + def __init__(self, min_trimmed_coverage: float = 0.0, min_identity: float = 0.0): env_vars = dotenv_values() try: n_threads = env_vars.get("GMAP_N_THREADS") @@ -27,6 +27,30 @@ def __init__(self): self._n_threads = n_threads self._batch_mode = batch_mode + try: + assert min_trimmed_coverage >= 0.0 and min_trimmed_coverage <= 1.0 + self._min_trimmed_coverage = min_trimmed_coverage + except AssertionError: + logging.info( + ( + "--min_trimmed_coverage option for GMAP must be between 0.0 and " + f"1.0, but was {str(min_trimmed_coverage)}. Using default of 0.0" + ) + ) + self._min_trimmed_coverage = 0.0 + + try: + assert min_identity >= 0.0 and min_identity <= 1.0 + self._min_identity = min_identity + except AssertionError: + logging.info( + ( + "--min_identity option for GMAP must be between 0.0 and " + f"1.0, but was {str(min_identity)}. Using default of 0.0" + ) + ) + self._min_identity = 0.0 + def build_index( self, files_to_index: List[Path], @@ -111,6 +135,10 @@ def produce_alignment( str(self._n_threads), "-f", "gff3_gene", + "--min-trimmed-coverage", + str(self._min_trimmed_coverage), + "--min-identity", + str(self._min_identity), ] alignment_command.extend( [ @@ -137,6 +165,10 @@ def _produce_alignment_for_single_file(self, path_to_sequences: Path) -> None: str(self._n_threads), "-f", "gff3_gene", + "--min-trimmed-coverage", + str(self._min_trimmed_coverage), + "--min-identity", + str(self._min_identity), path_to_sequences.absolute().as_posix(), ] with open( diff --git a/pepti_map/main.py b/pepti_map/main.py index 38d3d7d..3b3941f 100644 --- a/pepti_map/main.py +++ b/pepti_map/main.py @@ -205,8 +205,10 @@ def align_reads_to_genome( genome: Union[str, None], gmap_index: Union[str, None], output_dir: str, + min_trimmed_coverage: float, + min_identity: float, ) -> None: - gmap_wrapper = GmapWrapper() + gmap_wrapper = GmapWrapper(min_trimmed_coverage, min_identity) # TODO: How to automatically use previously generated index? if gmap_index is not None and gmap_index != "": gmap_index_path = Path(gmap_index) @@ -394,6 +396,24 @@ def concat_output(paths_to_subdirectories: List[Path], output_dir: str) -> None: "the '-g/--genome' option is ignored." ), ) +@click.option( + "-mtc", + "--min-trimmed-coverage", + required=False, + type=float, + default=0.0, + show_default=True, + help="Sets the '--min-trimmed-coverage' option for GMAP during alignment.", +) +@click.option( + "-mid", + "--min-identity", + required=False, + type=float, + default=0.0, + show_default=True, + help="Sets the '--min-identity' option for GMAP during alignment.", +) def main( peptide_file: str, rna_file: str, @@ -407,6 +427,8 @@ def main( min_contig_length: int, genome: Union[str, None], gmap_index: Union[str, None], + min_trimmed_coverage: float, + min_identity: float, ): _setup(output_dir) @@ -468,7 +490,14 @@ def main( if last_step < Step.ALIGNMENT.value: logging.info("Aligning assembled RNA-seq reads to the genome.") - align_reads_to_genome(trinity_results_paths, genome, gmap_index, output_dir) + align_reads_to_genome( + trinity_results_paths, + genome, + gmap_index, + output_dir, + min_trimmed_coverage, + min_identity, + ) else: logging.info("Using already generated alignments.") From 37f95ee0af9ff984ca46ee25a7126a77bb0e05ba Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Sun, 23 Jun 2024 23:37:56 +0200 Subject: [PATCH 02/10] Fix error when no GMAP alignments --- pepti_map/main.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/pepti_map/main.py b/pepti_map/main.py index 3b3941f..8418cd4 100644 --- a/pepti_map/main.py +++ b/pepti_map/main.py @@ -196,7 +196,7 @@ def generate_trinity_results( return trinity_results_paths -def load_trinity_results_paths() -> List[Path]: +def load_current_results_paths() -> List[Path]: return TrinityWrapper.load_results_filepaths() @@ -207,7 +207,7 @@ def align_reads_to_genome( output_dir: str, min_trimmed_coverage: float, min_identity: float, -) -> None: +) -> List[Path]: gmap_wrapper = GmapWrapper(min_trimmed_coverage, min_identity) # TODO: How to automatically use previously generated index? if gmap_index is not None and gmap_index != "": @@ -222,13 +222,30 @@ def align_reads_to_genome( logging.error(missing_option_message) raise ValueError(missing_option_message) + new_results_paths: List[Path] = [] for trinity_results_path in trinity_results_paths: gmap_wrapper.produce_alignment( [trinity_results_path], trinity_results_path.parent / "alignment_result.gff3", ) + # Check if actual output was produced + with open( + trinity_results_path.parent / "alignment_result.gff3", + "rt", + encoding="utf-8", + ) as gmap_output: + line_count = 0 + for _ in gmap_output: + line_count += 1 + if line_count == 4: + new_results_paths.append(trinity_results_path) + break + + # Save new output paths + TrinityWrapper.save_results_filepaths(new_results_paths) _write_last_step(Step.ALIGNMENT.value) logging.info("Generated alignment of assembled contigs with GMAP.") + return new_results_paths def generate_pogo_input(paths_to_subdirectories: List[Path], peptide_file: str) -> None: @@ -486,11 +503,11 @@ def main( ) else: logging.info("Using already generated Trinity output files.") - trinity_results_paths = load_trinity_results_paths() + trinity_results_paths = load_current_results_paths() if last_step < Step.ALIGNMENT.value: logging.info("Aligning assembled RNA-seq reads to the genome.") - align_reads_to_genome( + trinity_results_paths = align_reads_to_genome( trinity_results_paths, genome, gmap_index, @@ -500,6 +517,7 @@ def main( ) else: logging.info("Using already generated alignments.") + trinity_results_paths = load_current_results_paths() paths_to_subdirectories = [ trinity_results_path.parent for trinity_results_path in trinity_results_paths From d87a687a7d9d0cc780e90297dbc58dcd901aa69c Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Mon, 24 Jun 2024 00:04:52 +0200 Subject: [PATCH 03/10] Update gmap for mac os --- environment_macos.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/environment_macos.yml b/environment_macos.yml index c5f0df6..1da2d67 100644 --- a/environment_macos.yml +++ b/environment_macos.yml @@ -52,10 +52,10 @@ dependencies: - bwidget=1.9.14 - bzip2=1.0.8 - c-ares=1.20.1 - - ca-certificates=2023.11.17 + - ca-certificates=2024.6.2 - cairo=1.16.0 - cctools_osx-64=973.0.1 - - certifi=2023.11.17 + - certifi=2024.6.2 - charset-normalizer=3.3.1 - clang=16.0.1 - clang-16=16.0.1 @@ -87,7 +87,7 @@ dependencies: - gfortran_impl_osx-64=11.4.0 - gfortran_osx-64=11.4.0 - glpk=5.0 - - gmap=2023.10.10 + - gmap=2024.05.20 - gmp=6.2.1 - graphite2=1.3.13 - gsl=2.7 @@ -157,7 +157,7 @@ dependencies: - numpy=1.24.4 - openjdk=21.0.1 - openjpeg=2.5.0 - - openssl=3.2.0 + - openssl=3.3.1 - packaging=23.2 - pandas=2.0.3 - pandoc=3.1.3 From f3be503b2550e594156d0cc0302d3c36b3bbcd9a Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Tue, 2 Jul 2024 22:27:07 +0200 Subject: [PATCH 04/10] Update gmap for linux --- environment_linux.yml | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/environment_linux.yml b/environment_linux.yml index 90022ee..f8d61c8 100644 --- a/environment_linux.yml +++ b/environment_linux.yml @@ -56,9 +56,9 @@ dependencies: - bwidget=1.9.14 - bzip2=1.0.8 - c-ares=1.24.0 - - ca-certificates=2023.11.17 + - ca-certificates=2024.6.2 - cairo=1.16.0 - - certifi=2023.11.17 + - certifi=2024.6.2 - charset-normalizer=3.3.2 - click=8.1.7 - contourpy=1.1.1 @@ -82,7 +82,7 @@ dependencies: - gfortran_impl_linux-64=13.2.0 - giflib=5.2.1 - glpk=5.0 - - gmap=2023.10.10 + - gmap=2024.05.20 - gmp=6.3.0 - graphite2=1.3.13 - gsl=2.7 @@ -155,7 +155,7 @@ dependencies: - numpy=1.24.4 - openjdk=17.0.3 - openjpeg=2.5.0 - - openssl=3.2.0 + - openssl=3.3.1 - packaging=23.2 - pandas=2.0.3 - pandoc=3.1.3 From ac20ac149585b893f0e9ec932bb627a70b705190 Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Mon, 19 Aug 2024 21:55:12 +0200 Subject: [PATCH 05/10] Handle antisense direction during gtf creation --- .../output_generation/pogo_input_helper.py | 24 +++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index 93aca3e..c2191b5 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -132,6 +132,8 @@ def _write_new_feature_coordinates( ): raise ValueError("No start and end coordinates given.") + # Adapt attributes for start/stop of coverage in sequence to the whole length of + # the contig even if not the whole contig was matched to the genome if direction == ".": contig_id, start_exon_contig_end, contig_start, _ = start_exon.attributes[ "Target" @@ -179,6 +181,27 @@ def _write_new_feature_coordinates( cds_ids = [cds_entry.attributes["ID"][0] for cds_entry in cds] mrna_id = mrna.attributes["ID"][0] + # If direction is antisense, change alignment to be on the complementary strand + if direction == "-": + if strand == "+": + new_strand = "-" + else: + new_strand = "+" + + gene.strand = new_strand + mrna.strand = new_strand + for exon in exons: + exon.strand = new_strand + for cds_entry in cds: + cds_entry.strand = new_strand + + exons.reverse() + cds.reverse() + + start_exon_index = (len(exons) - 1) - start_exon_index + end_exon_index = (len(exons) - 1) - end_exon_index + + # Adapt the alignment genome coordinates to fit the whole contig length if ( (direction == "+" and strand == "+") or (direction == "-" and strand == "-") @@ -295,6 +318,7 @@ def generate_gtf_input_file( path_to_gff.absolute().as_posix(), (path_to_gff.parent / "gffutils_db.sqlite").absolute().as_posix(), ) + with open( output_directory / "pogo_gtf_in.gtf", "wt", encoding="utf-8" ) as output_gtf: From 7b9e8f865bcfb96b7c1bc456e9acc4c542b4b5ec Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Fri, 30 Aug 2024 01:35:38 +0200 Subject: [PATCH 06/10] WIP: Cut out contig parts based on alignment --- .../output_generation/pogo_input_helper.py | 514 ++++++++++-------- 1 file changed, 292 insertions(+), 222 deletions(-) diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index c2191b5..9310e94 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -108,6 +108,200 @@ def _write_feature_in_gtf_format( + "\n" ) + # @classmethod + # def _write_new_feature_coordinates_old( + # cls, + # output_gtf: TextIO, + # gene: gffutils.Feature, + # mrna: gffutils.Feature, + # exons: List[gffutils.Feature], + # cds: List[gffutils.Feature], + # start_exon_index: int, + # end_exon_index: int, + # strand: str, + # direction: str, + # contig_length: int, + # ) -> None: + # start_exon = exons[start_exon_index] + # end_exon = exons[end_exon_index] + # if ( + # not start_exon.start + # or not start_exon.end + # or not end_exon.start + # or not end_exon.end + # ): + # raise ValueError("No start and end coordinates given.") + + # # Adapt attributes for start/stop of coverage in sequence to the whole length of + # # the contig even if not the whole contig was matched to the genome + # if direction == ".": + # contig_id, start_exon_contig_end, contig_start, _ = start_exon.attributes[ + # "Target" + # ][0].split(" ") + # _, contig_end, end_exon_contig_start, _ = end_exon.attributes["Target"][ + # 0 + # ].split(" ") + # contig_start = int(contig_start) + # contig_end = int(contig_end) + + # if start_exon == end_exon: + # start_exon.attributes["Target"] = " ".join( + # [contig_id, str(contig_length), "1", direction] + # ) + # else: + # start_exon.attributes["Target"] = " ".join( + # [contig_id, start_exon_contig_end, "1", direction] + # ) + # end_exon.attributes["Target"] = " ".join( + # [contig_id, str(contig_length), end_exon_contig_start, direction] + # ) + # else: + # contig_id, contig_start, start_exon_contig_end, _ = start_exon.attributes[ + # "Target" + # ][0].split(" ") + # _, end_exon_contig_start, contig_end, _ = end_exon.attributes["Target"][ + # 0 + # ].split(" ") + # contig_start = int(contig_start) + # contig_end = int(contig_end) + + # if start_exon == end_exon: + # start_exon.attributes["Target"] = " ".join( + # [contig_id, "1", str(contig_length), direction] + # ) + # else: + # start_exon.attributes["Target"] = " ".join( + # [contig_id, "1", start_exon_contig_end, direction] + # ) + # end_exon.attributes["Target"] = " ".join( + # [contig_id, end_exon_contig_start, str(contig_length), direction] + # ) + + # exon_ids = [exon.attributes["ID"][0] for exon in exons] + # cds_ids = [cds_entry.attributes["ID"][0] for cds_entry in cds] + # mrna_id = mrna.attributes["ID"][0] + + # # If direction is antisense, change alignment to be on the complementary strand + # if direction == "-": + # if strand == "+": + # new_strand = "-" + # else: + # new_strand = "+" + + # gene.strand = new_strand + # mrna.strand = new_strand + # for exon in exons: + # exon.strand = new_strand + # for cds_entry in cds: + # cds_entry.strand = new_strand + + # exons.reverse() + # cds.reverse() + + # start_exon_index = (len(exons) - 1) - start_exon_index + # end_exon_index = (len(exons) - 1) - end_exon_index + + # # Adapt the alignment genome coordinates to fit the whole contig length + # if ( + # (direction == "+" and strand == "+") + # or (direction == "-" and strand == "-") + # or (direction == "." and strand == "+") + # ): + # new_start = start_exon.start - (contig_start - 1) + # start_exon.start = new_start + # mrna.start = new_start + # gene.start = new_start + # new_end = end_exon.end + (contig_length - contig_end) + # end_exon.end = new_end + # mrna.end = new_end + # gene.end = new_end + + # cls._write_feature_in_gtf_format( + # output_gtf, + # gene, + # gene.attributes["ID"][0], + # ) + # for frame in range(3): + # mrna.attributes["ID"] = mrna_id + "." + str(frame) + # cls._write_feature_in_gtf_format( + # output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] + # ) + # for exon_index, exon in enumerate(exons): + # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) + # exon.attributes["Parent"] = mrna.attributes["ID"] + # cls._write_feature_in_gtf_format( + # output_gtf, + # exon, + # gene.attributes["ID"][0], + # mrna.attributes["ID"][0], + # ) + # for cds_index, cds_entry in enumerate(cds): + # cds_entry.start = exons[cds_index].start + # cds_entry.end = exons[cds_index].end + # if cds_index == start_exon_index: + # cds_entry.start = cds_entry.start + frame # pyright: ignore + # if cds_index == end_exon_index: + # cds_entry.end = cds_entry.end - ( # pyright: ignore + # (contig_length - frame) % 3 + # ) + # cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) + # cds_entry.attributes["Parent"] = mrna.attributes["ID"] + # cls._write_feature_in_gtf_format( + # output_gtf, + # cds_entry, + # gene.attributes["ID"][0], + # mrna.attributes["ID"][0], + # ) + + # elif ( + # (direction == "+" and strand == "-") + # or (direction == "-" and strand == "+") + # or (direction == "." and strand == "-") + # ): + # new_end = start_exon.end + (contig_start - 1) + # start_exon.end = new_end + # mrna.end = new_end + # gene.end = new_end + # new_start = end_exon.start - (contig_length - contig_end) + # end_exon.start = new_start + # mrna.start = new_start + # gene.start = new_start + + # cls._write_feature_in_gtf_format(output_gtf, gene, gene.attributes["ID"][0]) + # for frame in range(3): + # mrna.attributes["ID"] = mrna_id + "." + str(frame) + # cls._write_feature_in_gtf_format( + # output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] + # ) + # for exon_index, exon in enumerate(exons): + # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) + # exon.attributes["Parent"] = mrna.attributes["ID"] + # cls._write_feature_in_gtf_format( + # output_gtf, + # exon, + # gene.attributes["ID"][0], + # mrna.attributes["ID"][0], + # ) + # for cds_index, cds_entry in enumerate(cds): + # cds_entry.start = exons[cds_index].start + # cds_entry.end = exons[cds_index].end + # if cds_index == start_exon_index: + # cds_entry.end = cds_entry.end - frame # pyright: ignore + # if cds_index == end_exon_index: + # cds_entry.start = cds_entry.start + ( # pyright: ignore + # (contig_length - frame) % 3 + # ) + # cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) + # cds_entry.attributes["Parent"] = mrna.attributes["ID"] + # cls._write_feature_in_gtf_format( + # output_gtf, + # cds_entry, + # gene.attributes["ID"][0], + # mrna.attributes["ID"][0], + # ) + # else: + # raise ValueError("Strand must be one of '+', '-'.") + @classmethod def _write_new_feature_coordinates( cls, @@ -115,72 +309,11 @@ def _write_new_feature_coordinates( gene: gffutils.Feature, mrna: gffutils.Feature, exons: List[gffutils.Feature], - cds: List[gffutils.Feature], - start_exon_index: int, - end_exon_index: int, strand: str, direction: str, contig_length: int, - ) -> None: - start_exon = exons[start_exon_index] - end_exon = exons[end_exon_index] - if ( - not start_exon.start - or not start_exon.end - or not end_exon.start - or not end_exon.end - ): - raise ValueError("No start and end coordinates given.") - - # Adapt attributes for start/stop of coverage in sequence to the whole length of - # the contig even if not the whole contig was matched to the genome - if direction == ".": - contig_id, start_exon_contig_end, contig_start, _ = start_exon.attributes[ - "Target" - ][0].split(" ") - _, contig_end, end_exon_contig_start, _ = end_exon.attributes["Target"][ - 0 - ].split(" ") - contig_start = int(contig_start) - contig_end = int(contig_end) - - if start_exon == end_exon: - start_exon.attributes["Target"] = " ".join( - [contig_id, str(contig_length), "1", direction] - ) - else: - start_exon.attributes["Target"] = " ".join( - [contig_id, start_exon_contig_end, "1", direction] - ) - end_exon.attributes["Target"] = " ".join( - [contig_id, str(contig_length), end_exon_contig_start, direction] - ) - else: - contig_id, contig_start, start_exon_contig_end, _ = start_exon.attributes[ - "Target" - ][0].split(" ") - _, end_exon_contig_start, contig_end, _ = end_exon.attributes["Target"][ - 0 - ].split(" ") - contig_start = int(contig_start) - contig_end = int(contig_end) - - if start_exon == end_exon: - start_exon.attributes["Target"] = " ".join( - [contig_id, "1", str(contig_length), direction] - ) - else: - start_exon.attributes["Target"] = " ".join( - [contig_id, "1", start_exon_contig_end, direction] - ) - end_exon.attributes["Target"] = " ".join( - [contig_id, end_exon_contig_start, str(contig_length), direction] - ) - - exon_ids = [exon.attributes["ID"][0] for exon in exons] - cds_ids = [cds_entry.attributes["ID"][0] for cds_entry in cds] - mrna_id = mrna.attributes["ID"][0] - + contig: str, + ) -> str: # If direction is antisense, change alignment to be on the complementary strand if direction == "-": if strand == "+": @@ -192,115 +325,70 @@ def _write_new_feature_coordinates( mrna.strand = new_strand for exon in exons: exon.strand = new_strand - for cds_entry in cds: - cds_entry.strand = new_strand exons.reverse() - cds.reverse() - - start_exon_index = (len(exons) - 1) - start_exon_index - end_exon_index = (len(exons) - 1) - end_exon_index - - # Adapt the alignment genome coordinates to fit the whole contig length - if ( - (direction == "+" and strand == "+") - or (direction == "-" and strand == "-") - or (direction == "." and strand == "+") - ): - new_start = start_exon.start - (contig_start - 1) - start_exon.start = new_start - mrna.start = new_start - gene.start = new_start - new_end = end_exon.end + (contig_length - contig_end) - end_exon.end = new_end - mrna.end = new_end - gene.end = new_end + # TODO: Is assumption that exons are already in correct order + # if dir=indeterminate true? + # Potentially only labeled "indeterminate" if only one exon? + + exon_starts: List[int] = [] + exon_ends: List[int] = [] + for exon in exons: + if direction == ".": + _, exon_end, exon_start, _ = exon.attributes["Target"][0].split(" ") + else: + _, exon_start, exon_end, _ = exon.attributes["Target"][0].split(" ") + exon_starts.append(int(exon_start)) + exon_ends.append(int(exon_end)) + + cut_contig_parts: List[str] = [] + for i in range(len(exon_starts)): + current_start = exon_starts[i] + current_end = exon_ends[i] + current_part = contig[current_start - 1 : current_end] # noqa: E203 + cut_contig_parts.append(current_part) + start_end_cut_contig = "".join(cut_contig_parts) + + if len(start_end_cut_contig) < (0.7 * len(contig)): + # TODO: Filter out alignment and report + pass + + # TODO: Not needed? + # exon_ids = [exon.attributes["ID"][0] for exon in exons] + # cds_ids = [exon_id.replace("exon", "cds") for exon_id in exon_ids] + mrna_id = mrna.attributes["ID"][0] + + cls._write_feature_in_gtf_format( + output_gtf, + gene, + gene.attributes["ID"][0], + ) + for frame in range(3): + mrna.attributes["ID"] = mrna_id + "." + str(frame) cls._write_feature_in_gtf_format( - output_gtf, - gene, - gene.attributes["ID"][0], + output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] ) - for frame in range(3): - mrna.attributes["ID"] = mrna_id + "." + str(frame) + for exon in exons: + # TODO: Not needed? + # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) + # exon.attributes["Parent"] = mrna.attributes["ID"] cls._write_feature_in_gtf_format( - output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] + output_gtf, + exon, + gene.attributes["ID"][0], + mrna.attributes["ID"][0], ) - for exon_index, exon in enumerate(exons): - exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) - exon.attributes["Parent"] = mrna.attributes["ID"] - cls._write_feature_in_gtf_format( - output_gtf, - exon, - gene.attributes["ID"][0], - mrna.attributes["ID"][0], - ) - for cds_index, cds_entry in enumerate(cds): - cds_entry.start = exons[cds_index].start - cds_entry.end = exons[cds_index].end - if cds_index == start_exon_index: - cds_entry.start = cds_entry.start + frame # pyright: ignore - if cds_index == end_exon_index: - cds_entry.end = cds_entry.end - ( # pyright: ignore - (contig_length - frame) % 3 - ) - cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) - cds_entry.attributes["Parent"] = mrna.attributes["ID"] - cls._write_feature_in_gtf_format( - output_gtf, - cds_entry, - gene.attributes["ID"][0], - mrna.attributes["ID"][0], - ) - - elif ( - (direction == "+" and strand == "-") - or (direction == "-" and strand == "+") - or (direction == "." and strand == "-") - ): - new_end = start_exon.end + (contig_start - 1) - start_exon.end = new_end - mrna.end = new_end - gene.end = new_end - new_start = end_exon.start - (contig_length - contig_end) - end_exon.start = new_start - mrna.start = new_start - gene.start = new_start - - cls._write_feature_in_gtf_format(output_gtf, gene, gene.attributes["ID"][0]) - for frame in range(3): - mrna.attributes["ID"] = mrna_id + "." + str(frame) + for exon in exons: + exon.featuretype = "CDS" cls._write_feature_in_gtf_format( - output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] + output_gtf, + exon, + gene.attributes["ID"][0], + mrna.attributes["ID"][0], ) - for exon_index, exon in enumerate(exons): - exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) - exon.attributes["Parent"] = mrna.attributes["ID"] - cls._write_feature_in_gtf_format( - output_gtf, - exon, - gene.attributes["ID"][0], - mrna.attributes["ID"][0], - ) - for cds_index, cds_entry in enumerate(cds): - cds_entry.start = exons[cds_index].start - cds_entry.end = exons[cds_index].end - if cds_index == start_exon_index: - cds_entry.end = cds_entry.end - frame # pyright: ignore - if cds_index == end_exon_index: - cds_entry.start = cds_entry.start + ( # pyright: ignore - (contig_length - frame) % 3 - ) - cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) - cds_entry.attributes["Parent"] = mrna.attributes["ID"] - cls._write_feature_in_gtf_format( - output_gtf, - cds_entry, - gene.attributes["ID"][0], - mrna.attributes["ID"][0], - ) - else: - raise ValueError("Strand must be one of '+', '-'.") + + return start_end_cut_contig @classmethod def generate_gtf_input_file( @@ -308,11 +396,17 @@ def generate_gtf_input_file( path_to_gff: Path, output_directory: Path, sequence_lengths_per_contig: List[int], - ) -> List[int]: + contig_sequences: List[Tuple[str, str]], + ) -> Tuple[List[int], List[List[str]]]: # Track number of transcripts to write protein FASTA with matching ids number_of_transcripts_per_contig: List[int] = [ 0 for _ in range(len(sequence_lengths_per_contig)) ] + # Per original contig, there can be several new contigs + # based on different cutoffs + new_contig_sequences: List[List[str]] = [ + [] for _ in range(len(contig_sequences)) + ] gffutils_db = gffutils.create_db( path_to_gff.absolute().as_posix(), @@ -334,41 +428,10 @@ def generate_gtf_input_file( strand = first_exon.strand target: str = first_exon.attributes["Target"][0] contig_id, _, _, direction = target.split(" ") - contig_length = sequence_lengths_per_contig[int(contig_id[-1])] - if direction == ".": # indeterminate - start_exon_index = exons.index( - min( - exons, - key=lambda exon: int( - exon.attributes["Target"][0].split(" ")[2] - ), - ) - ) - end_exon_index = exons.index( - max( - exons, - key=lambda exon: int( - exon.attributes["Target"][0].split(" ")[1] - ), - ) - ) - else: # sense or antisense - start_exon_index = exons.index( - min( - exons, - key=lambda exon: int( - exon.attributes["Target"][0].split(" ")[1] - ), - ) - ) - end_exon_index = exons.index( - max( - exons, - key=lambda exon: int( - exon.attributes["Target"][0].split(" ")[2] - ), - ) - ) + contig_id = int(contig_id.split("-")[-1]) + contig_length = sequence_lengths_per_contig[contig_id] + contig = contig_sequences[contig_id] + mrna = [ gene_child for gene_child in gene_children @@ -377,45 +440,48 @@ def generate_gtf_input_file( 0 ] # There can be only one mRNA per gene mrna_id = mrna.attributes["ID"][0] - number_of_transcripts_per_contig[int(mrna_id.split(".")[0][-1])] += 1 + # TODO: Unify with the one above? + contig_idx = int(mrna_id.split(".")[0].split("-")[-1]) + number_of_transcripts_per_contig[contig_idx] += 1 - cls._write_new_feature_coordinates( + new_contig = cls._write_new_feature_coordinates( output_gtf, gene_feature, mrna, exons, - [ - gene_child - for gene_child in gene_children - if gene_child.featuretype == "CDS" - ], - start_exon_index, - end_exon_index, strand, direction, contig_length, + contig[1], ) + new_contig_sequences[contig_idx].append(new_contig) - return number_of_transcripts_per_contig + return (number_of_transcripts_per_contig, new_contig_sequences) @staticmethod def generate_protein_fasta_input_file( - contig_sequences: List[Tuple[str, str]], + contig_ids: List[str], + contig_sequences: List[List[str]], output_directory: Path, number_of_transcripts_per_contig: List[int], ) -> None: + # TODO: Adapt to new separation of ids and seqs with open( output_directory / "pogo_fasta_in.fa", "wt", encoding="utf-8" ) as output_file: - for contig_index, (contig_id, contig_sequence) in enumerate( - contig_sequences - ): - for translation, frame in get_three_frame_translations( - contig_sequence, False + for contig_id, contig_cut_sequences in zip(contig_ids, contig_sequences): + # TODO: This relies on the assumption that all paths for one contig are + # reported in the GMAP alignment in ascending numerical order. + # Can we really be sure about this? + for transcript_index, contig_sequence in enumerate( + contig_cut_sequences ): - for transcript_index in range( - number_of_transcripts_per_contig[contig_index] + for translation, frame in get_three_frame_translations( + contig_sequence, False ): + # for transcript_index in range( + # number_of_transcripts_per_contig[contig_index] + # ): gene_id = f"{contig_id}_path{str(transcript_index + 1)}" transcript_id = ( f"{contig_id}_mrna{str(transcript_index + 1)}_{str(frame)}" @@ -449,13 +515,17 @@ def generate_gtf_and_protein_files_for_directory( contig_sequences = cls._get_contig_sequences( path_to_directory / "resulting_contigs.fa" ) - number_of_transcripts_per_contig = cls.generate_gtf_input_file( - path_to_directory / "alignment_result.gff3", - path_to_directory, - [len(contig_sequence[1]) for contig_sequence in contig_sequences], + number_of_transcripts_per_contig, updated_contig_sequences = ( + cls.generate_gtf_input_file( + path_to_directory / "alignment_result.gff3", + path_to_directory, + [len(contig_sequence[1]) for contig_sequence in contig_sequences], + contig_sequences, + ) ) cls.generate_protein_fasta_input_file( - contig_sequences, + [contig_sequence[0] for contig_sequence in contig_sequences], + updated_contig_sequences, path_to_directory, number_of_transcripts_per_contig, ) From 5364bdc5a3c1092f73b801b56b15e24e9b1d70a8 Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Thu, 12 Sep 2024 23:15:19 +0200 Subject: [PATCH 07/10] Rewrite CDS features based on exons --- .../output_generation/pogo_input_helper.py | 44 ++++++++++++++++--- 1 file changed, 37 insertions(+), 7 deletions(-) diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index 9310e94..788c81f 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -317,14 +317,14 @@ def _write_new_feature_coordinates( # If direction is antisense, change alignment to be on the complementary strand if direction == "-": if strand == "+": - new_strand = "-" + strand = "-" else: - new_strand = "+" + strand = "+" - gene.strand = new_strand - mrna.strand = new_strand + gene.strand = strand + mrna.strand = strand for exon in exons: - exon.strand = new_strand + exon.strand = strand exons.reverse() @@ -358,6 +358,8 @@ def _write_new_feature_coordinates( # exon_ids = [exon.attributes["ID"][0] for exon in exons] # cds_ids = [exon_id.replace("exon", "cds") for exon_id in exon_ids] mrna_id = mrna.attributes["ID"][0] + exon_start_coords = [exon.start for exon in exons] + exon_end_coords = [exon.end for exon in exons] cls._write_feature_in_gtf_format( output_gtf, @@ -369,18 +371,46 @@ def _write_new_feature_coordinates( cls._write_feature_in_gtf_format( output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] ) - for exon in exons: + for exon_idx, exon in enumerate(exons): # TODO: Not needed? # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) # exon.attributes["Parent"] = mrna.attributes["ID"] + exon.featuretype = "exon" + exon.start = exon_start_coords[exon_idx] + exon.end = exon_end_coords[exon_idx] cls._write_feature_in_gtf_format( output_gtf, exon, gene.attributes["ID"][0], mrna.attributes["ID"][0], ) - for exon in exons: + for exon_idx, exon in enumerate(exons): + # TODO: Is there a better solution, + # e.g. copying and modifying the feature? exon.featuretype = "CDS" + # TODO + # problem: overwrite of exon coords over the three iterations + # -> save them + # strand = +, dir = sense -> add frame to start of first CDS + # strand = +, dir = antisense -> subtract frame from end of first CDS (is first after reversing) + # strand = -, dir = sense -> subtract frame from end of first CDS + # strand = -, dir = antisense -> add frame to start of first CDS (is first after reversing) + # --> differentiation between +/- strand should suffice after reversing + if strand == "+": + if exon_idx == 0: + exon.start = exon.start + frame # pyright: ignore + if exon_idx == (len(exons) - 1): + exon.end = exon.end - ( # pyright: ignore + (len(start_end_cut_contig) - frame) % 3 + ) + else: + if exon_idx == 0: + exon.end = exon.end - frame # pyright: ignore + if exon_idx == (len(exons) - 1): + exon.start = exon.start + ( # pyright: ignore + (len(start_end_cut_contig) - frame) % 3 + ) + cls._write_feature_in_gtf_format( output_gtf, exon, From 2f4e2be3d27b6f6a2912df00f20ef97d790a7bde Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Fri, 20 Sep 2024 08:42:20 +0200 Subject: [PATCH 08/10] Sort exons before rewrite to gtf --- .../output_generation/pogo_input_helper.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index 788c81f..afb6e87 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -75,6 +75,10 @@ def generate_all_peptide_input_files( output_directory, merged_indexes[set_index] ) + @staticmethod + def _get_exon_feature_id(exon: gffutils.Feature) -> int: + return int(exon.attributes["ID"][0].split(".")[-1].replace("exon", "")) + @staticmethod def _write_feature_in_gtf_format( output_gtf: TextIO, @@ -302,6 +306,7 @@ def _write_feature_in_gtf_format( # else: # raise ValueError("Strand must be one of '+', '-'.") + # TODO: Remove unneeded arguments @classmethod def _write_new_feature_coordinates( cls, @@ -314,6 +319,15 @@ def _write_new_feature_coordinates( contig_length: int, contig: str, ) -> str: + # The exons need to be sorted to follow the same order as in the original GFF. + # This is necessary because gffutils does not necessarily return the children + # of a feature in order when calling children(). The parameter order_by + # cannot be used here because it does not allow us to select for an + # exon identifier or specify the same order as in the original. + # TODO: Is there an easier option for ordering + # that does not involve string splitting? + exons.sort(key=cls._get_exon_feature_id) + # If direction is antisense, change alignment to be on the complementary strand if direction == "-": if strand == "+": @@ -388,9 +402,6 @@ def _write_new_feature_coordinates( # TODO: Is there a better solution, # e.g. copying and modifying the feature? exon.featuretype = "CDS" - # TODO - # problem: overwrite of exon coords over the three iterations - # -> save them # strand = +, dir = sense -> add frame to start of first CDS # strand = +, dir = antisense -> subtract frame from end of first CDS (is first after reversing) # strand = -, dir = sense -> subtract frame from end of first CDS @@ -488,6 +499,7 @@ def generate_gtf_input_file( return (number_of_transcripts_per_contig, new_contig_sequences) + # TODO: Remove unneeded arguments @staticmethod def generate_protein_fasta_input_file( contig_ids: List[str], From 223bfc23dbedd6edb8ad1dccb6c30241d8daa3f3 Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Sat, 5 Apr 2025 19:21:15 +0200 Subject: [PATCH 09/10] WIP: Handle indels in alignment results --- pepti_map/main.py | 21 +- .../output_generation/pogo_input_helper.py | 231 ++---------------- 2 files changed, 44 insertions(+), 208 deletions(-) diff --git a/pepti_map/main.py b/pepti_map/main.py index 8418cd4..240f3f1 100644 --- a/pepti_map/main.py +++ b/pepti_map/main.py @@ -248,7 +248,11 @@ def align_reads_to_genome( return new_results_paths -def generate_pogo_input(paths_to_subdirectories: List[Path], peptide_file: str) -> None: +def generate_pogo_input( + paths_to_subdirectories: List[Path], + peptide_file: str, + no_indels: bool +) -> None: pogo_input_helper = PoGoInputHelper( Path(peptide_file), PATH_PEPTIDE_TO_CLUSTER_MAPPING_FILE ) @@ -256,7 +260,7 @@ def generate_pogo_input(paths_to_subdirectories: List[Path], peptide_file: str) paths_to_subdirectories, PATH_TO_MERGED_INDEXES ) PoGoInputHelper.generate_gtf_and_protein_files_for_multiple_directories( - paths_to_subdirectories + paths_to_subdirectories, no_indels ) _write_last_step(Step.POGO_INPUT.value) logging.info("Generated PoGo input files.") @@ -357,6 +361,7 @@ def concat_output(paths_to_subdirectories: List[Path], output_dir: str) -> None: "-pi", "--precompute-intersections", is_flag=True, + default=False, help=( "If used, the intersection sizes for the Jaccard Index " "calculation are precomputed during the matching phase." @@ -431,6 +436,15 @@ def concat_output(paths_to_subdirectories: List[Path], output_dir: str) -> None: show_default=True, help="Sets the '--min-identity' option for GMAP during alignment.", ) +@click.option( + "-ni", + "--no-indels", + required=False, + is_flag=True, + default=False, + help=("If set, contig alignments containing indels " + "are excluded from further processing.") +) def main( peptide_file: str, rna_file: str, @@ -446,6 +460,7 @@ def main( gmap_index: Union[str, None], min_trimmed_coverage: float, min_identity: float, + no_indels: bool ): _setup(output_dir) @@ -525,7 +540,7 @@ def main( if last_step < Step.POGO_INPUT.value: logging.info("Generating input files for PoGo.") - generate_pogo_input(paths_to_subdirectories, peptide_file) + generate_pogo_input(paths_to_subdirectories, peptide_file, no_indels) else: logging.info("Using already generated PoGo input files.") diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index afb6e87..8919666 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -1,4 +1,5 @@ from collections import defaultdict +from functools import partial import logging import multiprocessing from dotenv import dotenv_values @@ -112,200 +113,6 @@ def _write_feature_in_gtf_format( + "\n" ) - # @classmethod - # def _write_new_feature_coordinates_old( - # cls, - # output_gtf: TextIO, - # gene: gffutils.Feature, - # mrna: gffutils.Feature, - # exons: List[gffutils.Feature], - # cds: List[gffutils.Feature], - # start_exon_index: int, - # end_exon_index: int, - # strand: str, - # direction: str, - # contig_length: int, - # ) -> None: - # start_exon = exons[start_exon_index] - # end_exon = exons[end_exon_index] - # if ( - # not start_exon.start - # or not start_exon.end - # or not end_exon.start - # or not end_exon.end - # ): - # raise ValueError("No start and end coordinates given.") - - # # Adapt attributes for start/stop of coverage in sequence to the whole length of - # # the contig even if not the whole contig was matched to the genome - # if direction == ".": - # contig_id, start_exon_contig_end, contig_start, _ = start_exon.attributes[ - # "Target" - # ][0].split(" ") - # _, contig_end, end_exon_contig_start, _ = end_exon.attributes["Target"][ - # 0 - # ].split(" ") - # contig_start = int(contig_start) - # contig_end = int(contig_end) - - # if start_exon == end_exon: - # start_exon.attributes["Target"] = " ".join( - # [contig_id, str(contig_length), "1", direction] - # ) - # else: - # start_exon.attributes["Target"] = " ".join( - # [contig_id, start_exon_contig_end, "1", direction] - # ) - # end_exon.attributes["Target"] = " ".join( - # [contig_id, str(contig_length), end_exon_contig_start, direction] - # ) - # else: - # contig_id, contig_start, start_exon_contig_end, _ = start_exon.attributes[ - # "Target" - # ][0].split(" ") - # _, end_exon_contig_start, contig_end, _ = end_exon.attributes["Target"][ - # 0 - # ].split(" ") - # contig_start = int(contig_start) - # contig_end = int(contig_end) - - # if start_exon == end_exon: - # start_exon.attributes["Target"] = " ".join( - # [contig_id, "1", str(contig_length), direction] - # ) - # else: - # start_exon.attributes["Target"] = " ".join( - # [contig_id, "1", start_exon_contig_end, direction] - # ) - # end_exon.attributes["Target"] = " ".join( - # [contig_id, end_exon_contig_start, str(contig_length), direction] - # ) - - # exon_ids = [exon.attributes["ID"][0] for exon in exons] - # cds_ids = [cds_entry.attributes["ID"][0] for cds_entry in cds] - # mrna_id = mrna.attributes["ID"][0] - - # # If direction is antisense, change alignment to be on the complementary strand - # if direction == "-": - # if strand == "+": - # new_strand = "-" - # else: - # new_strand = "+" - - # gene.strand = new_strand - # mrna.strand = new_strand - # for exon in exons: - # exon.strand = new_strand - # for cds_entry in cds: - # cds_entry.strand = new_strand - - # exons.reverse() - # cds.reverse() - - # start_exon_index = (len(exons) - 1) - start_exon_index - # end_exon_index = (len(exons) - 1) - end_exon_index - - # # Adapt the alignment genome coordinates to fit the whole contig length - # if ( - # (direction == "+" and strand == "+") - # or (direction == "-" and strand == "-") - # or (direction == "." and strand == "+") - # ): - # new_start = start_exon.start - (contig_start - 1) - # start_exon.start = new_start - # mrna.start = new_start - # gene.start = new_start - # new_end = end_exon.end + (contig_length - contig_end) - # end_exon.end = new_end - # mrna.end = new_end - # gene.end = new_end - - # cls._write_feature_in_gtf_format( - # output_gtf, - # gene, - # gene.attributes["ID"][0], - # ) - # for frame in range(3): - # mrna.attributes["ID"] = mrna_id + "." + str(frame) - # cls._write_feature_in_gtf_format( - # output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] - # ) - # for exon_index, exon in enumerate(exons): - # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) - # exon.attributes["Parent"] = mrna.attributes["ID"] - # cls._write_feature_in_gtf_format( - # output_gtf, - # exon, - # gene.attributes["ID"][0], - # mrna.attributes["ID"][0], - # ) - # for cds_index, cds_entry in enumerate(cds): - # cds_entry.start = exons[cds_index].start - # cds_entry.end = exons[cds_index].end - # if cds_index == start_exon_index: - # cds_entry.start = cds_entry.start + frame # pyright: ignore - # if cds_index == end_exon_index: - # cds_entry.end = cds_entry.end - ( # pyright: ignore - # (contig_length - frame) % 3 - # ) - # cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) - # cds_entry.attributes["Parent"] = mrna.attributes["ID"] - # cls._write_feature_in_gtf_format( - # output_gtf, - # cds_entry, - # gene.attributes["ID"][0], - # mrna.attributes["ID"][0], - # ) - - # elif ( - # (direction == "+" and strand == "-") - # or (direction == "-" and strand == "+") - # or (direction == "." and strand == "-") - # ): - # new_end = start_exon.end + (contig_start - 1) - # start_exon.end = new_end - # mrna.end = new_end - # gene.end = new_end - # new_start = end_exon.start - (contig_length - contig_end) - # end_exon.start = new_start - # mrna.start = new_start - # gene.start = new_start - - # cls._write_feature_in_gtf_format(output_gtf, gene, gene.attributes["ID"][0]) - # for frame in range(3): - # mrna.attributes["ID"] = mrna_id + "." + str(frame) - # cls._write_feature_in_gtf_format( - # output_gtf, mrna, gene.attributes["ID"][0], mrna.attributes["ID"][0] - # ) - # for exon_index, exon in enumerate(exons): - # exon.attributes["ID"] = exon_ids[exon_index] + "." + str(frame) - # exon.attributes["Parent"] = mrna.attributes["ID"] - # cls._write_feature_in_gtf_format( - # output_gtf, - # exon, - # gene.attributes["ID"][0], - # mrna.attributes["ID"][0], - # ) - # for cds_index, cds_entry in enumerate(cds): - # cds_entry.start = exons[cds_index].start - # cds_entry.end = exons[cds_index].end - # if cds_index == start_exon_index: - # cds_entry.end = cds_entry.end - frame # pyright: ignore - # if cds_index == end_exon_index: - # cds_entry.start = cds_entry.start + ( # pyright: ignore - # (contig_length - frame) % 3 - # ) - # cds_entry.attributes["ID"] = cds_ids[cds_index] + "." + str(frame) - # cds_entry.attributes["Parent"] = mrna.attributes["ID"] - # cls._write_feature_in_gtf_format( - # output_gtf, - # cds_entry, - # gene.attributes["ID"][0], - # mrna.attributes["ID"][0], - # ) - # else: - # raise ValueError("Strand must be one of '+', '-'.") - # TODO: Remove unneeded arguments @classmethod def _write_new_feature_coordinates( @@ -438,6 +245,7 @@ def generate_gtf_input_file( output_directory: Path, sequence_lengths_per_contig: List[int], contig_sequences: List[Tuple[str, str]], + no_indels=False, ) -> Tuple[List[int], List[List[str]]]: # Track number of transcripts to write protein FASTA with matching ids number_of_transcripts_per_contig: List[int] = [ @@ -459,6 +267,22 @@ def generate_gtf_input_file( ) as output_gtf: for gene_feature in gffutils_db.features_of_type("gene"): gene_children = list(gffutils_db.children(gene_feature)) + + mrna = [ + gene_child + for gene_child in gene_children + if gene_child.featuretype == "mRNA" + ][ + 0 + ] # There can be only one mRNA per gene + + # TODO: If no indels allowed: Check if feature contains indels -> If so, exclude (make sure that returned number of transcripts matches up for creation of fasta file) + if no_indels: + mrna_indels = int(mrna.attributes["indels"][0]) + if mrna_indels != 0: + # TODO + pass + exons = [ gene_child for gene_child in gene_children @@ -473,13 +297,6 @@ def generate_gtf_input_file( contig_length = sequence_lengths_per_contig[contig_id] contig = contig_sequences[contig_id] - mrna = [ - gene_child - for gene_child in gene_children - if gene_child.featuretype == "mRNA" - ][ - 0 - ] # There can be only one mRNA per gene mrna_id = mrna.attributes["ID"][0] # TODO: Unify with the one above? contig_idx = int(mrna_id.split(".")[0].split("-")[-1]) @@ -552,7 +369,7 @@ def _get_contig_sequences(path_to_contig_sequences: Path) -> List[Tuple[str, str @classmethod def generate_gtf_and_protein_files_for_directory( - cls, path_to_directory: Path + cls, path_to_directory: Path, no_indels=False ) -> None: contig_sequences = cls._get_contig_sequences( path_to_directory / "resulting_contigs.fa" @@ -563,6 +380,7 @@ def generate_gtf_and_protein_files_for_directory( path_to_directory, [len(contig_sequence[1]) for contig_sequence in contig_sequences], contig_sequences, + no_indels, ) ) cls.generate_protein_fasta_input_file( @@ -574,8 +392,7 @@ def generate_gtf_and_protein_files_for_directory( @classmethod def generate_gtf_and_protein_files_for_multiple_directories( - cls, - paths_to_directories: List[Path], + cls, paths_to_directories: List[Path], no_indels=False ) -> None: # TODO: Refactor code duplication try: @@ -589,5 +406,9 @@ def generate_gtf_and_protein_files_for_multiple_directories( ) with multiprocessing.Pool(n_processes) as pool: pool.map( - cls.generate_gtf_and_protein_files_for_directory, paths_to_directories + partial( + cls.generate_gtf_and_protein_files_for_directory, + no_indels=no_indels, + ), + paths_to_directories, ) From 63a31eda488828e662ea93f5e1f8a3a9f88cb98c Mon Sep 17 00:00:00 2001 From: Clara Uktar Date: Sun, 6 Apr 2025 13:46:19 +0200 Subject: [PATCH 10/10] Exclude alignments with indels if parameter active --- .../output_generation/pogo_input_helper.py | 68 +++++++++---------- 1 file changed, 31 insertions(+), 37 deletions(-) diff --git a/pepti_map/output_generation/pogo_input_helper.py b/pepti_map/output_generation/pogo_input_helper.py index 8919666..58bf853 100644 --- a/pepti_map/output_generation/pogo_input_helper.py +++ b/pepti_map/output_generation/pogo_input_helper.py @@ -113,7 +113,6 @@ def _write_feature_in_gtf_format( + "\n" ) - # TODO: Remove unneeded arguments @classmethod def _write_new_feature_coordinates( cls, @@ -123,7 +122,6 @@ def _write_new_feature_coordinates( exons: List[gffutils.Feature], strand: str, direction: str, - contig_length: int, contig: str, ) -> str: # The exons need to be sorted to follow the same order as in the original GFF. @@ -209,10 +207,14 @@ def _write_new_feature_coordinates( # TODO: Is there a better solution, # e.g. copying and modifying the feature? exon.featuretype = "CDS" - # strand = +, dir = sense -> add frame to start of first CDS - # strand = +, dir = antisense -> subtract frame from end of first CDS (is first after reversing) - # strand = -, dir = sense -> subtract frame from end of first CDS - # strand = -, dir = antisense -> add frame to start of first CDS (is first after reversing) + # strand = +, dir = sense + # -> add frame to start of first CDS + # strand = +, dir = antisense + # -> subtract frame from end of first CDS (is first after reversing) + # strand = -, dir = sense + # -> subtract frame from end of first CDS + # strand = -, dir = antisense + # -> add frame to start of first CDS (is first after reversing) # --> differentiation between +/- strand should suffice after reversing if strand == "+": if exon_idx == 0: @@ -243,13 +245,12 @@ def generate_gtf_input_file( cls, path_to_gff: Path, output_directory: Path, - sequence_lengths_per_contig: List[int], contig_sequences: List[Tuple[str, str]], no_indels=False, - ) -> Tuple[List[int], List[List[str]]]: - # Track number of transcripts to write protein FASTA with matching ids - number_of_transcripts_per_contig: List[int] = [ - 0 for _ in range(len(sequence_lengths_per_contig)) + ) -> Tuple[List[List[int]], List[List[str]]]: + # Track contig alignment ids to write protein FASTA with matching ids + alignment_ids_per_contig: List[List[int]] = [ + [] for _ in range(len(contig_sequences)) ] # Per original contig, there can be several new contigs # based on different cutoffs @@ -276,12 +277,12 @@ def generate_gtf_input_file( 0 ] # There can be only one mRNA per gene - # TODO: If no indels allowed: Check if feature contains indels -> If so, exclude (make sure that returned number of transcripts matches up for creation of fasta file) + # If no indels allowed: Check if feature contains indels + # -> If so, gene feature is skipped if no_indels: mrna_indels = int(mrna.attributes["indels"][0]) if mrna_indels != 0: - # TODO - pass + continue exons = [ gene_child @@ -294,14 +295,8 @@ def generate_gtf_input_file( target: str = first_exon.attributes["Target"][0] contig_id, _, _, direction = target.split(" ") contig_id = int(contig_id.split("-")[-1]) - contig_length = sequence_lengths_per_contig[contig_id] contig = contig_sequences[contig_id] - mrna_id = mrna.attributes["ID"][0] - # TODO: Unify with the one above? - contig_idx = int(mrna_id.split(".")[0].split("-")[-1]) - number_of_transcripts_per_contig[contig_idx] += 1 - new_contig = cls._write_new_feature_coordinates( output_gtf, gene_feature, @@ -309,12 +304,16 @@ def generate_gtf_input_file( exons, strand, direction, - contig_length, contig[1], ) + + mrna_id = mrna.attributes["ID"][0] + contig_idx = int(mrna_id.split(".")[0].split("-")[-1]) + path_number = int(mrna_id.split(".")[1].replace("mrna", "")) new_contig_sequences[contig_idx].append(new_contig) + alignment_ids_per_contig[contig_idx].append(path_number) - return (number_of_transcripts_per_contig, new_contig_sequences) + return (alignment_ids_per_contig, new_contig_sequences) # TODO: Remove unneeded arguments @staticmethod @@ -322,28 +321,24 @@ def generate_protein_fasta_input_file( contig_ids: List[str], contig_sequences: List[List[str]], output_directory: Path, - number_of_transcripts_per_contig: List[int], + alignment_ids_per_contig: List[List[int]], ) -> None: # TODO: Adapt to new separation of ids and seqs with open( output_directory / "pogo_fasta_in.fa", "wt", encoding="utf-8" ) as output_file: - for contig_id, contig_cut_sequences in zip(contig_ids, contig_sequences): - # TODO: This relies on the assumption that all paths for one contig are - # reported in the GMAP alignment in ascending numerical order. - # Can we really be sure about this? - for transcript_index, contig_sequence in enumerate( - contig_cut_sequences + for contig_id, contig_cut_sequences, alignment_ids in zip( + contig_ids, contig_sequences, alignment_ids_per_contig + ): + for alignment_id, contig_sequence in zip( + alignment_ids, contig_cut_sequences ): for translation, frame in get_three_frame_translations( contig_sequence, False ): - # for transcript_index in range( - # number_of_transcripts_per_contig[contig_index] - # ): - gene_id = f"{contig_id}_path{str(transcript_index + 1)}" + gene_id = f"{contig_id}_path{str(alignment_id)}" transcript_id = ( - f"{contig_id}_mrna{str(transcript_index + 1)}_{str(frame)}" + f"{contig_id}_mrna{str(alignment_id)}_{str(frame)}" ) output_file.write( ( @@ -374,11 +369,10 @@ def generate_gtf_and_protein_files_for_directory( contig_sequences = cls._get_contig_sequences( path_to_directory / "resulting_contigs.fa" ) - number_of_transcripts_per_contig, updated_contig_sequences = ( + alignment_ids_per_contig, updated_contig_sequences = ( cls.generate_gtf_input_file( path_to_directory / "alignment_result.gff3", path_to_directory, - [len(contig_sequence[1]) for contig_sequence in contig_sequences], contig_sequences, no_indels, ) @@ -387,7 +381,7 @@ def generate_gtf_and_protein_files_for_directory( [contig_sequence[0] for contig_sequence in contig_sequences], updated_contig_sequences, path_to_directory, - number_of_transcripts_per_contig, + alignment_ids_per_contig, ) @classmethod