Pipeline to fix a common problem with TCGA somatic VCFs on the GDC Portal: a single case can have multiple VCFs for the same caller, because different tumor aliquots (vials/portions) were sequenced at different times. The pipeline:
- Detects how many aliquots/cases have more than one VCF per caller.
- Filters the manifest down to a single VCF per (case, caller), picking the one generated from the tumor BAM with higher coverage.
- Downloads only the filtered files.
- Generates, per case, a consensus VCF with variants present in at least 2 callers, using MuTect2 as the reference VCF (output format).
A single conda environment with everything needed (gdc-client,
bcftools, tabix, python) — see Installation.
You also need the following, already downloaded and complete, from the GDC Data Portal (download cart):
manifest.txt("Manifest" button)metadata.json("Metadata" button)sample_sheet.tsv("Sample Sheet" button)token.txt("Download Token" button, personal credential)
The pipeline makes no calls to the GDC API: everything is computed from these local files, so they must be complete (an interrupted metadata/sample sheet download will silently produce partial results, with no warning beyond the line count).
./scripts/00_setup_environment.shCreates the tcga-vcf-consensus conda environment (defined in
environment.yml) with gdc-client, bcftools/tabix and python.
The rest of the scripts activate it automatically
(scripts/_activate_env.sh); you don't need to activate it by hand.
# 1. Download the full manifest (optional, if you also want the raw data)
./scripts/01_download_manifest.sh manifest.txt token.txt gdc_downloads
# 2. See how many aliquots have >1 VCF per caller (optional diagnostic)
python3 scripts/02_summarize_multi_vcf_per_aliquot.py metadata.json
# 3 + 4. Filter the manifest by coverage and download only the filtered set
./scripts/04_download_filtered.sh manifest.txt metadata.json sample_sheet.tsv token.txt gdc_downloads_filtered
# -> generates gdc_manifest_filtered.txt and case_caller_map.tsv
# 5. Per-case consensus (>=2 callers, MuTect2 as reference)
./scripts/06_consensus_calling.sh case_caller_map.tsv gdc_downloads_filtered consensus_vcfsTo run it on a SLURM cluster, copy
scripts/submit_consensus.sbatch.example to submit_consensus.sbatch
(at the repo root), adjust the paths if you're not using the default
names, and run sbatch submit_consensus.sbatch.
When a case has several VCFs for the same caller (different tumor
aliquots), the one generated from the tumor BAM with the highest
mean_coverage in metadata.json is selected. Each file's case_id is
taken from sample_sheet.tsv (the Case ID column).
For each case, the VCFs from all available callers are compared with
bcftools isec -n+2 -w1, always listing MuTect2 first: the result
is the set of MuTect2 records that also appear in at least one other
caller (support from 2+ callers in total). Cases with no MuTect2 VCF, or
with only one caller available, are skipped (no consensus possible).
Note: some callers (e.g. SvABA) ship their VCFs as plain gzip instead of
bgzip, which bcftools rejects. The consensus script always
re-compresses to bgzip before indexing, to avoid that failure.
environment.yml # conda env: gdc-client + bcftools/tabix + python
scripts/
├── 00_setup_environment.sh # creates/updates the conda environment
├── _activate_env.sh # internal helper, activates the env (meant to be sourced)
├── 01_download_manifest.sh # generic download via gdc-client
├── 02_summarize_multi_vcf_per_aliquot.py # diagnostic: aliquots with >1 VCF/caller
├── 03_filter_manifest_by_coverage.py # filters manifest: 1 VCF per case+caller
├── 04_download_filtered.sh # orchestrates filtering + download
├── 05_consensus_calling.py # consensus, >=2 callers, MuTect2 reference
├── 06_consensus_calling.sh # wrapper around the script above
└── submit_consensus.sbatch.example # SLURM submission template
Manifests, GDC tokens, exported metadata and VCFs are specific to each
download and are not versioned (see .gitignore). A GDC token is a
personal credential: it must never be pushed to a repository.