Readgroup issue (#28)

* catch errors when rnaseq is not provided (exitron/altsplicing calling) and add reference genome index as input for indel calling (htcaller)

* added reference genome index on germline indel calling (which is required when only indel calling has been activated & remove -C from BWA mem call (on DNAseq data) which causes issues on Illumina identifier

CHANGELOG.md

@@ -16,6 +16,14 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Prioritization of neoantigens is now done separately for each variant type (speeds up the process)
 - NMD information (e.g., escape rule,...) is now also calculated for all variants
 
+## [0.2.6] - 2024-06-20
+
+### Fix 
+
+- Added routines to catch errors when rnaseq data is not provided but exitron/alternative splicing calling is activated
+- Added reference genome index as input to germline indel calling (necessary when only indel calling is activated)
+- removed -C from BWA mem call (on DNAseq data) to avoid error on Illumina identifiers
+
 ## [0.2.5] - 2024-06-19
 
 ### Fix 

workflow/rules/align.smk

@@ -231,10 +231,9 @@ if config['data']['dnaseq_filetype'] in ['.fq','.fastq']:
     threads: config['threads']
     shell:
       """
-        bwa mem -t{threads} -C resources/refs/bwa/genome {input.reads} \
-        | samtools addreplacerg -r ID:{wildcards.group} -r SM:{wildcards.sample} \
-        -r LB:{wildcards.sample} -r PL:ILLUMINA -r PU:{wildcards.group} - - \
-        | samtools sort -@ 6 -n -m1g - -o {output} > {log} 2>&1
+        bwa mem -t{threads} resources/refs/bwa/genome \
+            -R '@RG\\tID:{wildcards.group}\\tSM:{wildcards.sample}\\tLB:{wildcards.sample}\\tPL:ILLUMINA' \
+            {input.reads} | samtools sort -@ 6 -n -m1g - -o {output} > {log} 2>&1
       """
 
   rule dnaseq_postproc:

workflow/rules/common.smk

@@ -734,18 +734,22 @@ def get_prioritization_long_indels(wildcards):
 def get_prioritization_exitrons(wildcards):
   exitrons = []
   if config["exitronsplicing"]["activate"]:
-    exitrons += expand("results/{sample}/annotation/exitrons.vcf",
-                       sample=config["data"]["name"])
-
+    if len(config["data"]["rnaseq"]) != 0:
+      exitrons += expand("results/{sample}/annotation/exitrons.vcf",
+                         sample=config["data"]["name"])
+    else:
+      print('rnaseq data has not been specified in the config file, but exitron calling is activated - skipping...')
   return exitrons
 
 
 def get_prioritization_altsplicing(wildcards):
   altsplicing = []
   if config["altsplicing"]["activate"]:
-    altsplicing += expand("results/{sample}/annotation/altsplicing.vcf",
-                          sample=config["data"]["name"])
-
+    if len(config["data"]["rnaseq"]) != 0:
+      altsplicing += expand("results/{sample}/annotation/altsplicing.vcf",
+                            sample=config["data"]["name"])
+    else:
+      print('rnaseq data has not been specified in the config file, but alternative splicing calling is activated - skipping...')
   return altsplicing
 
 def get_prioritization_custom(wildcards):

workflow/rules/germline.smk

@@ -73,7 +73,8 @@ rule detect_variants_htc_first_round:
   input:
     bam="results/{sample}/{seqtype}/align/{group}_final_BWA_split/{chr}.bam",
     idx="results/{sample}/{seqtype}/align/{group}_final_BWA_split/{chr}.bam.bai",
-    ref="resources/refs/genome.fasta"
+    ref="resources/refs/genome.fasta",
+    ref_idx="resources/refs/genome.fasta.fai"
   output:
     vcf="results/{sample}/{seqtype}/indel/htcaller/{group}_variants.1rd/{chr}.vcf"
   message: