Mutation cohort (TCGA)

TL;DR: we turn a TCGA project's somatic mutation calls into one matrix, each column a distinct mutation and each row a tumor. JBrowse groups the rows by whichever clinical field you point it at, receptor status or stage, so the mutations a subtype shares line up.

Prerequisites

A JBrowse 2 instance to add tracks to (see the web quickstart, or the desktop quickstart, which loads these tracks by URL with nothing to host) and the JBrowse CLI
These files, hosted:

File	What
`https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz`	the cohort's somatic mutations
`https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv`	per-tumor histology, receptors, stage
`https://jbrowse.org/demos/tcga/tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz`	per-gene mutation rate per subtype

Where the data comes from

TCGA-BRCA open-access somatic mutation calls from the GDC (TCGA 2012).

primary-tumor Masked Somatic Mutation MAFs, queried and downloaded through the GDC API: https://api.gdc.cancer.gov/files
per-tumor clinical annotation, from harmonized case fields and each case's clinical XML: https://api.gdc.cancer.gov/cases
the cohort VCF, rehosted so the figures and their live links load without the GDC round trip: https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz
the recurrence track split by clinical group: https://jbrowse.org/demos/tcga/tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz
the clinical table those groups come from: https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv

The hg38 reference and gene track beside them are the hosted UCSC hub's own entries.

What the two files hold

The VCF is the GDC's per-tumor Masked Somatic Mutation calls merged into one multi-sample file, one column per tumor:

#CHROM POS       ID  REF ALT  INFO                          FORMAT    TCGA-A2-A0T2-01A  TCGA-A8-A07C-01A
chr3   179234297 .   A   G    GENE=PIK3CA;HGVSP=p.H1047R... GT:AD:DP  0/1:81,29:110     0/0

Two conventions bear on any figure read off that matrix:

0/0 marks a site the caller did not call, since a MAF carries no coverage record for one
every somatic call is written het, because a MAF gives no ploidy

Read counts are kept in AD/DP. INFO/CSQ re-encodes the MAF's own VEP columns (Consequence, IMPACT, HGVSp_Short, SIFT, PolyPhen), so the track can color cells by consequence impact without running an annotator.

The clinical TSV is one row per tumor barcode and one column per attribute:

name              histology  er        pr        her2      subtype    stage
TCGA-3C-AAAU-01A  lobular    positive  positive  negative  HR+/HER2-  X

Its columns come from three places:

histology and stage are the GDC's harmonized case fields
er/pr/her2 are read from each case's clinical XML, with in-situ hybridization taking precedence over immunohistochemistry for HER2
subtype is derived from those three; a tumor whose receptor calls do not resolve it stays unknown

The table names more tumors than either track carries, since a case with no mutation calls still has receptor status. JBrowse reports the unmatched tumors when the track loads.

Load the cohort VCF into JBrowse

The assembly comes first. The hosted FASTA calls its contigs bare (1) while the VCF uses chr1, so pass the alias file and both resolve.

export OUT=/var/www/html/jbrowse2

jbrowse add-assembly https://jbrowse.org/genomes/GRCh38/fasta/hg38.prefix.fa.gz \
  --name hg38 --type bgzipFasta \
  --refNameAliases https://s3.amazonaws.com/jbrowse.org/genomes/GRCh38/hg38_aliases.txt \
  --out $OUT

The cohort itself is a VariantTrack whose adapter carries the clinical TSV, with the matrix display on top of it:

Config file CLI (add-track-json) JBrowse Desktop

{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 1010,
      "featureColor": "jexl:impactColor(feature)"
    }
  ]
}

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 1010,
      "featureColor": "jexl:impactColor(feature)"
    }
  ]
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 1010,
      "featureColor": "jexl:impactColor(feature)"
    }
  ]
}

Three settings there:

The matrix display lays columns out by feature index, so a gene's mutations pack together however far apart they sit, with a connector band from each column to its position
featureColor with the impactColor helper colors each cell by its VEP impact tier from CSQ, the same as Color by... → Consequence impact in the track menu
samplesTsvLocation makes the clinical columns available to group and color rows by

Rows auto-fit by dividing height. A row goes below a pixel, so a band's mutation density reads as how dark it is.

Group the rows by clinical annotation

groupBy names a column of the samples TSV and makes each of its values a contiguous band of rows. colorBy puts the matching color strip in the gutter. colorBy is also Color by... → Samples in the track menu; groupBy is config only.

Config file CLI (add-track-json) JBrowse Desktop

{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 450,
      "lineZoneHeight": 130,
      "featureColor": "jexl:impactColor(feature)",
      "groupBy": "histology",
      "colorBy": "histology"
    }
  ]
}

jbrowse add-track-json '{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 450,
      "lineZoneHeight": 130,
      "featureColor": "jexl:impactColor(feature)",
      "groupBy": "histology",
      "colorBy": "histology"
    }
  ]
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "VariantTrack",
  "trackId": "tcga_brca_mutations",
  "name": "TCGA-BRCA somatic mutations (979 primary tumors)",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "VcfTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutations.vcf.gz",
    "samplesTsvLocation": {
      "uri": "https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv"
    }
  },
  "displays": [
    {
      "type": "LinearMultiSampleVariantMatrixDisplay",
      "height": 450,
      "lineZoneHeight": 130,
      "featureColor": "jexl:impactColor(feature)",
      "groupBy": "histology",
      "colorBy": "histology"
    }
  ]
}

With both slots at histology, the rows band by how the tumor was called under the microscope:

CDH1's exons with rows grouped and colored by histology and the gene's introns collapsed. The truncating (HIGH impact) cells crowd into the lobular band and the much larger ductal band above it is nearly empty. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

E-cadherin loss is the defining lesion of lobular breast cancer (Ciriello et al. 2015).

To get that window, right-click CDH1 in the gene lane, choose Collapse introns, and Replace current view (see Gene track):

The whole CDH1 transcript reshaped to its exons from the gene's own context menu, and the 979-tumor matrix redrawn over the coding sequence. Open the session this video starts in ↗

Two more things in that figure travel to any gene-scale matrix:

lineZoneHeight (or the handle under the band) opens the connector band, which shows where in the transcript a gene's calls fall. CDH1's fan lands in exon after exon, as a tumor suppressor's truncating calls do
A ClinVar track puts the germline record beside the somatic one on the same coordinates. It reads at hotspot or single-exon zoom; across sixteen collapsed exons the lane is a barcode

The GDC's open mutation calls are exome only, so these figures are all gene-scale.

Group by receptor subtype

With the slots at subtype, the rows band by receptor status (TCGA 2012):

TP53's exons with rows grouped and colored by receptor subtype, introns collapsed. The triple-negative band is visibly the dense one and the much larger HR+/HER2- band above it is sparse, with the calls spread the length of the coding sequence. Open this view in JBrowse ↗ (large dataset — this session takes a while to load)

The bottom band is the tumors whose receptor calls do not resolve a subtype. Hovering a column names its mutation and consequence; clicking opens the variant popup with per-tumor read counts.

PIK3CA at the same zoom piles on two hotspots in the HR+/HER2- band, H1047R in the kinase domain and E542K/E545K in the helical one.

Add a mutation recurrence track

The bands above differ in height, so their darkness does not compare as a rate. mutation_recurrence.py, the counterpart of the copy-number cohort's cnv_recurrence.py, takes the same SAMPLES.tsv:COLUMN group spec and writes one interval per gene valued in percent of the cohort:

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/mutation_recurrence.py
python3 mutation_recurrence.py tcga_brca_mutations.vcf.gz by_subtype.bedGraph \
  --groups tcga_brca_clinical.tsv:subtype

#chrom  start      end        HR+/HER2-  HER2+  triple-negative  unknown
chr3    179199065  179234302  40.56      30.18  11.19            30.71
chr17   7670683    7676564    19.44      39.64  80.42            32.28

BedGraphTabixAdapter takes every column past end as its own signal, and a MultiQuantitativeTrack draws one row per group.

Config file CLI (add-track-json) JBrowse Desktop

{
  "type": "MultiQuantitativeTrack",
  "trackId": "tcga_brca_mutation_recurrence_by_subtype",
  "name": "TCGA-BRCA mutation recurrence by receptor subtype",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "BedGraphTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz"
  },
  "displayDefaults": {
    "height": 260,
    "minScore": 0,
    "maxScore": 100,
    "showRowSeparators": true
  }
}

jbrowse add-track-json '{
  "type": "MultiQuantitativeTrack",
  "trackId": "tcga_brca_mutation_recurrence_by_subtype",
  "name": "TCGA-BRCA mutation recurrence by receptor subtype",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "BedGraphTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz"
  },
  "displayDefaults": {
    "height": 260,
    "minScore": 0,
    "maxScore": 100,
    "showRowSeparators": true
  }
}'

In JBrowse Desktop, or in any running JBrowse Web session, open a view on this track’s assembly, then File → Open track..., choose Add track from pasted JSON, and paste:

{
  "type": "MultiQuantitativeTrack",
  "trackId": "tcga_brca_mutation_recurrence_by_subtype",
  "name": "TCGA-BRCA mutation recurrence by receptor subtype",
  "assemblyNames": ["hg38"],
  "category": ["TCGA"],
  "adapter": {
    "type": "BedGraphTabixAdapter",
    "uri": "https://jbrowse.org/demos/tcga/tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz"
  },
  "displayDefaults": {
    "height": 260,
    "minScore": 0,
    "maxScore": 100,
    "showRowSeparators": true
  }
}

minScore/maxScore pin every row to one axis. Open it above the matrix and each band has its own rate over it. TP53 climbs toward the triple-negative group where PIK3CA falls.

--impact sets what counts as a hit, defaulting to the HIGH and MODERATE tiers, the ones impactColor paints. The rate has no background model, and gene length enters directly: TTN ranks near the top on passenger mutations alone.

Cluster the rows by genotype

Clustering → Cluster rows by genotype... in the track menu orders the rows by genotype (see Clustering rows), gathering every carrier into one block so a hotspot column becomes a solid bar. It replaces the clinical bands while on.

Thin the matrix down to recurrent mutations

Filter by... → Minor allele frequency in the track menu (or the minorAlleleFrequencyFilter slot) drops low-frequency columns, which on somatic data keeps only recurrent mutations: at PIK3CA the hotspots survive and the private columns go. See filtering by allele frequency and missingness.

Each somatic call is one alt allele out of two, so a mutation carried by 10% of the cohort sits at 0.05. CDH1's truncating calls are spread along the gene, so a threshold that isolates a hotspot empties that window.

Use your own cohort

maf_to_vcf.py takes any directory of MAFs (vcf2maf output, cBioPortal study downloads, your own caller) whose rows carry Chromosome, Start_Position, the two allele columns, Tumor_Sample_Barcode, and CONTEXT. For grouping, any TSV whose first column matches the VCF's sample names works.

Where to go next

The copy-number cohort paints the same tumors one row each, and its next steps apply here too.

Reproduce it end to end

One script builds every file above for any project id: build_tcga_cohort_mutations.sh, which merges the MAFs with maf_to_vcf.py, assembles the clinical table with tcga_clinical_tsv.py, and tallies the per-gene rates with mutation_recurrence.py. It needs curl, python3, and bgzip + tabix from htslib, which on Debian/Ubuntu is apt install curl python3 tabix.

curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/build_tcga_cohort_mutations.sh
bash build_tcga_cohort_mutations.sh TCGA-BRCA 20 # 20 tumors, to test the pipeline
bash build_tcga_cohort_mutations.sh TCGA-BRCA    # the full cohort, ~10 minutes
npx --yes serve jbrowse2                         # then open the printed URL

It writes tcga_brca_mutations.vcf.gz (+ .tbi), tcga_brca_clinical.tsv and tcga_brca_mutation_recurrence_by_subtype.bedGraph.gz (+ .tbi), then a jbrowse2/ opening on PIK3CA with the recurrence rows over the matrix. The assembly is the hosted UCSC hg38 hub's entry copied in, so the reference is never downloaded.

Swap in any other project id (TCGA-LUAD, TCGA-COAD, ...) for a different cohort, with --no-receptors to tcga_clinical_tsv.py for a non-breast project. A third argument names the clinical column the recurrence track splits on; subtype is breast only, while histology and stage work for any project.

Four steps decide whether the track is correct:

Open-access files only. The Masked Somatic Mutation MAFs are the aliquot-merged ensemble calls with germline sites masked out, and need no dbGaP application
The sample type comes off the barcode inside each MAF. A GDC file query filters on what a case has, so asking for Primary Tumor keeps a metastasis MAF too. The merge step filters on the barcode's sample-type code (01, primary solid tumor), the same tumors the copy-number cohort paints
Indels are anchored off the MAF's CONTEXT column. A MAF writes a deletion against a - alt, where VCF needs a shared preceding base. That base is in CONTEXT, so no reference FASTA is fetched
One MAF per sample barcode. A few cases were sequenced twice under one barcode, and the merge keeps one aliquot. Sample names are truncated to the sample barcode, so a tumor is one row name in both tracks

References

Feedback on this tutorial is welcome: contact us.

Mutation cohort (TCGA)

Prerequisites#

Where the data comes from#

What the two files hold#

Load the cohort VCF into JBrowse#

Group the rows by clinical annotation#

Group by receptor subtype#

Add a mutation recurrence track#

Cluster the rows by genotype#

Thin the matrix down to recurrent mutations#

Use your own cohort#

Where to go next#

Reproduce it end to end#

See also#

References#

Prerequisites

Where the data comes from

What the two files hold

Load the cohort VCF into JBrowse

Group the rows by clinical annotation

Group by receptor subtype

Add a mutation recurrence track

Cluster the rows by genotype

Thin the matrix down to recurrent mutations

Use your own cohort

Where to go next

Reproduce it end to end

See also

References