Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418

SHuang-Broad · 2023-08-18T15:45:30Z

Project-specific.

SHuang-Broad · 2023-08-18T15:46:42Z

docker/lr-custom-gatk/Dockerfile

+           /usr/share/man/??_* && \
+    samtools --version
+
+    # echo "deb [signed-by=/usr/share/keyrings/cloud.google.gpg] https://packages.cloud.google.com/apt cloud-sdk main" | tee -a /etc/apt/sources.list.d/google-cloud-sdk.list && \


Note to self: take this away in final push

SHuang-Broad · 2023-08-18T15:48:29Z

docker/lr-gcloud-samtools/Makefile

@@ -0,0 +1,12 @@
+VERSION = 0.1.1


Down the road, the intention is to replace lr-basic with this image.
You'll see here that the only meaningful difference with lr-basic is it installs google-cloud-cli not google-cloud-sdk.

Can we eliminate lr-basic now?

tedsharpe

There are a lot of suggestions, but they're mostly pretty trivial. Adopt the ones that seem in good taste to you.

tedsharpe · 2023-10-17T18:45:35Z

.dockstore.yml

@@ -1,17 +1,11 @@
 version: 1.2


I haven't looked at this in detail. I'm going to assume it's OK.

tedsharpe · 2023-10-17T18:47:53Z

docker/lr-bam-dedup/Dockerfile

@@ -0,0 +1,73 @@
+FROM python:3.9.16-slim-bullseye


These dockerfiles are so long, and have so much repeated content.
Wouldn't it be possible to start with your lr-gcloud-samtools image as a first stage?

tedsharpe · 2023-10-17T18:50:57Z

docker/lr-bam-dedup/Dockerfile

+                 libcurl4-openssl-dev \
+                 liblzma-dev \
+                 libncurses5-dev \
+                 autoconf \


Note that all this stuff -- automake, gcc, wget, etc. -- is left in the final image.
It would be great to leave all that behind in an earlier stage. (See above comment about using lr-basic as a first stage, and do all the compilation there.)

I now see that you remove most of it below. I still think it would be cleaner to remove it via staging.

tedsharpe · 2023-10-17T18:58:13Z

docker/lr-bam-dedup/remove_duplicate_ont_aln.py

+                sam_flag = read.flag
+                pos = read.reference_start
+                cigar = read.cigarstring
+                signature = f"{n}-{chrom}-{pos}-{mq}-{sam_flag}-{cigar}"


I note that you've added the (very lengthy!) cigar string to the signature. I believe that the signature is already unnecessarily lengthy.
It seems to me that we could filter out secondary and supplementary alignments, and then just drop duplicate queryNames. After all, according to the SAM spec:
"For each read/contig in a SAM file, it is required that one and only one line associated with the
read satisfies ‘FLAG & 0x900 == 0’. This line is called the primary line of the read."
You could still reset your dictionaries at each new position, assuming that the duplicates you're trying to eliminate always occur at the same position.
Edit: Maybe you're adding all these fields to the signature so that they appear in the output list of duplicate signatures. If so, so be it. But you could still do the comparison on just the query name.

tedsharpe · 2023-10-17T19:02:35Z

docker/lr-basic/Dockerfile

@@ -64,6 +64,7 @@ RUN apt-get -qqy update --fix-missing && \
                 gnupg \
                 less \
                 tabix \
+                 tree \


Might be good to add bsdmainutils, too. There are a number of WDLs that use the "column" command.

tedsharpe · 2023-10-19T15:23:09Z

wdl/tasks/Utility/BAMutils.wdl

+        num_rgs=$(wc -l one_rg_per_line.txt | awk '{pritn $1}')
+        if [[ ${num_rgs} -gt 1 ]]; then exit 1; fi
+
+        cat one_rg_per_line.txt | tr '\t' '\n' > rh_header.txt


useless cat. tr <one_rg_per_line.txt > rh_header.txt

tedsharpe · 2023-10-19T15:43:55Z

wdl/tasks/Utility/BAMutils.wdl

+        set -euxo pipefail
+
+        # assumption
+        sort_order=$(samtools view -H ~{qns_bam} | grep "^@HD" | tr '\t' '\n' | grep "^SO:" | awk -F ':' '{print $2}')


Consider sed rather than tr | grep | awk.

tedsharpe · 2023-10-19T18:18:23Z

wdl/tasks/Utility/BAMutils.wdl

+        # remote grab read names
+        echo "false" > samtools.failed.txt
+        samtools view ~{qns_bam} \
+        | awk -F '\t' '{print $1}' \


would've expected cut -f1 to be faster than awk, but it isn't.
shows that you always have to test, i guess.

tedsharpe · 2023-10-19T19:17:55Z

wdl/tasks/Utility/BAMutils.wdl

+
+        set +e
+        count=1
+        while true; do


I'm really surprised by two things:

That this is necessary. The token is good for an hour. We're just processing a shard of a BAM.

That it works. It appears to me that the token is only used once, at the time that the stream is opened.

Have you checked to see that re-auth is necessary and effective?
If so, I'd like to suggest that you put the re-auth while loop in the background rather than the main samtools command. You can then just abandon the background process when the main command exits: no testing counts or checking for exited processes.

tedsharpe · 2023-10-20T13:37:36Z

wdl/tasks/Utility/BAMutils.wdl

+
+    Int disk_size = 3 * ceil(size(qnorder_bam, "GB"))
+
+    command <<<


I'm wondering if you could vastly simplify this entire task with a simple awk program:

GCS_OAUTH_TOKEN=`gcloud auth application-default print-access-token` && export GCS_OAUTH_TOKEN touch duplicated.readnames.txt samtools view -h ~{qnorder_bam} | awk ' BEGIN{FS=OFS="\t"} {if(substr($0,1,1)=="@")print; else if($1==lastQName)print lastQName > "duplicated.readnames.txt"; else {print; lastQName = $1} }' | samtools view -bh - > "~{base}.dedup.bam"

tedsharpe · 2023-10-20T20:00:25Z

wdl/tasks/Utility/PBUtils.wdl

+        cat readgroup.movies.txt
+
+        # get primrose PG lines
+        grep "^@PG" header.txt | grep -v "^@SQ" | grep "^@PG" | grep -F 'ID:primrose' | tr '\t' '\n' | grep '^CL:' > primrose.pg.lines.txt


Forgot to note this one: the 2nd and 3rd greps look superfluous to me.

bshifaw

Some minor comments

bshifaw · 2023-10-25T14:52:54Z

wdl/pipelines/ONT/Preprocessing/RemoveDuplicateFromMergedONTBamAndSplitByReadgroup.wdl

+    input {
+        File bam
+    }
+    output {


Have the output block after the command but before runtime block

bshifaw · 2023-10-25T15:15:24Z

wdl/pipelines/ONT/Preprocessing/RemoveDuplicateFromMergedONTBamAndSplitByReadgroup.wdl

+        String gcs_out_root_dir
+    }
+
+    output {


Having the output block right after the input block helps quickly show what the workflow takes and what it gives. Just wanted to point out that miniwdl linting is throwing a warning because the tasks are being used before being called (WORKHORSE) in the workflow.

bshifaw · 2023-10-25T15:27:21Z

wdl/tasks/Utility/PBUtils.wdl

+        desciption: "Verify that a PacBio's BAM has primrose run on all its read groups"
+    }
+    input {
+        String bam


(Optional) You could set the bam variable as a File and avoid localization by using the localization_optional

bshifaw · 2023-10-25T15:35:16Z

wdl/pipelines/TechAgnostic/Utility/SplitBamByReadgroup.wdl

+        String platform_unit = GetReadGroupInfo.read_group_info['PU']
+
+        if (debug_mode) {
+            call Utils.CountBamRecords { input: bam = bam }


Similar comment made by Ted but for CountBamRecords calls

bshifaw · 2023-10-25T15:37:57Z

wdl/pipelines/TechAgnostic/Utility/SplitBamByReadgroup.wdl

+            call FF.FinalizeToFile as SaveFq {
+                input: file = BamToFastq.reads_fq, outdir = outdir
+            }
+            if (debug_mode) { call Utils.CountFastqRecords { input: fastq = BamToFastq.reads_fq } }


Same comment for CountFastqRecords, not being referenced

bshifaw · 2023-10-25T15:42:01Z

wdl/tasks/Utility/BAMutils.wdl

+        description: "Check several metadata of an input BAM (aliged? sort order? etc)"
+    }
+    input {
+        File bam


The localization_optional: true setting in the parameter_meta block of this task indicates to Cromwell that the file does not need to be downloaded

bshifaw · 2023-10-25T15:43:04Z

wdl/tasks/Utility/BAMutils.wdl

+        RuntimeAttr? runtime_attr_override
+    }
+
+    output {


output block after command block

bshifaw · 2023-10-25T15:43:31Z

wdl/tasks/Utility/BAMutils.wdl

+        RuntimeAttr? runtime_attr_override
+    }
+
+    output {


output after command block

bshifaw · 2023-10-25T15:44:11Z

wdl/tasks/Utility/BAMutils.wdl

+        Int? num_ssds
+        RuntimeAttr? runtime_attr_override
+    }
+    output {


output after command block

bshifaw · 2023-10-25T15:44:38Z

wdl/tasks/Utility/BAMutils.wdl

+        RuntimeAttr? runtime_attr_override
+    }
+
+    output {


output after command block

* update Hifiasm to version 0.19.5 * update how Hifiasm outputs are compressed (bgz replacing gz), also * monitor hifiasm resources usage

* update docker used in PBSV tasks to the version coming with official SMRTLink releases (2.9.0) * change how the 2-step PBSV process is done (following the recommended way now)

* to version 2.0.7 * using TRF bed * conditionally phase sv (requires phased bam) * generates its own vcf.gz and tbi

Overhaul how small variants are called in the WG pipelines * default to use DV to call small variants, Clair3 analysis needs to be requested explicitly * retire the Pepper toolchain completely from the CCS pipeline, using DV directly * for R10.4+ ONT data, also use DV directly * older ONT data would still use the PEPPER-DV-Margin pipeline * offers GPU version (though based on, it's not worth it yet) * update how bam haplotagging is done Cleanup structural variants calling * experiment with SNF2 phasing SV calls (implicitly depends on small variants calling now) * tune PBSV calling - discover now supports --hifi - output vcf.gz and tbi - less verbose logging by default Misc.: * optimizations to BAM merging and metrics workflow * updates coverage collection step * new R script to visualize log from vm_monitoring_script.sh

* organize dockstore.yml file a bit * make WDL validation shell script more usable * update pbmm2 and pbindex to versions in SMRTLink * update GeneralUtils.wdl - two bash-like new tasks [CoerceMapToArrayOfPairs, CoerceArrayOfPairsToMap] - cleanup task CollapseArrayOfStrings * update resource allocations to tasks - NanoplotFromBam (also changes docker) - MosDepthWGS

* incorporates gcloud cli (not just gsutil) * integrate libdeflate for more speedups

…sorted ONT bam

incorporate new tasks and optimize them * [CountMethylCallReads, GatherReadsWithoutMethylCalls] from sh_beans * [GetPileup, BamToRelevantPileup] from sh_more_atomic_qc * [GetReadGroupLines, GetSortOrder, SplitNameSortedUbam] from sh_ont_fc * [SamtoolsFlagStats, ParseFlagStatsJson] from sh_trvial_stats * [FilterBamByLen, InferSampleName] from sh_seqkit * [CountAlignmentRecords, StreamingBamErrored, CountAlignmentRecordsByFlag] from sh_maha_aln_metrics * [ResetSamplename] from sh_ingest_singlerg * [MergeBamsWithSamtools] from sh_ont_fc.Utils.wdl * [BamToFastq] from sh_more_bam_qcs and optimize it with sh_ingest_singlerg.Utils.wdl delete * GetSortOrder as that's now implemented in GatherBamMetadata * Drop2304Alignments as that's no longer used update dockers to the latest

CHERRY-PICK FROM VARIOUS QC/METRICS BRANCHES: * collect information about ML/MM tags in a long-read BAM (sh_beans) * a heuristic way to find peaks in a distribution (using dyst) (sh_dyst_peaker) * filter reads by length in a BAM * collect some read quality stats from (length-filtered) FASTQ/BAM (sh_seq_kit) * VerifyBamID2 (for contamination estimation) * naive sex-concordance check (sh_more_atomic_qc) * check fingerprint of a single BAM file (sh_sample_fp) * collect SAM flag stats (sh_trivial_stats)

* make BeanCounter finalization optional (wdl/pipelines/TechAgnostic/Utility/CountTheBeans.wdl) * custom struct for sub-workflow config using a JSON (wdl/pipelines/TechAgnostic/Utility/LongReadsContaminationEstimation.wdl) * make fingerprint checking subworkflow control size filtering (wdl/tasks/QC/FPCheckAoU.wdl) (wdl/pipelines/TechAgnostic/Utility/VerifyBamFingerprint.wdl) * fix a warning by IDE/miniwdl complaining WDL stdlib function length only applies to Array (wdl/tasks/Utility/BAMutils.wdl) * various updates to Finalize (wdl/tasks/Utility/Finalize.wdl) New tasks in (wdl/tasks/Utility/GeneralUtils.wdl) to * correctly convert Map to TSV * concatenate files

* SampleLevelAlignedMetrics.wdl * PBCLRWholeGenome.wdl

* new struct in AlignedBamQCandMetrics.wdl to facilicate as-sub-workflow calling * change parameters name for fingerprint workflows

* make saving of reads without methylation SAM tags optional * better parameter naming

(affects contamination estimation)

* make BAM split by readgroup more efficient * make the ONT duplicated reads verification step resilient to http transient errors * new task to verify primrose was run on PacBio BAM * new task to get basecall model from ONT BAM

* (for both PacBio, ONT) to split a BAM by read group, and reset the alignment * deduplicate and reset ONT aligned bam * optionally check primrose runs on PacBio input BAMs

…, samtools fastq, samtools reset]

SHuang-Broad requested review from kvg and bshifaw August 18, 2023 15:45

SHuang-Broad commented Aug 18, 2023

View reviewed changes

SHuang-Broad mentioned this pull request Aug 18, 2023

Fix bug in duplicate-removal code #406

Closed

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch 2 times, most recently from def61fc to 185abdb Compare September 1, 2023 14:09

SHuang-Broad requested a review from tedsharpe October 13, 2023 00:58

SHuang-Broad mentioned this pull request Oct 13, 2023

Follow up to PR #418 #430

Open

tedsharpe approved these changes Oct 20, 2023

View reviewed changes

tedsharpe reviewed Oct 20, 2023

View reviewed changes

bshifaw assigned SHuang-Broad Oct 25, 2023

bshifaw reviewed Oct 25, 2023

View reviewed changes

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from 12e954b to 185abdb Compare October 31, 2023 14:02

SHuang-Broad added 11 commits November 8, 2023 12:08

UPDATES TO THE Hifiasm pipeline:

e09dc65

* update Hifiasm to version 0.19.5 * update how Hifiasm outputs are compressed (bgz replacing gz), also * monitor hifiasm resources usage

For both CCS/ONT, update PBSV

d3afc4b

* update docker used in PBSV tasks to the version coming with official SMRTLink releases (2.9.0) * change how the 2-step PBSV process is done (following the recommended way now)

For both CCS/ONT, update Sniffles-2

3a2ac5c

* to version 2.0.7 * using TRF bed * conditionally phase sv (requires phased bam) * generates its own vcf.gz and tbi

CRITICAL: updating pbindex to that from the SMRTLink 12 release

29aa964

New docker that's intended to replace lr-basic:

1ff0912

* incorporates gcloud cli (not just gsutil) * integrate libdeflate for more speedups

New docker that offers a new mode remove duplicates from a queryname-…

4224ab4

…sorted ONT bam

New docker that updates the samtools in GATK to the latest (1.18)

a13a541

update some legacy code due to new interfaces in updated BAMutils.wdl

b16c619

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from d3c809c to 80e19c0 Compare December 4, 2023 03:45

SHuang-Broad added 4 commits December 19, 2023 17:07

New workflow to save files (avoid multiple calls to Finalize)

cbfed4e

New workflow to unify QC and metrics collection on aligned BAM

8039c88

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from a105f00 to f0367e9 Compare January 2, 2024 15:57

SHuang-Broad added 3 commits January 2, 2024 22:00

DEPRECATION:

b8e67ff

* SampleLevelAlignedMetrics.wdl * PBCLRWholeGenome.wdl

Update utility code:

e5a79c2

* new struct in AlignedBamQCandMetrics.wdl to facilicate as-sub-workflow calling * change parameters name for fingerprint workflows

a few tweaks to to AlignedBamQCandMetrics:

b5dc978

* make saving of reads without methylation SAM tags optional * better parameter naming

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch 3 times, most recently from c7a8e45 to e397755 Compare January 17, 2024 19:03

correct two mistakes: 1. a typo; 2. should append not replace

1abcc7a

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch 2 times, most recently from 28a2bde to bbbaa8c Compare January 19, 2024 17:11

make read-length util tasks pre-emptible

da5bfb0

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from bbbaa8c to e1322e8 Compare January 22, 2024 02:56

SHuang-Broad added 4 commits January 23, 2024 00:20

make fingerprint reads extraction more efficient/resilient

7975d6d

fix a stupid logical typo error (in saving files)

5620d74

option to force run fingerprintcheck on small bams

e910bac

make pbindex task more frugal

04b82f8

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from 9337e79 to ffd7624 Compare February 2, 2024 17:20

Safer and more efficient way to do targetted pileup conversion

14be4e2

(affects contamination estimation)

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from ffd7624 to 2657c15 Compare March 1, 2024 22:20

SHuang-Broad added 3 commits March 11, 2024 09:07

new task allowing one to send email in a WDL task

1f8dcae

lower resources for seqkit stats and read length gathering from fastq

1ac2a31

Chaning utility WDLs:

ce90a22

* make BAM split by readgroup more efficient * make the ONT duplicated reads verification step resilient to http transient errors * new task to verify primrose was run on PacBio BAM * new task to get basecall model from ONT BAM

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch 2 times, most recently from eb568d3 to aa0e08d Compare May 9, 2024 01:45

New workflows

3d3a508

* (for both PacBio, ONT) to split a BAM by read group, and reset the alignment * deduplicate and reset ONT aligned bam * optionally check primrose runs on PacBio input BAMs

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from aa0e08d to 86c6a57 Compare May 9, 2024 01:58

make several BAM ops more frugal: [GATK ValidateSamFile, GATK SortSam…

60517b4

…, samtools fastq, samtools reset]

SHuang-Broad force-pushed the sh_ingest_high_pass_inputs branch from 86c6a57 to 60517b4 Compare May 9, 2024 03:01

SHuang-Broad added 2 commits May 12, 2024 19:05

a typo bug fix

a3da0a6

seeing new format in ONT RG headers

6775bd3

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418

Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418

SHuang-Broad commented Aug 18, 2023

SHuang-Broad Aug 18, 2023

SHuang-Broad Aug 18, 2023

tedsharpe Oct 18, 2023

tedsharpe left a comment

tedsharpe Oct 17, 2023

tedsharpe Oct 17, 2023

tedsharpe Oct 17, 2023

tedsharpe Oct 23, 2023

tedsharpe Oct 17, 2023

tedsharpe Oct 17, 2023

tedsharpe Oct 19, 2023

tedsharpe Oct 19, 2023

tedsharpe Oct 19, 2023

tedsharpe Oct 19, 2023

tedsharpe Oct 20, 2023 •

edited

tedsharpe Oct 20, 2023

bshifaw left a comment

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

bshifaw Oct 25, 2023

@@ @@ -64,6 +64,7 @@ RUN apt-get -qqy update --fix-missing && \ @@
                                gnupg \
                                less \
                                tabix \
+                               tree \


		Int disk_size = 3 * ceil(size(qnorder_bam, "GB"))

		command <<<

Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418

Are you sure you want to change the base?

Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418

Conversation

SHuang-Broad commented Aug 18, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tedsharpe left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tedsharpe Oct 20, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

bshifaw left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tedsharpe Oct 20, 2023 •

edited