Remove constitutive k-mers

.github/workflows/ci.yml

@@ -61,18 +61,23 @@ jobs:
       NXF_VER: '20.07.1'
       NXF_ANSI_LOG: false
     strategy:
+      fail-fast: false
       matrix:
         profile_flags:
           - "test --sketch_scaled false --sketch_scaled_log2 2"
           - "test --sketch_scaled false --sketch_num_hashes 20"
           - "test --sketch_scaled false --sketch_num_hashes_log2 20"
+          - "test_bam"
           - "test_bam --barcodes_file false --rename_10x_barcodes false --save_fastas false --write_barcodes_meta_csv false"
           - "test_bam --rename_10x_barcodes false --write_barcodes_meta_csv false"
           - "test_bam --skip_sig_merge"
           - "test_bam --write_barcodes_meta_csv false"
           - "test_bam --barcodes_file false --rename_10x_barcodes false"
           - "test_bam --rename_10x_barcodes false"
           - "test_fastas"
+          - "test_constitutive_from_download_refseq"
+          - "test_constitutive_from_fasta"
+          - "test_constitutive_from_sig"
           - "test_protein_fastas"
           - "test_remove_ribo"
           - "test_sig_merge"

Dockerfile

@@ -4,14 +4,25 @@ LABEL authors="Olga Botvinnik" \
 
 # Install the conda environment
 COPY environment.yml /
-RUN conda env create --quiet -f /environment.yml && conda clean -a
+RUN conda install -c conda-forge mamba
+RUN mamba env create -f /environment.yml && mamba clean -a
 
 # Add conda installation dir to PATH (instead of doing 'conda activate')
 ENV PATH /opt/conda/envs/nf-core-kmermaid-0.1.0dev/bin:$PATH
 
 # Dump the details of the installed packages to a file for posterity
 RUN conda env export --name nf-core-kmermaid-0.1.0dev > nf-core-kmermaid-0.1.0dev.yml
 
+# Install super fast rust code to remove nuisance hashes (e.g. ribosomal) from signatures
+RUN git clone -b olgabot/mut-warning https://github.com/olgabot/2021-01-27-olga-remove-protein.git
+# Soft link all conda C-related libraries to their non-prefixed name 
+# for rust to be able to build the C libraries
+RUN for f in $(ls /opt/conda/envs/nf-core-kmermaid-0.1.0dev/bin/x86_64-conda_cos6-linux-gnu*); \
+      do g=$(echo $f | sed 's:x86_64-conda_cos6-linux-gnu-::') ; echo $g; ln -s $f $g ; done
+RUN cd 2021-01-27-olga-remove-protein && cargo build --release 
+# Add "subtract" command to path
+ENV PATH $HOME/2021-01-27-olga-remove-protein/target/release:$PATH
+
 # Instruct R processes to use these empty files instead of clashing with a local version
 RUN touch .Rprofile
 RUN touch .Renviron

bin/filter_fasta_regex.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python
+
+import argparse
+import re
+
+
+import screed
+
+
+def write_records_to_fasta(records, fasta):
+    with open(fasta, "w") as f:
+        for record in records:
+            f.write(f'>{record["name"]}\n{record["sequence"]}\n')
+
+
+def filter_records(fasta, pattern):
+    filtered_records = []
+    with screed.open(fasta) as records:
+        for record in records:
+            name = record["name"]
+            if re.findall(pattern, name, flags=re.I):
+                filtered_records.append(record)
+    return filtered_records
+
+
+def filter_fasta_with_regex(fasta_to_filter, out_fasta, regex):
+    record_subset = filter_records(fasta_to_filter, regex)
+    write_records_to_fasta(record_subset, out_fasta)
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(
+        description="""Extract sequences whose names match a pattern"""
+    )
+    parser.add_argument("--input-fasta", type=str, help="Sequence file to filter")
+    parser.add_argument("--output-fasta", type=str, help="File to write")
+    parser.add_argument(
+        "--regex-pattern",
+        type=str,
+        help="Regular expression pattern to match for the names of seuqences in the file",
+    )
+    args = parser.parse_args()
+
+    filter_fasta_with_regex(args.input_fasta, args.output_fasta, args.regex_pattern)

bin/scrape_software_versions.py

@@ -14,8 +14,10 @@
     "SKA": ["v_ska.txt", r"SKA Version: (\S+)"],
     "htslib": ["v_samtools.txt", r"htslib (\S+)"],
     "Sourmash": ["v_sourmash.txt", r"sourmash (\S+)"],
-    "SortMeRNA": ["v_sortmerna.txt", r"SortMeRNA version (\S+),"],
+    "Rsync": ["v_rsync.txt", r"rsync  version (\S+)"],
+    "Rsync (Protocol)": ["v_rsync.txt", r"protocol version (\S+)"],
     "orpheum": ["v_orpheum.txt", r"Version: (\S+)"],
+    "Python": ["v_python.txt", r"Python (\S+)"],
 }
 results = OrderedDict()
 results["nf-core/kmermaid"] = '<span style="color:#999999;">N/A</span>'
@@ -25,11 +27,13 @@
 results["bam2fasta"] = '<span style="color:#999999;">N/A</span>'
 results["fastp"] = '<span style="color:#999999;">N/A</span>'
 results["htslib"] = '<span style="color:#999999;">N/A</span>'
+results["orpheum"] = '<span style="color:#999999;">N/A</span>'
+results["Python"] = '<span style="color:#999999;">N/A</span>'
+results["Rsync"] = '<span style="color:#999999;">N/A</span>'
+results["Rsync (Protocol)"] = '<span style="color:#999999;">N/A</span>'
 results["Samtools"] = '<span style="color:#999999;">N/A</span>'
 results["SKA"] = '<span style="color:#999999;">N/A</span>'
 results["Sourmash"] = '<span style="color:#999999;">N/A</span>'
-results["SortMeRNA"] = '<span style="color:#999999;">N/A</span>'
-results["orpheum"] = '<span style="color:#999999;">N/A</span>'
 
 # Search each file using its regex
 for k, v in regexes.items():

bin/validate_sketch_value.py

@@ -20,7 +20,7 @@ def get_sketch_value(value, value_log2):
             if "," in value:
                 logger.exception(
                     f"Can only provide a single number to --sketch_num_hashes or"
-                    f" --sketch_scaled. Provided '{value}"
+                    f" --sketch_scaled. Provided '{value}'"
                 )
             sketch_value = int(value)
         else:

conf/base.config

@@ -54,6 +54,7 @@ process {
 
   withName: 'multiqc|get_software_versions' {
     memory = { check_max( 2.GB * task.attempt, 'memory' ) }
+    errorStrategy = "ignore"
     cache = false
   }
   withName: 'sourmash_compute_sketch_fastx_nucleotide|sourmash_compute_sketch_fastx_peptide' {

conf/test.config

@@ -17,7 +17,6 @@ params {
   // Input data
   // samples = 'testing/samples.csv'
   // fastas = 'testing/fastas/*.fasta'
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   // read_pairs = 'testing/fastqs/*{1,2}.fastq.gz'
   // sra = "SRP016501"
@@ -29,4 +28,8 @@ params {
     ['SRR4238351', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238351_subsamp.fastq.gz']],
     ['SRR4238355', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238355_subsamp.fastq.gz']],
   ]
+  // Remove constitutively expressed genes
+  test_mini_refseq_download = true
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_bam.config

@@ -19,7 +19,6 @@ params {
     'https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid--bam-unique-names/testdata/mouse_lung.bam',
     'https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid--bam-unique-names/testdata/mouse_brown_fat_ptprc_plus_unaligned.bam']
   // Sketch Parameters
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   read_pairs = false
   save_fastas = "fastas"
@@ -28,4 +27,6 @@ params {
   // For bam, each fasta record represents each barcode and each should have a signature
   // they should not be merged, For computation on bam file using sourmash, please set true for the below flag
   tenx_min_umi_per_cell = 2
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_constitutive_from_download_refseq.config

@@ -0,0 +1,32 @@
+/*
+ * -------------------------------------------------
+ *  Nextflow config file for running tests
+ * -------------------------------------------------
+ * Defines bundled input files and everything required
+ * to run a fast and simple test. Use as follows:
+ *   nextflow run nf-core/kmermaid -profile test
+ */
+
+params {
+  config_profile_name = 'Test profile'
+  config_profile_description = 'Minimal test dataset to check pipeline function'
+  // Limit resources so that this can run on Travis
+  max_cpus = 2
+  max_memory = 6.GB
+  max_time = 48.h
+
+  // Input data
+  input_paths = [
+    ['SRR4050379', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_2.fastq.gz']],
+    ['SRR4050380', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_2.fastq.gz']],
+    ['SRR4238351', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238351_subsamp.fastq.gz']],
+    ['SRR4238355', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238355_subsamp.fastq.gz']],
+  ]
+
+  // "Other" is the smallest refseq taxonomy subdirectory: ftp://ftp.ncbi.nlm.nih.gov/refseq/release/other/
+  // Protein fasta is 453 B
+  refseq_taxonomy = 'vertebrate_mammalian'
+  test_mini_refseq_download = true
+}

conf/test_constitutive_from_fasta.config

@@ -0,0 +1,32 @@
+/*
+ * -------------------------------------------------
+ *  Nextflow config file for running tests
+ * -------------------------------------------------
+ * Defines bundled input files and everything required
+ * to run a fast and simple test. Use as follows:
+ *   nextflow run nf-core/kmermaid -profile test
+ */
+
+params {
+  config_profile_name = 'Test profile'
+  config_profile_description = 'Minimal test dataset to check pipeline function'
+  // Limit resources so that this can run on Travis
+  max_cpus = 2
+  max_memory = 6.GB
+  max_time = 48.h
+
+  // Input data
+  input_paths = [
+    ['SRR4050379', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_2.fastq.gz']],
+    ['SRR4050380', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_2.fastq.gz']],
+    ['SRR4238351', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238351_subsamp.fastq.gz']],
+    ['SRR4238355', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238355_subsamp.fastq.gz']],
+  ]
+  constitutive_protein_fasta = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa.gz"
+  constitutive_rna_fasta = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa.gz"
+
+  translate_proteome_fasta = 'https://github.com/nf-core/test-datasets/raw/kmermaid/reference/gencode.v32.pc_translations.subsample5.randomseed0.fa'
+  bloomfilter_tablesize = '1e6'
+}

conf/test_constitutive_from_sig.config

@@ -0,0 +1,29 @@
+/*
+ * -------------------------------------------------
+ *  Nextflow config file for running tests
+ * -------------------------------------------------
+ * Defines bundled input files and everything required
+ * to run a fast and simple test. Use as follows:
+ *   nextflow run nf-core/kmermaid -profile test
+ */
+
+params {
+  config_profile_name = 'Test profile'
+  config_profile_description = 'Minimal test dataset to check pipeline function'
+  // Limit resources so that this can run on Travis
+  max_cpus = 2
+  max_memory = 6.GB
+  max_time = 48.h
+
+  // Input data
+  input_paths = [
+    ['SRR4050379', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050379_pass_2.fastq.gz']],
+    ['SRR4050380', ['https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_1.fastq.gz',
+                    'https://github.com/nf-core/test-datasets/raw/kmermaid/testdata/SRR4050380_pass_2.fastq.gz']],
+    ['SRR4238351', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238351_subsamp.fastq.gz']],
+    ['SRR4238355', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238355_subsamp.fastq.gz']],
+  ]
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+}

conf/test_fastas.config

@@ -17,7 +17,6 @@ params {
   // Input data
   // samples = 'testing/samples.csv'
   // fastas = 'testing/fastas/*.fasta'
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   // read_pairs = 'testing/fastqs/*{1,2}.fastq.gz'
   // sra = "SRP016501"
@@ -26,4 +25,6 @@ params {
     ['SRR4050380', ['https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid--bam-unique-names/testdata/SRR4050380_pass_concatenated.fasta']],
 
   ]
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_full.config

@@ -12,10 +12,12 @@ params {
   config_profile_description = 'Full test dataset to check pipeline function'
 
   // Input data for full size test
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   input_paths = [
     ['GM12878', ['ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/007/SRR3192657/SRR3192657_1.fastq.gz,ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/007/SRR3192657/SRR3192657_2.fastq.gz','ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/008/SRR3192658/SRR3192658_1.fastq.gz,ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/008/SRR3192658/SRR3192658_2.fastq.gz']],
     ['K562', ['ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/008/SRR3192408/SRR3192408_1.fastq.gz,ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/008/SRR3192408/SRR3192408_2.fastq.gz', 'ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/009/SRR3192409/SRR3192409_1.fastq.gz,ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR319/009/SRR3192409/SRR3192409_2.fastq.gz']]
   ]
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+
 }

conf/test_protein_fastas.config

@@ -26,8 +26,8 @@ params {
       ['https://github.com/czbiohub/test-datasets/raw/predictorthologs/testdata/bonobo_liver_ptprc__molecule-dayhoff__coding_reads_peptides.fasta']]]
 
   // Sketch Parameters
-  sketch_scaled = 2
   molecules = 'protein,dayhoff,hp'
   read_pairs = false
-
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_remove_ribo.config

@@ -17,7 +17,6 @@ params {
   // Input data
   // samples = 'testing/samples.csv'
   // fastas = 'testing/fastas/*.fasta'
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   // read_pairs = 'testing/fastqs/*{1,2}.fastq.gz'
   // sra = "SRP016501"
@@ -31,4 +30,6 @@ params {
     ['SRR4238351', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238351_subsamp.fastq.gz']],
     ['SRR4238355', ['https://github.com/nf-core/test-datasets/raw/rnaseq/testdata/SRR4238355_subsamp.fastq.gz']],
   ]
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_sig_merge.config

@@ -18,7 +18,6 @@ params {
   bam = ['https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid--bam-unique-names/testdata/mouse_lung.bam',
     'https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid--bam-unique-names/testdata/mouse_brown_fat_ptprc_plus_unaligned.bam']
   // Sketch Parameters
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   read_pairs = false
   save_fastas = "fastas"
@@ -29,4 +28,6 @@ params {
 
   reference_proteome_fasta = 'https://github.com/nf-core/test-datasets/raw/kmermaid/reference/ptprc_bam_translations.fa'
   bloomfilter_tablesize = '1e6'
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }

conf/test_tenx_tgz.config

@@ -20,14 +20,13 @@ params {
     'https://github.com/nf-core/test-datasets/raw/olgabot/kmermaid-unaligned-tgz-v3/testdata/mouse_brown_fat_ptprc_plus_unaligned.tgz'
   ]
   // Sketch Parameters
-  sketch_scaled = 2
   molecules = 'dna,protein,dayhoff'
   read_pairs = false
   save_fastas = "fastas"
   save_intermediate_files = "/tmp/"
   write_barcode_meta_csv = "metadata.csv"
   // For bam, each fasta record represents each barcode and each should have a signature
-  // they should not be merged, For computation on bam file using sourmash, please set true for the below flag
   tenx_min_umi_per_cell = 10
-  shard_size = 350
+  constitutive_protein_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.protein.fa__only_constitutive_genes.fa__molecule-protein,dayhoff__ksize-21,30,51__scaled-10__track_abundance-true.sig"
+  constitutive_rna_sig = "https://github.com/czbiohub/test-datasets/raw/olgabot/kmermaid--housekeeping-fasta/reference/vertebrate_mammalian--205--2021-03-15.rna.fa__only_constitutive_genes.fa__molecule-dna__ksize-21,30,51__scaled-10__track_abundance-true.sig"
 }