Bring the MAS-ISO-seq workflows into our production pipelines

.dockstore.yml

@@ -100,3 +100,11 @@ workflows:
   subclass: wdl
   primaryDescriptorPath: /wdl/ONTPfTypeDrugResistanceMarkers.wdl
   testParameterFiles:
+- name: PBMASIsoSeqQuantify
+  subclass: wdl
+  primaryDescriptorPath: /wdl/PBMASIsoSeqQuantify.wdl
+  testParameterFiles:
+- name: PBMASIsoSeqDemultiplex
+  subclass: wdl
+  primaryDescriptorPath: /wdl/PBMASIsoSeqDemultiplex.wdl
+  testParameterFiles:

docker/lr-10x/Dockerfile

@@ -5,10 +5,9 @@ RUN apt-get install -y --no-install-recommends git ssh ca-certificates autoconf
     && rm -rf /var/lib/apt/lists/*
 
 # install htslib
-RUN git clone https://github.com/samtools/htslib.git \
-    && cd htslib \
-    && autoheader \
-    && autoconf \
+RUN wget https://github.com/samtools/htslib/releases/download/1.11/htslib-1.11.tar.bz2 \
+    && tar -jxf htslib-1.11.tar.bz2 \
+    && cd htslib-1.11 \
     && ./configure \
     && make \
     && make install
@@ -26,6 +25,27 @@ WORKDIR /lrma
 RUN wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O /miniconda.sh \
     && bash /miniconda.sh -b -p /miniconda
 
+# Install samtools 1.11:
+RUN apt-get update
+RUN apt-get install -y bzip2 curl gnupg2
+RUN apt-get install -y libc-dev ncurses-dev
+
+#### Specific for google cloud support
+RUN echo "deb [signed-by=/usr/share/keyrings/cloud.google.gpg] http://packages.cloud.google.com/apt cloud-sdk main" | tee -a /etc/apt/sources.list.d/google-cloud-sdk.list \
+		&& curl https://packages.cloud.google.com/apt/doc/apt-key.gpg | apt-key --keyring /usr/share/keyrings/cloud.google.gpg  add - \
+		&& apt-get update -y \
+		&& apt-get install google-cloud-sdk -y
+
+RUN apt-get install -y libcurl4-openssl-dev libssl-dev
+
+# Get samtools source:
+RUN wget https://github.com/samtools/samtools/releases/download/1.11/samtools-1.11.tar.bz2 \
+		&& tar -xjf samtools-1.11.tar.bz2 \
+		&& cd samtools-1.11 \
+		&& ./configure \
+		&& make install
+RUN rm -rf /samtools-1.11 /samtools-1.11.tar.bz2
+
 # install conda packages
 ADD ./environment.yml /lrma/environment.yml
 RUN pip install --upgrade pip
@@ -58,8 +78,22 @@ ADD reverse_adapter_sequence.fasta.bwt /lrma/reverse_adapter_sequence.fasta.bwt
 ADD reverse_adapter_sequence.fasta.pac /lrma/reverse_adapter_sequence.fasta.pac
 ADD reverse_adapter_sequence.fasta.sa /lrma/reverse_adapter_sequence.fasta.sa
 
+RUN pip3 install pysam biopython
+RUN pip3 install tqdm
+
 ADD tool.py /lrma/tool.py
 ADD tool_rle.py /lrma/tool_rle.py
+ADD tool_starcode_seeded.py /lrma/tool_starcode_seeded.py
+ADD extract_ilmn_bc_conf_scores.py /lrma/extract_ilmn_bc_conf_scores.py
+
+ADD restore_annotations_to_aligned_bam.py /lrma/restore_annotations_to_aligned_bam.py
+ADD extract_cbc_and_umi_from_annotated_read.py /lrma/extract_cbc_and_umi_from_annotated_read.py 
+ADD copy_contig_name_to_tag.py /lrma/copy_contig_name_to_tag.py 
+ADD tag_mas_sirv_umi_positions.py /lrma/tag_mas_sirv_umi_positions.py
+ADD update_umi_positions.py /lrma/update_umi_positions.py
+ADD update_umi_positions_2.py /lrma/update_umi_positions_2.py
+
+RUN echo -n "set number\nsyntax on\nset hlsearch" > ~/.vimrc
 
 RUN echo "source activate 10x_tool" > ~/.bashrc
-RUN pip3 install pysam biopython
+

docker/lr-10x/Makefile

@@ -1,16 +1,17 @@
-VERSION = 0.1.9
-TAG1 = quay.io/broad-long-read-pipelines/lr-10x:$(VERSION)
-TAG2 = quay.io/broad-long-read-pipelines/lr-10x:latest
-TAG3 = us.gcr.io/broad-dsp-lrma/lr-10x:$(VERSION)
-TAG4 = us.gcr.io/broad-dsp-lrma/lr-10x:latest
+VERSION = 0.1.18
+NAME = lr-10x
 
-all: build push
+TAG1 = us.gcr.io/broad-dsp-lrma/$(NAME):$(VERSION)
+TAG2 = us.gcr.io/broad-dsp-lrma/$(NAME):latest
+
+all: | build push
 
 build:
-	docker build -t $(TAG1) -t $(TAG2) -t $(TAG3) -t $(TAG4) .
+	docker build -t $(TAG1) -t $(TAG2) .
+
+build_no_cache:
+	docker build --no-cache -t $(TAG1) -t $(TAG2) .
 
 push:
 	docker push $(TAG1)
 	docker push $(TAG2)
-	docker push $(TAG3)
-	docker push $(TAG4)

docker/lr-10x/copy_contig_name_to_tag.py

@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+import argparse
+
+import pysam
+import tqdm
+
+gGENE_TAG_NAME = "XG"
+
+
+def main(input_bam, out_bam_name, gene_tag_name):
+
+    print("Verifying input files exist...")
+    files_ok = True
+    for f in [input_bam]:
+        if not os.path.exists(f):
+            print(f"ERROR: Input file does not exist: {f}")
+            files_ok = False
+    if not files_ok:
+        sys.exit(1)
+    print("Input files verified.")
+
+    reads_processed = 0
+
+    pysam.set_verbosity(0)  # silence message about the .bai file not being found
+    with pysam.AlignmentFile(
+            input_bam, "rb", check_sq=False, require_index=False
+    ) as bam_file, tqdm.tqdm(
+        desc="Progress",
+        unit=" read",
+        file=sys.stderr
+    ) as pbar:
+
+        # Get our header from the input bam file:
+        out_bam_header_dict = bam_file.header.to_dict()
+
+        # Add our program group to it:
+        pg_dict = {
+            "ID": f"copy-contig-name-to-tag-0.0.1",
+            "PN": "copy-contig-name-to-tag",
+            "VN": f"0.0.1",
+            "DS": "Copies the name of the contig to which a read is aligned to a tag within that read.",
+            "CL": " ".join(sys.argv),
+        }
+        if "PG" in out_bam_header_dict:
+            out_bam_header_dict["PG"].append(pg_dict)
+        else:
+            out_bam_header_dict["PG"] = [pg_dict]
+        out_header = pysam.AlignmentHeader.from_dict(out_bam_header_dict)
+
+        # Open our output file so we can write to it:
+        with pysam.AlignmentFile(out_bam_name, "wb", header=out_header) as out_bam_file:
+            for read in bam_file:
+
+                # Set our tag for our contig:
+                read.set_tag(gene_tag_name, read.reference_name)
+
+                # Write out our read:
+                out_bam_file.write(read)
+
+                pbar.update(1)
+                reads_processed += 1
+    print("Done!")
+    print(f"Reads processed: {reads_processed}")
+
+
+if __name__ == "__main__":
+
+    parser = argparse.ArgumentParser(
+        description=f"Copies the name of the contig to which a read is aligned to a tag within that read.  "
+                    f"(tag name: {gGENE_TAG_NAME})",
+    )
+
+    requiredNamed = parser.add_argument_group('required named arguments')
+    requiredNamed.add_argument('-b', '--bam',
+                               help='Aligned bam file from which to extract the contig information into a tag.',
+                               required=True)
+    requiredNamed.add_argument('-o', '--out-name',
+                               help='Output bam file name',
+                               required=True)
+
+    args = parser.parse_args()
+    main(args.bam, args.out_name, gGENE_TAG_NAME)

docker/lr-10x/environment.yml

@@ -7,10 +7,10 @@ dependencies:
   - python=3.6
   - cython
   - numpy
-  - samtools
   - bwapy
   - cigar
   - pip
   - pip:
     - pysam
     - biopython
+  - tqdm