feat: Allow bam input files (#94)

This PR allows bam files as input via the fields `bam_single` or
`bam_paired` in `units.tsv`.

---------

Co-authored-by: Johannes Köster <[email protected]>

.test/config/units.tsv

@@ -1,4 +1,4 @@
-sample unit fragment_len_mean fragment_len_sd fq1 fq2
+sample unit fragment_len_mean fragment_len_sd fq1 fq2 bam_single bam_paired
 A 1 ngs-test-data/reads/a.chr21.1.fq ngs-test-data/reads/a.chr21.2.fq
 B 1 ngs-test-data/reads/b.chr21.1.fq ngs-test-data/reads/b.chr21.2.fq
 B 2 300 14 ngs-test-data/reads/b.chr21.1.fq 

.test/three_prime/config/units.tsv

@@ -1,4 +1,4 @@
-sample unit fragment_len_mean fragment_len_sd fq1 fq2
+sample unit fragment_len_mean fragment_len_sd fq1 fq2 bam_single bam_paired
 SRR8309096 u1 430 43 quant_seq_test_data/SRR8309096.fastq.gz 
 SRR8309094 u1 430 43 quant_seq_test_data/SRR8309094.fastq.gz 
 SRR8309095 u1 430 43 quant_seq_test_data/SRR8309095.fastq.gz 

config/units.tsv

@@ -1,4 +1,4 @@
-sample unit fragment_len_mean fragment_len_sd fq1 fq2
+sample unit fragment_len_mean fragment_len_sd fq1 fq2 bam_single bam_paired
 A 1 raw/a.chr21.1.fq raw/a.chr21.2.fq
 B 1 raw/b.chr21.1.fq raw/b.chr21.2.fq
 B 2 300 14 raw/b.chr21.1.fq 

workflow/Snakefile

@@ -25,6 +25,7 @@ include: "rules/diffexp.smk"
 include: "rules/diffsplice.smk"
 include: "rules/enrichment.smk"
 include: "rules/datavzrd.smk"
+include: "rules/bam.smk"
 
 
 rule all:

workflow/rules/bam.smk

@@ -0,0 +1,33 @@
+rule bam_paired_to_fastq:
+ input:
+ lookup(
+ query="sample == '{sample}' & unit == '{unit}'",
+ within=units,
+ cols="bam_paired",
+ ),
+ output:
+ "results/fastq/{sample}-{unit}.1.fq.gz",
+ "results/fastq/{sample}-{unit}.2.fq.gz",
+ log:
+ "logs/fastq/{sample}-{unit}.separate.log",
+ params:
+ fastq="-n",
+ threads: 3
+ wrapper:
+ "v3.10.2/bio/samtools/fastq/separate"
+
+
+rule bam_single_to_fastq:
+ input:
+ lookup(
+ query="sample == '{sample}' & unit == '{unit}'",
+ within=units,
+ cols="bam_single",
+ ),
+ output:
+ "results/fastq/{sample}-{unit}.fq.gz",
+ log:
+ "logs/fastq/{sample}-{unit}.interleaved.log",
+ threads: 3
+ wrapper:
+ "v3.10.2/bio/samtools/fastq/interleaved"

workflow/rules/common.smk

@@ -92,13 +92,20 @@ def get_model(wildcards):
 
 def is_single_end(sample, unit):
  """Determine whether unit is single-end."""
+ bam_paired_not_present = pd.isnull(units.loc[(sample, unit), "bam_paired"])
  fq2_not_present = pd.isnull(units.loc[(sample, unit), "fq2"])
- return fq2_not_present
+ return fq2_not_present and bam_paired_not_present
 
 
 def get_fastqs(wildcards):
  """Get raw FASTQ files from unit sheet."""
- if is_single_end(wildcards.sample, wildcards.unit):
+ if not pd.isnull(units.loc[(wildcards.sample, wildcards.unit), "bam_single"]):
+ return f"results/fastq/{wildcards.sample}-{wildcards.unit}.fq.gz"
+ elif not pd.isnull(units.loc[(wildcards.sample, wildcards.unit), "bam_paired"]):
+ fqfrombam1 = f"results/fastq/{wildcards.sample}-{wildcards.unit}.1.fq.gz"
+ fqfrombam2 = f"results/fastq/{wildcards.sample}-{wildcards.unit}.2.fq.gz"
+ return [fqfrombam1, fqfrombam2]
+ elif is_single_end(wildcards.sample, wildcards.unit):
  return units.loc[(wildcards.sample, wildcards.unit), "fq1"]
  else:
  u = units.loc[(wildcards.sample, wildcards.unit), ["fq1", "fq2"]].dropna()

workflow/rules/quant.smk

@@ -1,8 +1,10 @@
 rule kallisto_index:
  input:
- fasta="resources/transcriptome.cdna.without_poly_a.fasta"
- if is_3prime_experiment
- else "resources/transcriptome.cdna.fasta",
+ fasta=(
+ "resources/transcriptome.cdna.without_poly_a.fasta"
+ if is_3prime_experiment
+ else "resources/transcriptome.cdna.fasta"
+ ),
  output:
  index="results/kallisto_cdna/transcripts.cdna.idx",
  log:

workflow/schemas/units.schema.yaml

@@ -10,11 +10,16 @@ properties:
  description: unit id
  fq1:
  type: string
- description: path to FASTQ file
+ description: path to FASTQ file (leave empty in case usage of bam_single or bam_paired)
  fq2:
  type: string
- description: path to second FASTQ file (leave empty in case of single-end)
+ description: path to second FASTQ file (leave empty in case of single-end or usage of bam_single or bam_paired)
+ bam_single:
+ type: string
+ description: path to single bam file (leave empty in case of usage of fastq files)
+ bam_paired:
+ type: string
+ description: path to paired bam file (leave empty in case of usage of fastq files)
 required:
  - sample
  - unit
- - fq1