Running Your Own Whole-Genome Analysis in the Cloud

September 27, 2026 Oak

genomics wgs cloud variant-calling bioinformatics

A feathered reptile specimen with white vapor-like plumage on a black studio plinth, spine barbs glowing in sequence with a few amber.

By the end of this you will have, in your own object storage bucket: a CRAM aligned to GRCh38, a small-variant VCF with per-sample quality annotations, a structural variant VCF, a QC report you can read, and a tab-separated table of variants annotated with gene, consequence, gnomAD population frequency, and ClinVar assertion. Total compute cost for one 30x genome on spot instances runs roughly $8 to $25 depending on how much you parallelize and how many times you restart. You need: your raw data (FASTQ or an unaligned CRAM), an AWS or GCP account with a quota for at least 32 vCPUs of a memory-heavy instance family, Docker, and about 500 GB of scratch disk. Everything below assumes Linux. It also assumes you are doing this for your own understanding, not for clinical decisions: nothing produced here is a diagnostic result, and a variant that looks important should go to a clinician and a CLIA-certified lab for independent confirmation before anyone acts on it.

1. Get the raw data out of the vendor portal and into a bucket

Ask for FASTQ, not just a VCF. A VCF is someone else’s opinion about your reads, filtered with someone else’s thresholds, aligned to a reference they chose. If the vendor only offers an aligned BAM or CRAM, take it and keep it, but ask for the FASTQs anyway. For 30x paired-end 150 bp WGS expect roughly 100 GB total as gzipped FASTQ, or 15-20 GB as CRAM.

Most vendors hand you either a presigned S3/GCS URL or an Aspera/FTP endpoint. Do the transfer on a cloud VM in the same region as your bucket, not from your laptop. A t3.medium with a big EBS volume is fine, and it can stream directly:

# Direct S3-to-S3 if they give you a presigned bucket path
aws s3 cp --recursive s3://vendor-delivery/ABC123/ s3://my-genome/raw/ \
  --request-payer requester

# Otherwise stream over HTTP without staging to disk
curl -sL "$PRESIGNED_URL" | aws s3 cp - s3://my-genome/raw/sample_R1.fastq.gz

Verify checksums before you delete anything upstream. Vendors ship an md5sum.txt; check it, because a truncated gzip will not fail until the aligner is three hours in.

aws s3 cp s3://my-genome/raw/ . --recursive --exclude "*" --include "*.md5"
md5sum -c md5sum.txt
# And a structural check that costs nothing:
gzip -t sample_R1.fastq.gz && echo "gzip stream intact"

Set a lifecycle rule now, while you remember: raw FASTQs go to Glacier Deep Archive after 30 days. At about $0.001 per GB-month that is roughly $1.20 a year to keep 100 GB of raw reads forever, versus $28 a year in S3 Standard. Restore takes 12 hours, which is fine for data you will touch once a year.

2. Pick one reference and never silently change it

Use GRCh38 with the analysis set that includes decoy contigs and excludes the alt contigs from primary alignment unless your caller is alt-aware. The practical choice is the Broad’s Homo_sapiens_assembly38.fasta or the GRCh38_no_alt_analysis_set from NCBI. We use the no-alt analysis set: alt-aware alignment adds complexity and most downstream callers, including DeepVariant’s published models, were trained without it.

T2T-CHM13 is a better assembly. It closes centromeres, resolves segmental duplications, and fixes false variants that come from GRCh38’s collapsed repeats. The problem is the ecosystem: gnomAD frequencies, ClinVar coordinates, and most annotation caches are still GRCh38-first. Our recommendation is to run GRCh38 as the primary and, if you care about a specific hard region, run a second alignment to CHM13 and compare. Do not mix coordinates in one spreadsheet.

REF=s3://my-genome/ref/GRCh38_full_analysis_set_plus_decoy_hla.fa
aws s3 cp $REF . && aws s3 cp $REF.fai .

Record the exact reference MD5 in your run metadata. CRAM does not store the reference sequence, only a checksum pointer to it. If you lose the exact FASTA you used, your CRAM becomes very expensive to read.

3. Choose the machine before you choose the pipeline

For a single genome, a single large instance beats a distributed pipeline. Batch orchestration earns its keep at hundreds of samples. At n=1 it adds IAM debugging, job queue tuning, and a class of failure where a task dies and you find out forty minutes later.

The binding constraint is bwa-mem2’s index, which is memory-hungry: roughly 25-30 bytes per reference base, so plan on 60-80 GB resident just for the index on GRCh38. We use r6i.8xlarge (32 vCPU, 256 GB RAM) with a 1 TB gp3 volume provisioned to 500 MB/s and 12,000 IOPS. On spot this is typically $0.35-0.55 per hour in us-east-1, versus about $2.02 on demand.

aws ec2 run-instances \
  --instance-type r6i.8xlarge \
  --instance-market-options 'MarketType=spot,SpotOptions={SpotInstanceType=one-time}' \
  --block-device-mappings '[{"DeviceName":"/dev/sda1","Ebs":{"VolumeSize":1000,"VolumeType":"gp3","Throughput":500,"Iops":12000}}]' \
  --iam-instance-profile Name=genome-runner \
  --image-id ami-xxxxxxxx --key-name mykey

If you want the alignment step faster and are willing to pay for it, the alternative is a DRAGEN-style FPGA instance (f1) or Illumina’s hosted DRAGEN, which does align-plus-call in well under an hour. We stay with open tools because the outputs are reproducible five years from now on commodity hardware.

Provisioned IOPS matters more than people expect. Sorting a 100 GB BAM is an I/O-bound operation, and a default 3,000-IOPS gp3 volume will turn a 40-minute sort into three hours.

4. Align, mark duplicates, and write CRAM

Run everything in containers with pinned digests. Version drift in an aligner is not visible in the output file, and six months later you will not remember which minor release you used.

# Pull pinned images
docker pull quay.io/biocontainers/bwa-mem2:2.2.1--hd03093a_5
docker pull quay.io/biocontainers/samtools:1.20--h50ea8bc_0

# Build the index once, cache it in S3 (~85 GB of index files, ~1 hour)
bwa-mem2 index GRCh38_full_analysis_set_plus_decoy_hla.fa

# Align. -K fixes the batch size so output is deterministic across thread counts.
bwa-mem2 mem -t 30 -K 100000000 -Y \
  -R '@RG\tID:S1\tSM:SAMPLE1\tLB:lib1\tPL:ILLUMINA\tPU:HXXXX.1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -u -@ 8 -m 4G -T /scratch/sort - \
| samtools markdup -@ 8 --reference GRCh38_full_analysis_set_plus_decoy_hla.fa \
    --output-fmt CRAM --write-index -f markdup_stats.txt - sample.cram

Three things in that command are deliberate. -K 100000000 makes the alignment deterministic regardless of thread count, so a rerun produces a byte-identical result and you can diff runs. -R with a real read group is mandatory: GATK and DeepVariant both refuse to run without one, and fixing it afterward means rewriting the whole file. Piping through samtools avoids ever writing an intermediate BAM, which saves about 150 GB of disk writes and maybe 25 minutes.

Expect 3-5 hours wall clock on 30 vCPUs for a 30x genome. CRAM output should land around 15-18 GB. If it is 40 GB, your duplicate rate or your error profile is off, and you should look at step 7 before going further.

Skip BQSR. For a single modern Illumina sample with DeepVariant downstream, base quality score recalibration costs an hour and changes almost nothing, because DeepVariant learns the error profile from the pileup image rather than trusting the reported qualities. If you are running GATK HaplotypeCaller instead, keep BQSR.

5. Call small variants

We use DeepVariant for germline SNVs and indels on Illumina WGS. It outperforms HaplotypeCaller on indels in homopolymer regions and needs less parameter care. GATK remains the right choice if you want a joint-genotyped cohort, a gVCF that merges cleanly with other people’s gVCFs, or somatic calling.

docker run -v $PWD:/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --reads=/data/sample.cram \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/scratch/dv_tmp

On 32 vCPUs this runs in roughly 2-4 hours. make_examples is CPU-bound and shards well, call_variants is the TensorFlow inference step and will use every core you give it. Keep the intermediate directory on instance store or a fast local volume, not on a network mount.

Then structural variants, which the small-variant caller cannot see at all:

manta_config.py --bam sample.cram --referenceFasta GRCh38...fa \
  --runDir manta_run && manta_run/runWorkflow.py -j 30

docker run -v $PWD:/data brentp/smoove smoove call \
  --name sample --fasta /data/GRCh38...fa -p 8 --genotype /data/sample.cram

SV calls from short reads are noisy. Deletions above 300 bp are reasonably reliable, insertions are not, and anything inside a segmental duplication should be treated as a hypothesis, not a result. If SVs matter to you, long reads are the answer, not more short-read callers.

6. Wrap it in Nextflow once it works

Get the pipeline right by hand first. Then, if you are going to run it more than twice, port it. nf-core/sarek already implements exactly the above with sensible defaults, and the AWS Batch executor handles spot reclamation with retries.

// nextflow.config
process {
  executor = 'awsbatch'
  queue    = 'genomics-spot'
  errorStrategy = { task.exitStatus in [143,137,104,134,139] ? 'retry' : 'finish' }
  maxRetries = 3
  withName: 'BWAMEM2_MEM' { cpus = 32; memory = '120 GB' }
}
aws.batch.maxParallelTransfers = 8
aws.region = 'us-east-1'
workDir = 's3://my-genome/work'

The exit codes in errorStrategy matter: 143 and 137 are SIGTERM and SIGKILL, which is what a reclaimed spot instance looks like from inside the container. Retrying those and failing hard on everything else means a spot interruption costs you one task, not the run.

nextflow run nf-core/sarek -r 3.4.4 -profile docker \
  --input samplesheet.csv --genome GATK.GRCh38 \
  --tools deepvariant,manta,vep --wes false -resume

Use -resume religiously. Nextflow hashes task inputs and skips completed work, which turns a failed 8-hour run into a 40-minute recovery.

7. Check the QC before you look at a single variant

This is the step people skip and it is the one that changes conclusions. Run these and read the numbers.

# Coverage, insert size, error rate, duplicate rate
samtools stats -@ 8 --reference GRCh38...fa sample.cram > sample.stats
mosdepth --by 1000 --fast-mode -t 8 sample sample.cram

# Cross-sample contamination
VerifyBamID --SVDPrefix resource/1000g.phase3.100k.b38.vcf.gz.dat \
  --Reference GRCh38...fa --BamFile sample.cram

# Variant-level summary
bcftools stats sample.dv.vcf.gz > sample.vcfstats
multiqc .

The numbers to look at, with the thresholds we use:

Mean coverage at or above 30x, with at least 90% of the callable genome covered at 20x or more. Below that, heterozygous calls start dropping out and absence of a variant stops meaning anything.
VerifyBamID FREEMIX below 0.02. Above 0.03 and you have contamination or a sample mixture, and no amount of filtering fixes it.
Duplicate rate under about 10%. High duplication means a low-complexity library, and your effective coverage is lower than your nominal coverage.
Ti/Tv ratio around 2.0-2.1 genome-wide. A value near 1.5 means you are calling a lot of noise.
Total SNVs in the 3.5 to 4.5 million range against GRCh38, plus roughly 500,000-800,000 indels. Well outside that range and something upstream is wrong, or your ancestry is far from the reference in a way worth understanding before you filter.

If you want a hard accuracy number, run the same pipeline on the GIAB HG002 FASTQs and compare against the v4.2.1 truth set with hap.py or rtg vcfeval restricted to the high-confidence BED. A sane short-read pipeline lands above 99.5% F1 for SNVs and somewhere in the high 99s for indels inside confident regions. That gives you a per-pipeline error rate to quote to yourself when a call looks surprising.

8. Annotate, then filter to something a human can read

Run VEP with a local cache, not the REST API. The API is rate-limited and a 4-million-variant VCF will take days.

docker run -v $PWD:/data -v $HOME/.vep:/vep ensemblorg/ensembl-vep:release_112.0 \
  vep --cache --offline --dir_cache /vep --assembly GRCh38 \
      --fasta /data/GRCh38...fa --species homo_sapiens \
      --input_file /data/sample.dv.vcf.gz --format vcf \
      --output_file /data/sample.vep.vcf.gz --vcf --compress_output bgzip \
      --everything --pick_allele_gene --fork 16 \
      --custom /data/gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_grpmax \
      --custom /data/clinvar_20240819.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT

Then reduce to a table:

bcftools +split-vep sample.vep.vcf.gz -f \
  '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%SYMBOL\t%Consequence\t%IMPACT\t%gnomADg_AF\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT[\t%GT\t%GQ\t%DP]\n' \
  -d -A tab -i 'FILTER="PASS"' > sample.annotated.tsv

From 4 million variants, filtering to PASS, IMPACT=HIGH or a ClinVar pathogenic assertion with at least two-star review status, and gnomAD allele frequency below 1% typically leaves a few hundred rows. Most are in genes with no relevance to you, or are recessive carrier status for conditions you already know are not present. ClinVar assertions with zero or one star are frequently wrong and get reclassified: check CLNREVSTAT before you believe anything.

Two categories need dedicated tools and will be wrong if you read them straight out of the VCF. Pharmacogenetic star alleles, especially CYP2D6, require copy-number-aware callers like Aldy, Cyrus, or StellarPGx, because CYP2D6 has hybrid genes and deletions that short-read diploid calling misrepresents. HLA typing needs a dedicated typer such as HLA-LA or OptiType. Anything here that looks medically significant is a starting point for a conversation with a genetic counselor and a confirmatory clinical test, not a conclusion.

Common problems

The run dies at 90% and you lose everything. Spot reclamation with no checkpointing. Either use Nextflow with an S3 work directory and -resume, or checkpoint by hand: write the sorted CRAM to S3 before starting variant calling, so a lost instance costs you one stage.

Disk fills during sorting. samtools sort writes temporary chunks equal to roughly the size of the uncompressed alignment. Point -T at a volume with 300 GB free and set -m conservatively: -m 4G with -@ 8 means 32 GB of RAM buffers, and if you set -m 16G with 32 threads you will get an OOM kill that looks like a mysterious exit 137.

CRAM will not open. CRAM decodes against the reference by MD5. If samtools view hangs or errors with a reference mismatch, you either lost the exact FASTA or REF_PATH is pointing at the EBI reference server and your instance has no outbound internet. Set REF_PATH=/local/cache/%2s/%2s/%s and REF_CACHE explicitly, or pass --reference on every command.

bwa-mem2 gets OOM-killed immediately. The index does not fit. Check with free -g that you have 80 GB or more available, or fall back to bwa mem, which is about 1.5-2x slower but needs roughly 6 GB.

Read group mismatch between FASTQ pairs. If the vendor split the run across lanes, align each lane separately with a distinct ID and PU but the same SM, then merge before marking duplicates. Marking duplicates across improperly merged lanes either misses real duplicates or flags legitimate reads.

Ti/Tv is fine but the variant count is 6 million. Usually means alignment against a reference with alt contigs using a non-alt-aware setup, which scatters reads and produces false heterozygous calls in paralogous regions. Realign against the no-alt analysis set.

Egress costs surprise you. Moving 100 GB out of AWS is about $9. Moving the same data between regions costs similarly. Keep compute, storage, and reference data in one region and pull down only the VCF and the QC report, which together are under 2 GB.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.