fixed filtering SE reads (#52)

* fixed filtering SE reads

Author: riasc

CHANGELOG.md

@@ -8,10 +8,19 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [0.5.x] -2023-xx-xx (unreleased)
 
+
 ### Features
 
 - NMD information (e.g., escape rule,...) is now also calculated for all variants
 
+## [0.3.4] - 2025-04-04
+
+### Fix 
+
+- Moved script for preparing hla input to workflow/scripts/genotyping
+- Added missing pipe to log file in rule hlatyping_mhcII
+- Added wrong input file to filtering of mhcII reads on SE
+
 ## [0.3.3] - 2025-04-02 
 
 ### Fix 

workflow/rules/common.smk

@@ -255,6 +255,17 @@ def get_all_mhcI_alleles(wildcards):
 
 
 ##### MHC CLASS II #####
+def get_input_filter_reads_mhcII_SE(wildcards):
+  if config["preproc"]["activate"]:
+    return expand("results/{sample}/{seqtype}/reads/{group}_preproc.fq.gz",
+                  sample=wildcards.sample,
+                  seqtype = "dnaseq" if wildcards.nartype == "DNA" else "rnaseq",
+                  group=wildcards.group)
+  else:
+    seqtype = "dnaseq" if wildcards.nartype == "DNA" else "rnaseq"
+    return config["data"][f"{seqtype}"][wildcards.group]
+
+
 def get_input_filter_reads_mhcII_PE(wildcards):
   seqtype = "dnaseq" if wildcards.nartype == "DNA" else "rnaseq"
   """ if the reads are in BAM format, we consider them as processed and forward
@@ -269,6 +280,7 @@ def get_input_filter_reads_mhcII_PE(wildcards):
     if config["preproc"]["activate"]:
       return expand("results/{sample}/{seqtype}/reads/{group}_{readpair}_preproc.fq.gz",
                     sample=wildcards.sample,
+                    seqtype= "dnaseq" if wildcards.nartype == "DNA" else "rnaseq",
                     group = wildcards.group,
                     readpair=["R1", "R2"]
       )

workflow/rules/hlatyping_mhcII.smk

@@ -1,7 +1,7 @@
 ######### MHC-II HLA GENOTYPING ###########
 rule filter_reads_mhcII_SE:
   input:
-    sample=["results/{sample}/hla/reads/{group}_{nartype}_SE.fq"], 
+    sample=get_input_filter_reads_mhcII_SE,
     idx=multiext(
       "resources/hla/bowtie2_index",
       ".1.bt2",
@@ -72,12 +72,29 @@ rule finalize_reads_mhcII:
   threads: 1
   shell:
     """
-      python workflow/scripts/finalize_mhcII_input.py \
+      python workflow/scripts/genotyping/finalize_mhcII_input.py \
           {input} \
           {output.fwd} \
           {output.rev}
     """
 
+rule readcounts_mhcII:
+  input: 
+    fwd = "results/{sample}/hla/mhc-II/reads/{group}_{nartype}_final_R1.fq",
+    rev = "results/{sample}/hla/mhc-II/reads/{group}_{nartype}_final_R2.fq"
+  output:
+    "results/{sample}/hla/mhc-II/reads/{group}_{nartype}_final_counts.txt"
+  log:
+    "logs/{sample}/genotyping/final_readcounts_mhcII_{group}_{nartype}.log"
+  conda:
+    "../envs/basic.yml"
+  threads: 1
+  shell:
+    """
+      python workflow/scripts/countreads.py \
+          {input.fwd} {output} > {log} 2>&1 | exit 0
+    """
+
 
 rule hlatyping_mhcII:
   input:
@@ -104,7 +121,7 @@ rule hlatyping_mhcII:
           {input.fwd} {input.rev} \
           {params.split} {params.dic} \
           {wildcards.group}_{wildcards.nartype} \
-          results/{wildcards.sample}/hla/mhc-II/genotyping/
+          results/{wildcards.sample}/hla/mhc-II/genotyping/ > {log} 2>&1
     """
 
 rule merge_predicted_mhcII_allels: