DEV Community: W. Lee Pang, PhD

Tips for scalable workflows on AWS

W. Lee Pang, PhD — Fri, 01 Jul 2022 17:45:15 +0000

In May I had the opportunity to give a talk at the National Cancer Institute’s Containers and Workflows Interest Group (NCI-CWIG) on “Scalable and Reproducible Genomics Data Analysis on AWS”. I’ve been giving talks like this for a while and the themes haven’t changed much over the years. I try to refresh my slides every go round, and this time through I added the following one:

I spent about a minute or two talking through the above points during the session, but I think each warrants a bit more of a deep dive.

Keep your container images small

Containers are a means to package tools and their dependencies so that they can run anywhere. That being said, I’ve seen containers of many different sizes from miniscule (5MB) to gargantuan (>6GB).

On the large end, there are tooling distribution containers. They include many tools with the intention of being a single container to either run your entire workflow within with one compute step, or to be used for all steps of a workflow over distributed compute. I don’t recommend this for the following reasons.

ease managing tooling versions

Containers with multiple in-built tools make it difficult to manage individual tooling versions. Say you have a workflow that does the following:

tool-a (version n) → tool-b (version m) → tool-c (version p)

One day you need to update tool-b to version m+1 either to patch a critical bug or add a valuable new feature. Otherwise, you want to keep the versions of tool-a and tool-c the same. To manage this, you might try rebuilding the whole container, upgrading just tool-b. However, you might run into clashes if tool-b has shared dependencies with tool-a and tool-c. Alternatively, if you still need to run workflows that use the old version of tool-b, you now have two large container images to store and maintain.

Instead I recommend having one container image per tool. This helps to isolate concerns and reduces the amount of stuff that needs to be in each container. With each tool in their own container in the example above, updating tool-b to version m+1 is simply a matter of swapping out the old container image for tool-b with a new one, leaving tool-a and tool-c untouched.

less network bandwidth when pulling containers

Big containers will use a lot of network bandwidth to pull down. This can be problematic if (a) the compute instance a job is placed on doesn’t have much network throughput, (b) there are multiple copies of the job on the same compute instance pulling the container image, (c) the upstream image registry has limits on how much container data you can pull.

Bioinformatics workflows commonly include scatter-gather steps to leverage parallelism - e.g. processing many small intervals of a larger sequence at once using many small compute jobs. Scatters can be large - according to latest research, there are ~64,000 genes in the human genome (roughly 20,000 that code proteins) and one can imagine each being an interval to process. Now imagine the network and throttling issues with pulling a 6GB container 20,000 to 64,000 times. There are ways to partially mitigate this, like enabling container image caching on instances, using Amazon ECR pull through caching, and using instance types with more network bandwidth for compute jobs. Of course, all of these strategies are made even more effective if the job container is small to begin with.

... but avoid caveats of going too small

There are a lot of tools built with C/C++ using glibc shared libraries. The AWS CLI v2 is one of these tools. It is common for workflow engines running on AWS to bind mount the AWS CLI from the host instance into the container so that it is available for interacting with other AWS services like staging data from Amazon S3. Challenges arise when a tooling container is based on an image without glibc shared libraries as is the case with ultra-minimal base images like alpine and busybox. You can still use these ultra-minimal images, but you need to take extra steps to ensure that glibc shared libraries are available. For example, the AWS CLI v2 is distributed with the shared libraries it needs, and to make it work on an alpine based container, you can modify the LD_LIBRARY_PATH environment variable in the container environment to point to where these shared libraries are installed.

Leverage embarrassing parallelism

When running workflows in the cloud you have lots of distributed compute resources available to you. While you can use one large instance to run all steps of a workflow, this is an anti-pattern that can result in higher than expected costs.

To run an entire workflow on a single instance, you’d need to size it for the biggest step in your workflow. In the case of genomics, this is likely to be at the start of a workflow - e.g. lot’s of memory to do a reference alignment - where the rest of the workflow has relatively modest cpu and memory requirements. While the alignment step gets the resources it needs, the instance will be under utilized for the remainder of the workflow. The non-resource heavy parts of the workflow will cost more to run than they have to.

Instead, break up the workflow - run each step on its own compute instance with resources right sized for each step. There are a lot of AWS instance types to choose from, and if you connect a workflow engine to AWS Batch, AWS Batch will manage picking the right one for you based on CPU and memory (and GPU) requirements.

Write infrastructure agnostic workflows

Portability and reproducibility are key issues in scientific research. There was a whole mini-symposium at BioIT World this year devoted to FAIR principles. What this boils down to is workflows should be able to run at multiple scales and on multiple compute infrastructure - e.g. your laptop, your local HPC, your favorite cloud provider, ... and repeat for your collaborators. This matches a workflow’s deployment life cycle: initial development and small scale testing on a laptop or local cluster, production scale execution in the cloud, and sharing with others once it gets published in a scientific journal.

One should avoid hard coding infrastructure specific details into a workflow definition. Workflow languages are meant to abstract the workflow execution logic from the compute resources it runs on. For example, most languages distill a workflow task to the following:

container image for a tool
compute requirements like cpu, memory
the command the container will run
the data sources the task needs to read
the data the task is supposed to produce

This is all fairly generic. That said, most languages also allow you to provide special per task options like:

what job queue to run on
how to handle retries
special user permissions or roles to access external resources

These are infrastructural specifics and limit workflow portability if hard coded in a workflow definition. Instead, one should parameterize these as much as possible and use execution “profiles” to adapt these values to the infrastructure used at workflow runtime. For example, for a Nextflow workflow, you can do this with a nextflow.config file that looks like the following:

profiles {
    stub {  # use for local development and rapid iteration
        process {
            memory = '1 MB'
            errorStrategy = 'terminate'
        }
    }

    hpc {
        process {
            executor = 'sge'
            queue = 'long'
            clusterOptions = '-pe smp 10 -l virtual_free=64G,h_rt=30:00:00'
        }
    }

    aws {
        process {
            executor = 'awsbatch'
            cliPath = '/opt/aws-cli/bin/aws'
            queue = 'my-spot-queue'
        }
    }
}

You can then invoke the workflow using one of the profiles using:

nextflow run ./myworkflow -profile aws

The profile configuration can also be easily modified to suit someone else’s environment without directly changing the workflow definition.

Use cloud native integration where possible

In addition to a wide variety of compute, running workflows on AWS also provides the ability to interact with other AWS services like Amazon RDS databases, key value stores like Amazon DynamoDB, purpose built AI/ML services like Amazon Comprehend Medical, and many more.

One common pattern to integrate with AWS from a workflow job is to call additional services using the AWS CLI. Overall, this works well, but there are a few considerations one should note when doing so. First and foremost, a workflow job needs to know where the AWS CLI installed and how to use it. You can do this by either installing the AWS CLI on the host compute and bind mounting it into the container job, or including the AWS CLI as part of the container image. That said, see my notes above on keeping container images small for associated caveats. Second, while the AWS CLI is great for scripting, for more complex operations direct integration via the AWS SDK is a better fit.

For example, a typical integration is accessing data in Amazon S3. There are lots of genomics data in Amazon S3. A common processing pattern is to stage input data from Amazon S3 onto a compute node, process it with a containerized tool, and stage the output data back to Amazon S3 ... wash and repeat for each step of a workflow. While it sounds counterintuitive, this approach is fairly economical and functions well overall. However, this approach can result in longer compute times since tools would need to wait for input data to stage in before processing can begin. This is more evident if a tool only needs to read a small portion of a larger data file (e.g. scatters over genomic intervals).

In contrast, processing can start immediately and only transfer what is necessary if tooling can read bytes of data directly from Amazon S3. Tools based on htslib can do this, so you can run something like:

samtools view -H s3://1000genomes/data/NA12878/alignment/NA12878.alt_bwamem_GRCh38DH.20150718.CEU.low_coverage.cra

Overall, this can help your workflows run faster and more efficiently (e.g. less use of node-attached storage).

Closing

There are a number of things you can do to ensure that your genomics workflows run at scale. Each of the points above are simple to implement on their own and can be combined as needed to best suit your workflow and tooling. I recommend giving each of these a try, and if you have tips of your own do share them in the comments!

Profiling workflows with the Amazon Genomics CLI

W. Lee Pang, PhD — Tue, 26 Apr 2022 16:15:46 +0000

Profiling is an essential part of developing genomics workflows. By identifying and eliminating expensive bottlenecks, users can make workflows performant and cost efficient. A common way to profile a workflow is to generate a timing chart of all tasks executed. The Amazon Genomics CLI (AGC) provides a unified experience for running workflows across multiple workflow engines. In doing so, it also enables a common way to generate timing charts and profile workflows. In this blog I’ll show you how this can be done with AGC and a little Python coding.

As a prerequisite, you’ll need to install AGC and activate your AWS account region for running workflows. This is covered in AGC’s documentation.

To run a workflow with AGC you use commands like the following - in this case I’m running the GATK data processing workflow that is included as an example when you install AGC:

cd ~/amazon-genomics-cli/examples/gatk-best-practices-project

agc context deploy onDemandCtx
agc workflow run -c onDemandCtx gatk4-data-processing

# returns a workflow run-id like:
# 20ceffe1-9664-4f41-a5e9-66ab5eaf57ac

Once the workflow completes, you can get task execution information from the workflow logs using:

agc logs workflow gatk4-data-processing -r 20ceffe1-9664-4f41-a5e9-66ab5eaf57ac

which looks like:

RunId: 20ceffe1-9664-4f41-a5e9-66ab5eaf57ac
State: COMPLETE
Tasks: 
    Name                        JobId                   StartTime               StopTime                ExitCode
    HaplotypeCallerGvcf_GATK4.MergeGVCFs        12d27ad2-df4e-44f7-90ba-6140bde14cf0    2022-03-24 17:50:06.575 +0000 UTC   2022-03-24 18:00:34.439 +0000 UTC   0
    HaplotypeCallerGvcf_GATK4.HaplotypeCaller   39e61999-079a-4671-9625-9a1fada62c6b    2022-03-24 17:29:12.955 +0000 UTC   2022-03-24 17:48:22.475 +0000 UTC   0
    HaplotypeCallerGvcf_GATK4.HaplotypeCaller   ce9a6f75-9c9d-4e55-b4f9-1f46c516a3b5    2022-03-24 17:29:12.961 +0000 UTC   2022-03-24 17:47:25.444 +0000 UTC   0
...

The output from agc logs workflow provides start and stop times for each task run by the workflow. We can use this to generate a simple Gantt chart of the workflow tasks. After dropping the above output into a file, you can parse out the tasks with the following Python code:

#!/bin/env python3

from datetime import datetime
import re

DT_FMTS = [
    '%Y-%m-%dT%H:%M:%SZ', '%Y-%m-%dT%H:%M:%S.%fZ', 
    '%Y-%m-%dT%H:%M:%S+00:00', '%Y-%m-%dT%H:%M:%S.%f+00:00',
    '%Y-%m-%d %H:%M:%S.%f %z %Z', '%Y-%m-%d %H:%M:%S %z %Z',
]

def parse_datetime(s, fmts = DT_FMTS):
    for fmt in fmts:
        try:
            dt = datetime.strptime(s, fmt)
            return dt
        except ValueError as e:
            pass
    else:
        raise ValueError(f"unable to parse datetime from '{s}'")


def parse_agc_log(log_file):
    with open(log_file, 'r') as f:
        lines = f.readlines()
        lines = list(filter(lambda x: x.strip(), lines))

    data = {
        k: lines[i].split(":")[1].strip()
        for i, k in enumerate(['runid', 'state'])
    }

    data['tasks'] = [
        dict(zip(
            ['name', 'job_id', 'start_time', 'end_time', 'exit_code'], 
            re.split('\t+', line.strip())
        )) 
        for line in lines[4:]
    ]

    for i, task in enumerate(data['tasks']):
        for k in ('start_time', 'end_time'):
            task[k] = parse_datetime(task[k])
        data['tasks'][i] = task

    data['tasks'] = sorted(data['tasks'], key=lambda x: (x['start_time'], x['name']))

    return data

Using the parse_agc_log() function above looks like:

parse_agc_log('logs/gatk4-data-processing__20ceffe1-9664-4f41-a5e9-66ab5eaf57ac.agc.log')

which will return something like:

{'runid': '20ceffe1-9664-4f41-a5e9-66ab5eaf57ac',
 'state': 'COMPLETE',
 'tasks': [{'name': 'HaplotypeCallerGvcf_GATK4.MergeGVCFs',
   'job_id': '26fcfea4-cfbc-4f04-95ef-115eb309f3df',
   'start_time': datetime.datetime(2022, 4, 11, 1, 55, 53, 734000, tzinfo=datetime.timezone(datetime.timedelta(0), 'UTC')),
   'end_time': datetime.datetime(2022, 4, 11, 2, 6, 40, 418000, tzinfo=datetime.timezone(datetime.timedelta(0), 'UTC')),
   'exit_code': '0'},
  {'name': 'HaplotypeCallerGvcf_GATK4.HaplotypeCaller',
   'job_id': 'a9b3e5f1-cbd8-482b-87ba-df34e0252571',
   'start_time': datetime.datetime(2022, 4, 11, 1, 36, 5, 426000, tzinfo=datetime.timezone(datetime.timedelta(0), 'UTC')),
   'end_time': datetime.datetime(2022, 4, 11, 1, 52, 42, 457000, tzinfo=datetime.timezone(datetime.timedelta(0), 'UTC')),
   'exit_code': '0'},
   ...

This works to parse workflow tasks for most cases. However, when the task name is long it can chomp the separator between the task name and the task job-id. For example:

# ok line
    PreProcessingForVariantDiscovery_GATK4.MergeBamAlignment    7e244b7b-e02d-4bb9-bee7-9005204382f3    2022-03-24 10:17:20.673 +0000 UTC   2022-03-24 10:23:00.67 +0000 UTC    0

# bad line
    PreProcessingForVariantDiscovery_GATK4.CreateSequenceGroupingTSV85c5a61b-867e-49ed-a11d-951910444c03    2022-03-24 10:05:53.583 +0000 UTC   2022-03-24 10:11:55.023 +0000 UTC   0

This makes it difficult for steps we want to do later that involve the job-id. The JSON formatter for agc is not currently used by the agc logs workflow command. I created an issue to address this to help the team prioritize it. I could probably solve this in the short term with some regular expressions, but I’d rather not. There is another way to consistently get machine parseable output.

Under the hood, agc uses the GA4GH Workflow Execution Schema (WES) to communicate with AWS resources to run workflows. Each AGC context deploys a GA4GH WES endpoint which you can send requests to directly. The endpoint requires AWS SigV4 signed requests, so the easiest way to talk to it is with a tool like awscurl, which does the request signing for you.

To get the WES endpoint for a specific context you can use the agc context describe command:

agc context describe onDemandCtx --format json

which returns something like:

{
    "Name": "onDemandCtx",
    "Status": "STARTED",
    "StatusReason": "",
    "MaxVCpus": 256,
    "RequestSpotInstances": false,
    "InstanceTypes": null,
    "Output": {
        "Url": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtx"
    },
    "WesEndpoint": {
        "Url": "https://yrtxeqnhkb.execute-api.us-east-1.amazonaws.com/prod/"
    }
}

Note the WES endpoint Url needs to have “ga4gh/wes/v1” appended to it. We can extract this value using jq.

wes=$(agc context describe onDemandCtx --format json | jq -r '.WesEndpoint.Url + "ga4gh/wes/v1"')

We can now use the WesEndpoint to send a GA4GH WES GetRunLog request for the workflow:

awscurl ${wes}/runs/20ceffe1-9664-4f41-a5e9-66ab5eaf57ac

which will return something like the following:

{
  "outputs": {
    "id": "20ceffe1-9664-4f41-a5e9-66ab5eaf57ac",
    "outputs": {
      "HaplotypeCallerGvcf_GATK4.output_vcf": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-MergeGVCFs/NA12878.g.vcf.gz",
      "HaplotypeCallerGvcf_GATK4.output_vcf_index": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-MergeGVCFs/NA12878.g.vcf.gz.tbi"
    }
  },
  "request": {
    "workflow_params": {},
    "workflow_type": "WDL",
    "workflow_type_version": "1.0"
  },
  "run_id": "0612738f-06cd-4d96-9438-a69c432dc59e",
  "state": "COMPLETE",
  "task_logs": [
    {
      "cmd": [
        "set -e\n\n/gatk/gatk --java-options \"-Xmx2G\"  \\\nMergeVcfs \\\n--INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-0/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-1/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-2/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-3/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-4/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-5/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-6/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-7/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-733263974272-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-8/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-9/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-10/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-11/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-12/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-13/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-14/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-15/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-16/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-17/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-18/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-19/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-20/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-21/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-22/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-23/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-733263974272-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-24/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-25/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-26/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-27/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-28/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-29/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-30/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-31/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-32/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-33/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-34/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-35/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-36/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-37/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-38/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-39/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-40/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-41/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-42/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-43/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-44/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-45/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-46/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-47/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-48/NA12878.g.vcf.gz --INPUT /cromwell_root/agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-49/NA12878.g.vcf.gz \\\n--OUTPUT NA12878.g.vcf.gz"
      ],
      "end_time": "2022-04-11T02:06:40.418Z",
      "exit_code": "0",
      "name": "HaplotypeCallerGvcf_GATK4.MergeGVCFs|26fcfea4-cfbc-4f04-95ef-115eb309f3df",
      "start_time": "2022-04-11T01:55:53.734Z",
      "stderr": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-MergeGVCFs/MergeGVCFs-stderr.log",
      "stdout": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-MergeGVCFs/MergeGVCFs-stdout.log"
    },
    {
      "cmd": [
        "set -e\n\n/gatk/gatk --java-options \"-Xmx6G -XX:GCTimeLimit=50 -XX:GCHeapFreeLimit=10\" \\\nHaplotypeCaller \\\n-R /cromwell_root/broad-references/hg38/v0/Homo_sapiens_assembly38.fasta \\\n-I s3://agc-733263974272-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-CramToBamTask/NA12878.bam \\\n-L /cromwell_root/gatk-test-data/intervals/hg38_wgs_scattered_calling_intervals/temp_0001_of_50/scattered.interval_list \\\n-O NA12878.g.vcf.gz \\\n-contamination 0 \\\n-G StandardAnnotation -G StandardHCAnnotation -G AS_StandardAnnotation \\\n-GQB 10 -GQB 20 -GQB 30 -GQB 40 -GQB 50 -GQB 60 -GQB 70 -GQB 80 -GQB 90 \\\n-ERC GVCF \\\n \\\n\n\n# Cromwell doesn't like optional task outputs, so we have to touch this file.\ntouch NA12878.g.vcf.gz.bamout.bam"
      ],
      "end_time": "2022-04-11T01:52:42.457Z",
      "exit_code": "0",
      "name": "HaplotypeCallerGvcf_GATK4.HaplotypeCaller|a9b3e5f1-cbd8-482b-87ba-df34e0252571",
      "start_time": "2022-04-11T01:36:05.426Z",
      "stderr": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-0/HaplotypeCaller-0-stderr.log",
      "stdout": "s3://agc-111122223333-us-east-1/project/GATKwdl/userid/user33IfmN/context/onDemandCtxCromwell1/cromwell-execution/HaplotypeCallerGvcf_GATK4/0612738f-06cd-4d96-9438-a69c432dc59e/call-HaplotypeCaller/shard-0/HaplotypeCaller-0-stdout.log"
    },
    ...
  ],
  ...
}

If we capture this output to a file, we can process it with the following Python code:

#!/bin/env python3

import json

def parse_wes_runlog(log_file):
    with open(log_file, 'r') as f:
        runlog = json.load(f)
        tasks = runlog['task_logs']

    if runlog.get('run_log'):
        for prop in ('start_time', 'end_time'):
            runlog["run_log"][prop] = parse_datetime(runlog["run_log"][prop])

    for i, task in enumerate(tasks):
        for k in ('start_time', 'end_time'):
            tasks[i][k] = parse_datetime(task[k])

        name, job_id = task['name'].split("|")
        tasks[i]['job_name'] = name
        tasks[i]['job_id'] = job_id

    runlog['task_logs'] = sorted(tasks, key=lambda x: (x['start_time'], x['name']))

    return runlog

Now that we can get task timing information in a consistent manner, let’s do some plotting. For this, I’m going to use Bokeh which generates nice interactive plots.

#!/bin/env python3

from bokeh.models import ColumnDataSource
from bokeh.plotting import figure, output_notebook, show

TIME_SCALE_FACTORS = {
    "sec": 1, "min": 1/60, "hr": 1/3600
}

def plot_tasks_bokeh(tasks, time_units='min'):
    """plots task timing using engine reported timing (wes runlog)"""
    time_scale_factor = TIME_SCALE_FACTORS[time_units]

    tare = min([task['start_time'] for task in tasks])
    wall_time = (max([task['end_time'] for task in tasks]) - tare).total_seconds() * time_scale_factor

    stats = {
        'num_tasks': len(tasks),
        'timing': {
            'wallclock': wall_time,
            'units': time_units
        }
    }

    stats_text = (
        f"tasks: {stats['num_tasks']}, "
        f"wall_time: {stats['timing']['wallclock']:.2f} {stats['timing']['units']}"
    )

    source = ColumnDataSource(data=dict(
        y=range(len(tasks)),
        names=[task['job_name'] for task in tasks],
        running_left=[(task['start_time'] - tare).total_seconds() * time_scale_factor for task in tasks],
        running_right=[(task['end_time'] - tare).total_seconds() * time_scale_factor for task in tasks],
        text_x=[((task['end_time'] - tare).total_seconds() + 30) * time_scale_factor for task in tasks],
    ))

    output_notebook()

    p = figure(width=1280, height=800)
    p.hbar(y='y', left='running_left', right='running_right', height=0.8, source=source)
    p.text(x='text_x', y='y', text='names', alpha=0.4, text_baseline='middle', text_font_size='1.5ex', source=source)

    p.y_range.flipped = True
    p.xaxis.axis_label = f"task execution time ({time_units})"
    p.title.text = stats_text

    show(p)

    return p

Using this plotting function:

#!/bin/env python3
runlog = parse_wes_runlog()
p = plot_tasks_bokeh(runlog['task_logs'])

we get a plot that looks like this:

This gives us a great view of what order tasks are run in, how long each task takes, how many shards a parallel task creates, and how each task contributes to the overall workflow run time. This is a good start, but we can go deeper.

AGC schedules workflow tasks to AWS Batch. Looking at each task name, we see that AWS Batch job-ids are encoded, which we’ve already split out in parse_wes_runlog() as job_name and job_id. Using these job-ids, we can query AWS Batch for additional job details. To do that we’ll add the following Python functions to our script:

#!/bin/env python3

import boto3
from concurrent.futures import ThreadPoolExecutor, as_completed

def get_job_descriptions(batch, job_ids):
    jobs = []
    job_ids_sets = chunks(job_ids, 100)

    with ThreadPoolExecutor(max_workers=10) as executor:
        future_jobs = {
            executor.submit(batch.describe_jobs, jobs=job_ids_set): job_ids_set
            for job_ids_set in job_ids_sets
        }

        for future in as_completed(future_jobs):
            try:
                response = future.result()
                jobs += response["jobs"]
            except Exception as e:
                print(f"error retrieving job descriptions: {e}")

    descriptions = {job['jobId']:job for job in jobs}
    return descriptions


def chunks(l, n):
    """split list l into chunks of size n"""
    for i in range(0, len(l), n):
        yield l[i : i + n]

These functions help to gather AWS Batch job descriptions concurrently in groups of 100 to accommodate the limits of the AWS Batch GetJobDescriptions API call.

We can now update our parse_wes_runlog() function to include retrieving the corresponding AWS Batch job description for each task by adding:

def parse_wes_runlog(log_file, batch=boto3.client('batch')):
    # ...
    descriptions = get_job_descriptions(batch, [task['job_id'] for task in tasks])
    for i, task in enumerate(tasks):
        tasks[i]['job_description'] = descriptions.get(task['job_id'])
    # ...

This gives us the following additional metrics for each task:

The time interval between when the job was created in AWS Batch to when the job actually started running. This is the amount of time the job spent waiting in a job queue.
The completion status for each task - i.e. “SUCCEEDED” or “FAILED”.
The type and amount of compute resources the job requested, specifically number of vCPUs and GBs of memory

Note: AWS Batch only stores job information for 48hrs after it has reached a “completed” state (e.g. SUCCEEDED or FAILED). It is a good idea to get this information as soon as you can and cache them somewhere.

Combining task start and stop times with vCPU and memory details, you can get a rough estimate of how much each task costs. If you go to the AWS EC2 console, you can export a table that contains instance types and their Linux On-demand per hour cost. Doing this for non-(a, d, n) C, M, and R instances types with x86_64 architecture, and fitting a multilinear model to the data, I get the following formula which I’ve encapsulated in a Python function:

#!/bin/env python3

def approx_instance_cost_per_hr(vcpu, mem_gib):
    # rough OLS fit model based on c, m, r instance pricing filtered for x86_64 and non-(a,d,n)
    # returns linux on-demand per hour cost in USD
    usd = 0.03550716*vcpu + 0.003633091*mem_gib + 0.00718333

    return usd

Note: AWS prices vary by region and are updated regularly, so your model constants may differ.

We can now update our plotting function in the following ways. First we’ll add more data to our ColumnDataSource:

    millis_to_sec = 1/1000
    mebi_to_gebi = 1/1024

    jobs = [task['job_description'] for task in tasks]
    tare = min([job['createdAt'] for job in jobs])

    colors={'SUCCEEDED': 'cornflowerblue', 'FAILED': 'orangered'}

    source = ColumnDataSource(data=dict(
        y = range(len(jobs)),
        names=[job['jobName'] for job in jobs],
        color=[colors[job['status']] for job in jobs],
        cost=[
            approx_instance_cost_per_hr(
                get_job_resources(job)['vcpus'], 
                get_job_resources(job)['memory'] * mebi_to_gebi
            ) * ((job['stoppedAt'] - job['startedAt']) * millis_to_sec / 3600)
            for job in jobs],
        cpus=[get_job_resources(job)['vcpus'] for job in jobs],
        memory=[get_job_resources(job)['memory'] * mebi_to_gebi for job in jobs],
        queued_left=[(job['createdAt'] - tare) * millis_to_sec * time_scale_factor for job in jobs],
        queued_right=[(job['startedAt'] - tare) * millis_to_sec * time_scale_factor for job in jobs],
        running_left=[(job['startedAt'] - tare) * millis_to_sec * time_scale_factor for job in jobs],
        running_right=[(job['stoppedAt'] - tare) * millis_to_sec * time_scale_factor for job in jobs],
        text_x=[((job['stoppedAt'] - tare) * millis_to_sec + 30) * time_scale_factor for job in jobs],
    ))

A few things to note. AWS Batch reports times in millisecond timestamps and memory in Mebibytes, hence the need for the millis_to_sec and mebi_to_gebi conversion factors. I’ve added data columns to associate with task completion state (color) as well as data to show task time spent in queue (queued_left and queued_right).

Our main task timing plot is the same as before with the addition of an extra set of bars to show how long tasks were in queue:

p.hbar(y='y', left='queued_left', right='queued_right', height=0.8, color='lightgrey', source=source, legend_label="queued")

We can also add linked plots for per task cpu requested:

p_cpu = figure(width=160, y_range=p.y_range, sizing_mode="stretch_height")
p_cpu.hbar(y='y', right='cpus', height=0.8, color="darkgrey", source=source)
p_cpu.x_range = Range1d(-1, 32)
p_cpu.xaxis.axis_label = "vcpus"
p_cpu.title.text = f"max cpus: {max(source.data['cpus'])}"
p_cpu.yaxis.visible = False

per task memory requested:

p_mem = figure(width=160, y_range=p.y_range, sizing_mode="stretch_height")
p_mem.hbar(y='y', right='memory', height=0.8, color="darkgrey", source=source)
p_mem.x_range = Range1d(-1, 32)
p_mem.xaxis.axis_label = "memory (GiB)"
p_mem.title.text = f"max mem: {max(source.data['memory']):.2f} GiB"
p_mem.yaxis.visible = False

and per task cost estimates:

p_usd = figure(width=160, y_range=p.y_range, sizing_mode="stretch_height")
p_usd.hbar(y='y', right='cost', height=0.8, color='limegreen', source=source)
x_max = 0.1
x_min = -(x_max * 0.05)
p_usd.x_range = Range1d(x_min, x_max*1.05)
p_usd.xaxis.axis_label = "task cost (USD)"
p_usd.title.text = f"total: {total_cost:.2f} USD"

We can also compute some overall workflow summary statistics like total_cpu_time and importantly total_cost:

cpus = [get_job_resources(job)['vcpus'] for job in jobs]
total_cpu_time = sum([(job['stoppedAt'] - job['startedAt']) * cpu for job, cpu in zip(jobs, cpus)]) * millis_to_sec * time_scale_factor
total_cost = sum(source.data['cost'])

Pulling the above together using a Bokeh gridplot layout, we get a plot that looks like the following:

Note: The input for the above workflow is a small dataset used for testing and demo purposes. The overall cost shown is not representative of a workflow processing a full 30x human whole genome sample.

With the plot above, you can see at a glance which task(s) impact the workflow the most in terms of execution time and cost. You can also see where tasks might be spending too much time in queue waiting for available resources. In both cases, you now have some insight into how the workflow can be optimized - e.g by tuning individual task resource requirements or adjusting overall compute availability via the Max vCPUs allowed in AWS Batch Compute Environments via AGC’s context configuration.

Importantly, the method above can be used for any workflow that AGC can run so you don’t have to remember engine specific details of how to get such profiling information. This opens up the possibility of comparing the performance a workflow across workflow engines, which is something I’ll explore further in a future blog post. In the meantime, I recommend generating some plots of your own! All the Python code I used above is available in Github.

Introducing myself

W. Lee Pang, PhD — Tue, 26 Apr 2022 16:14:36 +0000

Hi there! My name is Lee and I’m starting a new role as a Principal Developer Advocate at Amazon Web Services (AWS) on our HealthAI services team. I’ve had several roles at AWS like a Technical Business Developer, a Specialist Solutions Architect, and a Product Manager. All of them filled with interesting and fun challenges to solve on behalf of AWS customers.

Like several of my developer advocate colleagues, I want to help developers make the best use of the cloud. In particular, I want to do this for the bioinformatics and genomics community. There is so much genomics data being generated in the world - analyzing it at scale in the cloud is how we’re going to make sense of it all.

Working at AWS has been an exciting journey filled with an incredible amount of learning, and I don’t see that stopping any time soon. As a developer advocate, I’m expecting to learn a lot from the bioinformatics and genomics community about what works well and how things can be improved. Everything I learn, I plan to share here.

I’ve got lots of ideas and projects I’m working on, from how to use new AWS products and services for genomics, to highlighting cool bioinformatics tools that work well in the cloud. First in the queue is how to profile workflows using the Amazon Genomics CLI. I hope you enjoy it, and am looking forward having you dive in with me for what comes next!