Add variant calling from bam files workflow

.github/workflows/ci.yaml

@@ -83,6 +83,7 @@ jobs:
           - mapping-bwa-mem-1
           - mapping-bwa-mem2-1
           - mapping-bwa-mem2-2
+          - mappings-table
           - pileup-1
           - trimming-adapterremoval-1
           - trimming-adapterremoval-2

.gitignore

@@ -15,6 +15,7 @@ test/vep-cache.zip
 test/prep.log
 test/reference
 test/reference.tar.gz
+test/mappings.tsv
 test/samples.tsv
 test/samples-pe.tsv
 test/out-*

config/config.yaml

@@ -7,11 +7,12 @@
 # Adjust the below dummy file paths and other options as needed.
 data:
 
-  # Needed. Samples table that lists all samples with their units (re-sequencing runs) and fastq files.
+  # Needed for the standard use case of calling variants from `fastq` input files.
+  # The samples table lists all samples with their units (re-sequencing runs) and fastq files.
   # The path to the table file itself, as well as the paths to the fastq files within the table need
   # to be absolute! No relative paths (e.g., `../data`), and no paths using shortcuts such as `~`.
   # See the wiki (https://github.com/lczech/grenepipe/wiki) for the expected format.
-  # The file name needs to end in ".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fas", ".fas.gz".
+  # The file names need to end in ".fasta", ".fasta.gz", ".fa", ".fa.gz", ".fna", ".fas", ".fas.gz".
   samples-table: "/path/to/data/samples.tsv"
 
   # Needed. Path to the reference genome in fasta format, uncompressed, or gzipped.
@@ -48,6 +49,12 @@ data:
   # column are 1, the number of rows in the table corresponds to the number of samples.
   samples-count: 0
 
+  # Advanced usage: Skip all fastq related steps, i.e., the trimming, mapping, and the bam file
+  # processing, and instead start the whole pipeline with a given set of bam files, running only
+  # the variant calling and filtering.
+  # See https://github.com/moiexpositoalonsolab/grenepipe/wiki/Advanced-Usage#running-the-variant-calling-from-a-set-of-bam-files
+  mappings-table: ""
+
   # If the file provided by `reference-genome` above do not already exist, we can automatically
   # download it from ensembl (or any other URL), using the following specification.
   # When downloading from ensembl, try filling in the specificiation of your reference.
@@ -593,6 +600,9 @@ params:
     # processing all (optional) steps (filtering, clipping, dedup, base recalibration).
     # However, we can also run on the raw merged samples, where reads just have been mapped, sorted,
     # and merged per sample (all units of a sample into one bam file), but before any other processing.
+    # When running the pipeline from a set of bam files, that is, when processing a mappings-table
+    # instead of mapping the the fastq files of the `samples-table` in the pipeline itself,
+    # the bam files from that table are used instead, and these settings here are ignored.
     # Valid values: "processed", "merged"
     stats-bams: "processed"
     flagstat-bams: "processed"

example/README.md

@@ -42,6 +42,11 @@ Furthermore, for testing, we provide a simple `regions.bed` file that restricts
 to just some of the regions in the chromosomes. Note that this is mostly an experimental feature
 at the moment (see the `config.yaml`, under `settings: restrict-regions`, for details).
 
+Lastly, the `mapping` directory contains mapped `bam` files of the samples, also for testing
+purposes. This can for instance be used in order to run the pipeline for variant calling
+starting from a set of given `bam` files, instead of mapping `fastq` files first.
+See [here](https://github.com/moiexpositoalonsolab/grenepipe/wiki/Advanced-Usage#running-only-parts-of-the-pipeline) for details on this.
+
 
 Pipeline
 ==============

test/cases/mappings-table.txt

@@ -0,0 +1,2 @@
+samples-count: 0	samples-count: 3
+mappings-table: ""	mappings-table: "#BASEPATH#/test/mappings.tsv"

test/mappings-template.tsv

@@ -0,0 +1,4 @@
+sample	bam
+S1	#BASEPATH#/example/mapping/S1.bam
+S2	#BASEPATH#/example/mapping/S2.bam
+S3	#BASEPATH#/example/mapping/S3.bam

test/run.sh

@@ -41,9 +41,10 @@ if [[ ! -d ./test/reference ]]; then
     cp ./example/regions.bed               ./test/reference/regions.bed
 fi
 
-# Copy the samples table, so that we can change the paths without changing the original,
+# Copy the samples tables, so that we can change the paths without changing the original,
 # We need to use a different sed separator here that cannot occur in paths, to avoid conflict.
-cat ./test/samples_template.tsv | sed "s?#BASEPATH#?${BASEPATH}?g" > ./test/samples.tsv
+cat ./test/samples-template.tsv  | sed "s?#BASEPATH#?${BASEPATH}?g" > ./test/samples.tsv
+cat ./test/mappings-template.tsv | sed "s?#BASEPATH#?${BASEPATH}?g" > ./test/mappings.tsv
 
 # For seqprep, we also want a samples table with only paired-end reads.
 # So let's make a copy and remove the line that contains the single-end file.

workflow/rules/initialize-bam.smk

@@ -0,0 +1,155 @@
+# =================================================================================================
+#     Read Mapping Table
+# =================================================================================================
+
+# Similar setup to what we do in `initialize-fastq.smk`, see there for details.
+if "global" in config:
+    raise Exception("Config key 'global' already defined. Someone messed with our setup.")
+else:
+    config["global"] = {}
+
+# Read mappings table, and enforce to use strings in the index
+config["global"]["samples"] = pd.read_csv(
+    config["data"]["mappings-table"], sep="\t", dtype=str
+).set_index(["sample"], drop=False)
+
+# Get a list of all samples names, in the same order as the input sample table.
+# We cannot use a simple approach here, as this messes up the sample
+# order, which we do not want... (good that we noticed that bug though!)
+# So instead, we iterate, and add sample names incrementally.
+config["global"]["sample-names"] = list()
+for index, row in config["global"]["samples"].iterrows():
+    s = row["sample"]
+    if s not in config["global"]["sample-names"]:
+        config["global"]["sample-names"].append(s)
+
+
+# Helper function for user output to summarize the samples found in the table.
+# In the fastq case, we here print the units as well.
+# For bam files, it's just simply the number of samples.
+def get_sample_units_print():
+    return str(len(config["global"]["sample-names"]))
+
+
+# Wildcard constraints: only allow sample names from the spreadsheet to be used
+wildcard_constraints:
+    sample="|".join(config["global"]["sample-names"]),
+
+
+# =================================================================================================
+#     Sample and File Validity Checks
+# =================================================================================================
+
+
+# Check that the number of samples fits with the expectation, if provided by the user.
+if config["data"].get("samples-count", 0) > 0:
+    if config["data"]["samples-count"] != len(config["global"]["sample-names"]):
+        raise Exception(
+            "Inconsistent number of samples found in the sample table ("
+            + str(len(config["global"]["sample-names"]))
+            + ") and the samples count consistency check provided in the config file under "
+            + "`data: samples-count` ("
+            + str(config["data"]["samples-count"])
+            + ")."
+        )
+
+# We recommend to use absolute paths. Check that for the samples table.
+if not os.path.isabs(config["data"]["mappings-table"]):
+    logger.warning(
+        "Path to the samples table as provided in the config file is not an absolute path. "
+        "We recommend using absolute paths for all files.\n"
+    )
+
+# Run integrity checks similar to what we do for the fastq files.
+uniq_sample_names = list()
+problematic_filenames = 0
+relative_filenames = 0
+for index, row in config["global"]["samples"].iterrows():
+    if row["sample"] in uniq_sample_names:
+        raise Exception(
+            "Multiple rows with identical sample name found in samples table: "
+            + str(row["sample"])
+        )
+    uniq_sample_names.append(row["sample"])
+
+    # Do a check that the sample names are valid file names.
+    # They are used for file names, and would cause weird errors if they contain bad chars.
+    if not valid_filename(row["sample"]):
+        raise Exception(
+            "Invalid sample name name found in samples table that contains characters "
+            + "which cannot be used as sample names for naming output files: "
+            + str(row["sample"])
+            + "; for maximum robustness, we only allow alpha-numerical, dots, dashes, "
+            "and underscores. "
+            + "Use for example the script `tools/copy-samples.py --mode link [...] --clean` "
+            + "to create a new samples table and symlinks to the existing fastq files to solve this."
+        )
+
+    # Do a check of the fastq file names.
+    if not os.path.isfile(row["bam"]):
+        raise Exception(
+            "Input bam files listed in the input files mappings table "
+            + config["data"]["mappings-table"]
+            + " not found: "
+            + str(row["bam"])
+        )
+    if not valid_filepath(row["bam"]):
+        problematic_filenames += 1
+    if not os.path.isabs(row["bam"]):
+        relative_filenames += 1
+
+# Warning about input names and files.
+if problematic_filenames > 0:
+    logger.warning(
+        str(problematic_filenames)
+        + " of the "
+        + str(len(config["global"]["sample-names"]))
+        + " input samples listed in the input files mappings table "
+        + config["data"]["mappings-table"]
+        + " contain problematic characters. We generally advise to only use alpha-numeric "
+        "characters, dots, dashes, and underscores. "
+        "We will try to continue running with these files, but it might lead to errors.\n"
+    )
+if relative_filenames > 0:
+    logger.warning(
+        str(relative_filenames)
+        + " of the "
+        + str(len(config["global"]["sample-names"]))
+        + " input samples listed in the input files mappings table "
+        + config["data"]["mappings-table"]
+        + " use relative file paths. We generally advise to only use absolute paths. "
+        "We will try to continue running with these files, but it might lead to errors.\n"
+    )
+
+
+# Check if a given string can be converted to a number, https://stackoverflow.com/q/354038/4184258
+def is_number(s):
+    try:
+        float(s)  # for int, long, float
+    except ValueError:
+        return False
+    return True
+
+
+# We check if any sample names are only numbers, and warn about this. Bioinformatics is messy...
+numeric_sample_names = 0
+for sn in config["global"]["sample-names"]:
+    if is_number(sn):
+        numeric_sample_names += 1
+
+if numeric_sample_names > 0:
+    logger.warning(
+        str(numeric_sample_names)
+        + " of the "
+        + str(len(config["global"]["sample-names"]))
+        + " input sample names listed in the input files mappings table "
+        + config["data"]["mappings-table"]
+        + " are just numbers, or can be converted to numbers. We generally advise to avoid this, "
+        "as it might confuse downstream processing such as Excel and R. The same applies for names "
+        'that can be converted to dates ("Dec21"), but we do not check this here. '
+        "We will now continue running with these files, as we can work with this here, "
+        "but recommend to change the sample names.\n"
+    )
+
+del uniq_sample_names, problematic_filenames, relative_filenames
+del numeric_sample_names

workflow/rules/initialize-fastq.smk

@@ -39,12 +39,6 @@ config["global"]["unit-names"] = list(
 )
 
 
-# Wildcard constraints: only allow sample names from the spreadsheet to be used
-wildcard_constraints:
-    sample="|".join(config["global"]["sample-names"]),
-    unit="|".join(config["global"]["unit-names"]),
-
-
 # Helper function to get a list of all units of a given sample name.
 def get_sample_units(sample):
     res = list()
@@ -65,6 +59,12 @@ def get_sample_units_print():
         )
 
 
+# Wildcard constraints: only allow sample names from the spreadsheet to be used
+wildcard_constraints:
+    sample="|".join(config["global"]["sample-names"]),
+    unit="|".join(config["global"]["unit-names"]),
+
+
 # =================================================================================================
 #     Sample and File Validity Checks
 # =================================================================================================
@@ -76,7 +76,8 @@ if config["data"].get("samples-count", 0) > 0:
         raise Exception(
             "Inconsistent number of samples found in the sample table ("
             + str(len(config["global"]["sample-names"]))
-            + ") and the samples count consistency check provided in the config file ("
+            + ") and the samples count consistency check provided in the config file under "
+            + "`data: samples-count` ("
             + str(config["data"]["samples-count"])
             + ")."
         )
@@ -124,7 +125,7 @@ for index, row in config["global"]["samples"].iterrows():
         not pd.isnull(row["fq2"]) and not os.path.isfile(row["fq2"])
     ):
         raise Exception(
-            "Input fastq files listed in the input files table "
+            "Input fastq files listed in the input files samples table "
             + config["data"]["samples-table"]
             + " not found: "
             + str(row["fq1"])
@@ -146,7 +147,7 @@ if problematic_filenames > 0:
         str(problematic_filenames)
         + " of the "
         + str(len(config["global"]["sample-names"]))
-        + " input samples listed in the input files table "
+        + " input samples listed in the input files samples table "
         + config["data"]["samples-table"]
         + " contain problematic characters. We generally advise to only use alpha-numeric "
         "characters, dots, dashes, and underscores. "
@@ -159,15 +160,14 @@ if relative_filenames > 0:
         str(relative_filenames)
         + " of the "
         + str(len(config["global"]["sample-names"]))
-        + " input samples listed in the input files table "
+        + " input samples listed in the input files samples table "
         + config["data"]["samples-table"]
         + " use relative file paths. We generally advise to only use absolute paths. "
         "Use for example the script `tools/generate-table.py` "
         "to create a samples table with absolute paths. "
         "We will try to continue running with these files, but it might lead to errors.\n"
     )
 
-del problematic_filenames, relative_filenames
 
 
 # Check if a given string can be converted to a number, https://stackoverflow.com/q/354038/4184258
@@ -190,7 +190,7 @@ if numeric_sample_names > 0:
         str(numeric_sample_names)
         + " of the "
         + str(len(config["global"]["sample-names"]))
-        + " input sample names listed in the input files table "
+        + " input sample names listed in the input files samples table "
         + config["data"]["samples-table"]
         + " are just numbers, or can be converted to numbers. We generally advise to avoid this, "
         "as it might confuse downstream processing such as Excel and R. The same applies for names "
@@ -199,4 +199,5 @@ if numeric_sample_names > 0:
         "but recommend to change the sample names.\n"
     )
 
+del problematic_filenames, relative_filenames
 del numeric_sample_names

workflow/rules/initialize.smk

@@ -34,10 +34,14 @@ snakemake.utils.validate(config, schema="../schemas/config.schema.yaml")
 report: os.path.join(workflow.basedir, "reports/workflow.rst")
 
 
-# Include the functions neeed to initialize the pipeline for analysing a set of fastq samples.
-# This reads the samples table, and provides validation and user output functions for it.
+# Include the functions neeed to initialize the pipeline for analysing a set of fastq samples,
+# or, if the mappings table is given, starting from there.
+# This reads the samples/mappings table, and provides validation and user output functions for it.
 include: "initialize-reference.smk"
-include: "initialize-fastq.smk"
+if "mappings-table" in config["data"] and config["data"]["mappings-table"]:
+    include: "initialize-bam.smk"
+else:
+    include: "initialize-fastq.smk"
 
 
 # =================================================================================================