Merge pull request #22 from Juke34/reditools3

add reditools3

Author: eascarrunz

bin/stats/pluviometer.py

@@ -13,7 +13,7 @@
 import numpy as np
 from numpy.typing import NDArray
 from typing import TextIO, NamedTuple, Optional, Generator
-from site_variant_readers import RNAVariantReader, ReditoolsReader
+from site_variant_readers import RNAVariantReader, Reditools2Reader, Reditools3Reader
 import argparse
 from contextlib import nullcontext
 import sys
@@ -52,8 +52,8 @@ def parse_cli_input() -> argparse.Namespace:
         "--format",
         "-f",
         type=str,
-        choices=["reditools", "jacusa2", "sapin"],
-        default="reditools",
+        choices=["reditools2", "reditools3", "jacusa2", "sapin"],
+        default="reditools3",
         help="Sites file format",
     )
     parser.add_argument(
@@ -437,7 +437,8 @@ def update_queues(self, pos: int) -> QueueUpdates:
             if item.feature.location.end < pos:
                 skipped += 1
                 item.manager.update_progress(self.output_handle)
-                self.progress_bar.next()
+                if self.use_progress_bar:
+                    self.progress_bar.next()
             else:
                 self.load_active_queue(item)
                 activated += 1
@@ -562,9 +563,12 @@ def scan_and_count(self, reader: RNAVariantReader, header=True) -> None:
         else nullcontext(sys.stdout) as output_handle,
     ):
         match args.format:
-            case "reditools":
-                print("Reditools format\n")
-                sv_reader: RNAVariantReader = ReditoolsReader(sv_handle)
+            case "reditools2":
+                print("Reditools2 format\n")
+                sv_reader: RNAVariantReader = Reditools2Reader(sv_handle)
+            case "reditools3":
+                print("Reditools3 format\n")
+                sv_reader: RNAVariantReader = Reditools3Reader(sv_handle)
             case _:
                 raise Exception(f'Unimplemented format "{args.format}"')
 

bin/stats/site_variant_readers.py

@@ -58,7 +58,7 @@ def close(self) -> None:
         pass
 
 
-class ReditoolsReader(RNAVariantReader):
+class ReditoolsXReader(RNAVariantReader):
     header_strings = (
         "Region",
         "Position",
@@ -112,18 +112,12 @@ def _validate_header(self) -> None:
                 )
 
         return None
+    
+    def parse_strand(self):
+        pass
 
     def _parse_parts(self) -> SiteVariantData:
-        strand = int(self.parts[REDITOOLS_FIELD_INDEX["Strand"]])
-        match strand:
-            case 0:
-                strand = -1
-            case 1:
-                strand = 1
-            case 2:
-                strand = 0
-            case _:
-                raise Exception(f"Invalid strand value: {strand}")
+        strand = self.parse_strand()
 
         reference_nuc_str: str = self.parts[REDITOOLS_FIELD_INDEX["Reference"]]
 
@@ -153,3 +147,30 @@ def close(self) -> None:
         """Close the file"""
         self.file_handle.close()
 
+
+class Reditools2Reader(ReditoolsXReader):
+    def parse_strand(self) -> int:
+        strand = int(self.parts[REDITOOLS_FIELD_INDEX["Strand"]])
+        match strand:
+            case 0:
+                return -1
+            case 1:
+                return 1
+            case 2:
+                return 0
+            case _:
+                raise Exception(f"Invalid strand value: {strand}")
+
+class Reditools3Reader(ReditoolsXReader):
+    def parse_strand(self) -> int:
+        strand_str = self.parts[REDITOOLS_FIELD_INDEX["Strand"]]
+        match strand_str:
+            case "-":
+                return -1
+            case "+":
+                return 1
+            case "*":
+                return 0
+            case _:
+                raise Exception(f"Invalid strand value: {strand_str}")
+

bin/stats/tiny_pluviometer.py

@@ -2,7 +2,7 @@
 from BCBio import GFF
 from Bio import SeqIO, SeqRecord
 from Bio.SeqFeature import SeqFeature
-from site_variant_readers import ReditoolsReader
+from site_variant_readers import Reditools2Reader
 from typing import Optional
 from utils import SiteVariantData, NUC_STR_TO_IND
 import numpy as np
@@ -61,7 +61,7 @@ def find_features(feature_id: str, feature: SeqFeature, found: list[SeqFeature])
 
         for feature in target_features:
             with open(args.sites) as sites_handle:
-                reader = ReditoolsReader(sites_handle)
+                reader = Reditools2Reader(sites_handle)
                 variant_data: SiteVariantData = reader.read()
 
                 while variant_data.position <= feature.location.end:

config/softwares.config

@@ -36,6 +36,9 @@ process {
     withLabel: "reditools2" {
         container = singularity.enabled ? "${sifPath}/reditools2.sif" : "reditools2"
     }
+    withLabel: "reditools3" {
+        container = singularity.enabled ? "${sifPath}/reditools3.sif" : "reditools3"
+    }
     withLabel: 'rsem' {
         container = 'quay.io/biocontainers/rsem:1.3.3--pl5321h0033a41_7'
     }

docker/reditools3/Dockerfile

@@ -0,0 +1,10 @@
+FROM python:slim-bullseye
+
+RUN apt update
+
+# Install procps to provide ps, required by Nextflow
+RUN apt install -y procps
+
+RUN pip install REDItools3
+
+CMD [ "bash" ]

modules/reditools3.nf

@@ -0,0 +1,34 @@
+process reditools3 {
+    label "reditools3"
+    publishDir("${params.outdir}/reditools3", mode: "copy")
+    tag "${meta.id}"
+
+    input:
+        tuple(val(meta), path(bam), path(bamindex))
+        path genome
+
+    output:
+        tuple(val(meta), path("edit_table.txt"), emit: tuple_sample_serial_table)
+
+    script:
+        // Set the strand orientation parameter from the library type parameter
+        // Terms explained in https://salmon.readthedocs.io/en/latest/library_type.html
+        if (meta.libtype in ["ISR", "SR"]) {
+            // First-strand oriented
+            strand_orientation = "2"
+        } else if (meta.libtype in ["ISF", "SF"]) {
+            // Second-strand oriented
+            strand_orientation = "1"
+        } else if (meta.libtype in ["IU", "U"]) {
+            // Unstranded
+            strand_orientation = "0"
+        } else {
+            // Unsupported: Pass the library type string so that it's reported in
+            // the reditools error message
+            strand_orientation = meta.libtype
+        }
+
+        """
+        python -m reditools analyze ${bam} --reference ${genome} --strand ${strand_orientation} --output-file edit_table.txt
+        """
+}

rain.nf

@@ -26,8 +26,8 @@ params.library_type = "auto" // can be 'U', 'IU', 'MU', 'OU', 'ISF', 'ISR', 'MSF
 params.bam_library_type = null // can be 'U', 'IU', 'MU', 'OU', 'ISF', 'ISR', 'MSF', 'MSR', 'OSF', 'OSR', 'auto' - see https://github.com/Juke34/AliNe for more information
 
 // Edit counting params
-edit_site_tools = ["reditools2", "jacusa2", "sapin"]
-params.edit_site_tool = "reditools2"
+edit_site_tools = ["reditools2", "reditools3", "jacusa2", "sapin"]
+params.edit_site_tool = "reditools3"
 params.edit_threshold = 1
 params.aggregation_mode = "all"
 
@@ -140,6 +140,7 @@ include {multiqc} from './modules/multiqc.nf'
 include {fasta_uncompress} from "$baseDir/modules/pigz.nf"
 include {samtools_index; samtools_fasta_index; samtools_sort_bam} from './modules/samtools.nf'
 include {reditools2} from "./modules/reditools2.nf"
+include {reditools3} from "./modules/reditools3.nf"
 include {jacusa2} from "./modules/jacusa2.nf"
 include {sapin} from "./modules/sapin.nf"
 include {normalize_gxf} from "./modules/agat.nf"
@@ -465,7 +466,12 @@ workflow rain {
             case "reditools2":
                 reditools2(samtools_index.out.tuple_sample_bam_bamindex, genome, params.region)
                 normalize_gxf(annnotation)
-                pluviometer(reditools2.out.tuple_sample_serial_table, normalize_gxf.out.gff, "reditools")
+                pluviometer(reditools2.out.tuple_sample_serial_table, normalize_gxf.out.gff, "reditools2")
+                break
+            case "reditools3":
+                reditools3(samtools_index.out.tuple_sample_bam_bamindex, genome)
+                normalize_gxf(annnotation)
+                pluviometer(reditools3.out.tuple_sample_serial_table, normalize_gxf.out.gff, "reditools3")
                 break
             default:
                 exit(1, "Wrong edit site tool was passed")

singularity/reditools3/reditools3.def

@@ -0,0 +1,19 @@
+Bootstrap: docker
+From: python:slim-bullseye
+
+%post
+    apt-get update && apt-get install -y procps
+    pip install REDItools3
+
+%environment
+    export PATH=/usr/local/bin/:$PATH
+
+%runscript
+    exec bash "$@"
+
+%labels
+    Author eascarrunz.dev@mailfence.com
+    Tool REDItools3
+
+%help
+    Container for REDItools3.