CNV cohort (TCGA)
TL;DR: tumors from different patients tend to gain and lose the same regions, because those regions carry a gene driving the cancer. We stack copy-number segment calls for 1104 TCGA breast tumors, one row per tumor colored by gain or loss, so a recurrent event reads as a vertical stripe down the stack.
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; the whole 1104-tumor cohort is a few MB of segment calls:
| File | What |
|---|---|
https://jbrowse.org/demos/tcga/tcga_brca_cnv.bed.gz | the segment stack |
https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence.bedGraph.gz | cohort gain/loss frequencies |
https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence_by_subtype.bedGraph.gz | the same, split by clinical group |
https://jbrowse.org/demos/tcga/tcga_brca_clinical.tsv | per-tumor histology, receptors |
Where the data comes from
TCGA-BRCA, from the GDC's open-access Masked Copy Number Segment files (Affymetrix SNP 6.0, harmonized to GRCh38), so no dbGaP application or token is needed.
- primary-tumor segment calls for 1104 tumors, 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 segment stack, rehosted so the figures and their live links load without the GDC round trip: https://jbrowse.org/demos/tcga/tcga_brca_cnv.bed.gz
- the cohort recurrence track and the same split by clinical group: https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence.bedGraph.gz and https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence_by_subtype.bedGraph.gz
- the clinical table the stack is grouped by: 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 files hold
Reproduce it end to end below builds these files from the GDC for any project id.
The BED is one segment call per line, with a #-prefixed header naming the
columns past end:
#chrom start end name sample segmean
chr1 3301764 30796057 +0.15 TCGA-3C-AAAU-01A 0.1480
chr1 3301764 7589655 -0.98 TCGA-3C-AALI-01A -0.9761sample is a TCGA barcode and is what splits the rows; segmean is the
caller's log2 tumor/normal ratio and is what colors them.
Every primary tumor in the project stacks as its own row, so a copy-number change shared across the cohort reads as a vertical stripe running down hundreds of samples.
Load the segments into JBrowse
Two commands set the whole thing up, an assembly and a track. Start with hg38,
where the one thing to watch is naming: the hosted FASTA calls its contigs bare
(1) while the BED uses chr1, so pass the alias file alongside it 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 $OUTThe track goes in as JSON rather than a second command line, since add-track's
flags do not reach inside a display:
{
"type": "FeatureTrack",
"trackId": "tcga_brca_cnv",
"name": "TCGA-BRCA copy number (1104 primary tumors)",
"assemblyNames": ["hg38"],
"category": ["TCGA"],
"adapter": {
"type": "BedTabixAdapter",
"uri": "https://jbrowse.org/demos/tcga/tcga_brca_cnv.bed.gz"
},
"displays": [
{
"type": "LinearMultiRowFeatureDisplay",
"partitionField": "sample",
"color": "jexl:feature.segmean<-1?'#2166ac':feature.segmean<-0.3?'#92c5de':feature.segmean<0.3?'#f7f7f7':feature.segmean<1?'#f4a582':'#b2182b'",
"legend": [
{ "label": "Deep loss (log2 < -1)", "color": "#2166ac" },
{ "label": "Loss", "color": "#92c5de" },
{ "label": "Balanced", "color": "#f7f7f7" },
{ "label": "Gain", "color": "#f4a582" },
{ "label": "Amplification (log2 > 1)", "color": "#b2182b" }
]
}
]
}The adapter's uri shorthand resolves the .tbi beside the file, and
rowHeight
auto-fits, which at this row count leaves every tumor a single pixel line. Three
settings do the rest:
partitionFieldsplits the one file into one labeled row persample. A thousand barcodes gives a thousand rows.coloris a jexl expression binningsegmeanonto a diverging blue-to-red scale, since this BED carries noitemRgb.legendnames each bin's log2 cutoff and color.
Cluster the stack
Open the track at whole-genome zoom, then run Clustering → Cluster rows by similarity from the track menu (see Clustering rows), which turns a noisy stack of 1104 tumors into blocks of shared copy-number profile.
A vertical stripe is one locus called the same way across many rows, blue for recurrent loss and red for recurrent gain. A whole row tending red or blue is one heavily aneuploid tumor, and clustering pulls those rows together into a band. Zooming to a single locus turns the stripe back into per-tumor calls, and clustering on that window alone sorts the cohort into its copy-number classes there.
The stack starts in barcode order, which encodes nothing, so the sorted state is the one every figure below is in.
At 1104 rows in a few hundred pixels each row is well under one pixel tall, so rows alias together and the saturated colors crowd out the neutral ones. The stack maps where the events are, and the track below counts how many rows carry them.
Add a recurrence track
The same calls collapsed to per-bin frequencies give that count its own axis.
Each 100 kb bin of tcga_brca_cnv_recurrence.bedGraph.gz carries the percent of
the cohort gained and the percent lost, on the same log2 cutoffs the stack
colors by (gain above 0.3, loss below -0.3), so a stripe and its peak count the
same tumors:
#chrom start end gain loss
chr1 204700000 204800000 58.88 -1.36
chr8 127600000 127800000 49.73 -0.91
chr16 89200000 89300000 3.26 -46.38BedGraphTabixAdapter reads every column past end as its own signal, so one
file carries both. Loss is written negative so that a wiggle's bicolorPivot at
0 draws gains up in posColor and losses down in negColor.
{
"type": "QuantitativeTrack",
"trackId": "tcga_brca_cnv_recurrence",
"name": "TCGA-BRCA recurrence (% of 1104 tumors)",
"assemblyNames": ["hg38"],
"category": ["TCGA"],
"adapter": {
"type": "BedGraphTabixAdapter",
"uri": "https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence.bedGraph.gz"
},
"displayDefaults": {
"height": 120,
"posColor": "#b2182b",
"negColor": "#2166ac",
"minScore": -100,
"maxScore": 100
}
}minScore and maxScore
(display options) pin
the axis to the whole cohort, so a bar means the same fraction wherever you
navigate. posColor/negColor reuse the stack's amplification and deep-loss
colors, so the two tracks agree by eye.
Placed above the stack, each peak sits over a stripe and puts a number on it:
Each bar is the fraction of the cohort carrying a call past the cutoff. There is no background model, no significance test and no peak calling, and amplitude enters only through the gain/loss cutoff; GISTIC is the tool for the significance test.
Split the recurrence by clinical group
The recurrence track above pools every tumor into one average.
cnv_recurrence.py --groups runs the same tally once per value of a clinical
column and writes each group its own gain and loss column:
curl -fO https://raw.githubusercontent.com/GMOD/jbrowse-components/main/scripts/cnv_recurrence.py
python3 cnv_recurrence.py tcga_brca_cnv.bed.gz by_subtype.bedGraph \
--groups tcga_brca_clinical.tsv:subtypeThe --groups file is the same
clinical TSV
the mutation cohort groups its matrix rows by, so a tumor falls in the same
group in both tracks.
The eight columns arrive as eight signals from that one file, and a
MultiQuantitativeTrack draws
one row each:
{
"type": "MultiQuantitativeTrack",
"trackId": "tcga_brca_cnv_recurrence_by_subtype",
"name": "TCGA-BRCA recurrence by receptor subtype",
"assemblyNames": ["hg38"],
"category": ["TCGA"],
"adapter": {
"type": "BedGraphTabixAdapter",
"uri": "https://jbrowse.org/demos/tcga/tcga_brca_cnv_recurrence_by_subtype.bedGraph.gz"
},
"displayDefaults": {
"height": 620,
"posColor": "#b2182b",
"negColor": "#2166ac",
"minScore": -70,
"maxScore": 70,
"showRowSeparators": true
}
}Row order follows the file's column order, so the four gain rows come first, in red, and the four loss rows under them, in blue. The bottom row of each block holds the tumors whose receptor calls do not resolve a subtype.
minScore/maxScore
pin the rows to one axis so they read against each other, at 70 here where the
pooled track above uses 100, since each row only ever fills half its axis.
Gain and loss stay separate columns, since at the edge of the 17q amplicon the HER2+ group is gained and lost at nearly the same rate.
--min-group sets how many tumors a subtype needs before it is plotted, since
at a handful of tumors a percentage moves in visible steps. The script names
each group it dropped and how big it was.
Point --groups at any other column for a different split; histology and
stage come from harmonized GDC fields and so work for any TCGA project, while
subtype is breast specific.
Use your own cohort
Nothing here is TCGA-specific. Any caller that emits per-sample segments works; the track config above only needs a BED with a sample column and a numeric column to color by:
#chrom start end name sample segmeanCNVkit .call.cns, ASCAT, and
PURPLE segments
all reshape into that shape with the same concatenate-and-tag step.
Where to go next
The same one-row-per-sample pattern carries to other GDC open-access data:
- Allele-specific copy number (ASCAT, open access at the GDC) is the closest
extension. It reports major and minor allele copy number separately, so it
shows loss-of-heterozygosity: a copy-neutral LOH region is balanced by total
copy number and has lost a parental allele. Same
.segshape, same display, only the coloring expression changes. - Methylation (Beta Value arrays, open access) is probe-level with genomic coordinates, and loads the same way with beta as the color field.
Splice junction quantification is not open access at the GDC, so a cohort splicing view needs controlled-access RNA-seq and a dbGaP application. For open splicing data, look outside TCGA (GTEx and recount3 publish junction summaries).
Reproduce it end to end
One script builds every file above for any project id:
build_tcga_cohort_cnv.sh,
which summarizes recurrence with
cnv_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_cnv.sh
bash build_tcga_cohort_cnv.sh TCGA-BRCA 20 # 20 tumors, to test the pipeline
bash build_tcga_cohort_cnv.sh TCGA-BRCA # the full cohort, ~20 minutes
npx --yes serve jbrowse2 # then open the printed URLIt writes tcga_brca_cnv.bed.gz (+ .tbi), the two recurrence bedGraphs (+
.tbi) and tcga_brca_clinical.tsv, then a jbrowse2/ opening on ERBB2 with
both recurrence lanes over one row per tumor. The assembly is the hosted UCSC
hg38 hub's own entry copied in, so the reference is never downloaded.
The full run is almost entirely downloading. Swap in any other project id
(TCGA-OV, TCGA-LUAD, ...) for a different cohort, and pass a third argument
to group the recurrence by a different clinical column, since subtype is
breast specific.
Three of its steps decide whether the resulting track loads correctly:
- Open-access files only. The GDC's Masked Copy Number Segment files
(Affymetrix SNP 6.0, already harmonized to GRCh38, germline CNV probes
removed) need no dbGaP application. The query filters to
Primary Tumor, leaving out each case's matched blood normal. - The
.segto BED reshape..segnames contigs bare (1), so the script adds thechrprefix, and.segstarts are 1-based inclusive against BED's 0-based half-open, so it subtracts 1. It keeps one file per barcode, so a tumor's row is one aliquot's calls. Segment_Meanis carried through unchanged.
The recurrence step is separately runnable as
cnv_recurrence.py,
if you have a cohort BED already and want only the frequency file. It skips bins
where fewer than half the cohort has any call, which here trims only the
chromosome tips. That coverage mask is taken over the whole cohort even when
--groups is set, so the grouped file has the same gaps as the pooled one.
The clinical table comes from
tcga_clinical_tsv.py,
the same helper and the same arguments the
mutation cohort uses, so one file
serves both tracks.
See also
- Multi-row feature track
- Quantitative track
- Mutation cohort (TCGA)
- QTL mapping (BXD mice)
- ChromHMM chromatin states
- CNV across a population (1000 Genomes)
- Structural variants (Cancer GIAB)
- Using jexl callbacks
References
Feedback on this tutorial is welcome: contact us.