diff --git a/README.md b/README.md
index ddc640ef..a5171968 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
# Output directory containing the formatted manuscript
The [`gh-pages`](https://github.com/AlexsLemonade/OpenPBTA-manuscript/tree/gh-pages) branch hosts the contents of this directory at
This manuscript -(permalink) +(permalink) was automatically generated -from AlexsLemonade/OpenPBTA-manuscript@cfc4c7d +from AlexsLemonade/OpenPBTA-manuscript@b2753bb on April 6, 2023.
All original code was developed within the following repositories and is publicly available as follows. Primary data analyses can be found at https://github.com/d3b-center/OpenPBTA-workflows. @@ -2230,19 +2231,19 @@
We used BWA-MEM
77 to align paired-end DNA-seq reads to the version 38 patch release 12 of the Homo sapiens genome reference, obtained as a FASTA file from UCSC (see Key Resources Table).
-Next, we used the Broad Institute’s Best Practices78 to process Binary Alignment/Map files (BAMs) in preparation for variant discovery.
-We marked duplicates using SAMBLASTER
79, and we merged and sorted BAMs using Sambamba
80
-We used the BaseRecalibrator
submodule of the Broad’s Genome Analysis Tool Kit GATK
81 to process BAM files.
-Lastly, for normal/germline input, we used the GATK
HaplotypeCaller
82 submodule on the recalibrated BAM to generate a genomic variant call format (GVCF) file.
+
We used BWA-MEM
78 to align paired-end DNA-seq reads to the version 38 patch release 12 of the Homo sapiens genome reference, obtained as a FASTA file from UCSC (see Key Resources Table).
+Next, we used the Broad Institute’s Best Practices79 to process Binary Alignment/Map files (BAMs) in preparation for variant discovery.
+We marked duplicates using SAMBLASTER
80, and we merged and sorted BAMs using Sambamba
81
+We used the BaseRecalibrator
submodule of the Broad’s Genome Analysis Tool Kit GATK
82 to process BAM files.
+Lastly, for normal/germline input, we used the GATK
HaplotypeCaller
83 submodule on the recalibrated BAM to generate a genomic variant call format (GVCF) file.
This file is used as the basis for germline calling, described in the SNV calling for B-allele Frequency (BAF) generation section.
We obtained references from the Broad Genome References on AWS bucket with a general description of references at https://s3.amazonaws.com/broad-references/broad-references-readme.html.
To confirm sample matches and remove mis-matched samples from the dataset, we performed NGSCheckMate
83 on matched tumor/normal CRAM files.
+
To confirm sample matches and remove mis-matched samples from the dataset, we performed NGSCheckMate
84 on matched tumor/normal CRAM files.
Briefly, we processed CRAMs using BCFtools
to filter and call 20k common single nucleotide polymorphisms (SNPs) using default parameters.
We used the resulting VCFs to run NGSCheckMate
.
Per NGSCheckMate
author recommendations, we used <= 0.61 as a correlation coefficient cutoff at sequencing depths > 10 to predict mis-matched samples.
-We determined RNA-Seq read strandedness by running the infer_experiment.py
script from RNA-SeQC
84 on the first 200k mapped reads.
+We determined RNA-Seq read strandedness by running the infer_experiment.py
script from RNA-SeQC
85 on the first 200k mapped reads.
We removed any samples whose calculated strandedness did not match strandedness information provided by the sequencing center.
We required that at least 60% of RNA-Seq reads mapped to the human reference for samples to be included in analysis.
During OpenPBTA analysis, we identified some samples which were mis-identified or potentially swapped.
@@ -2256,37 +2257,37 @@
For patients with hypermutant samples, we first added population frequency of germline variants using ANNOVAR
85 and pathogenicity scoring from ClinVar86 using SnpSift
87.
+
For patients with hypermutant samples, we first added population frequency of germline variants using ANNOVAR
86 and pathogenicity scoring from ClinVar87 using SnpSift
88.
We then filtered for variants with read depth >= 15, variant allele fraction >= 0.20, and which were observed at < 0.1% allele frequency across each population in the Genome Aggregation Database (see Key Resources Table).
Finally, we retained variants in genes included in the KEGG MMR gene set (see Key Resources Table), POLE, and/or TP53 which were ClinVar-annotated as pathogenic (P) or likely pathogenic (LP) with review status of >= 2 stars.
All P/LP variants were manually reviewed by an interdisciplinary team of scientists, clinicians, and genetic counselors.
This workflow is available in the D3b GitHub repository.
We used four variant callers to call SNVs and indels from paired tumor/normal samples with Targeted Panel, WXS, and/or WGS data: Strelka2
88, Mutect2
89, Lancet
90, and VarDictJava
91.
+
We used four variant callers to call SNVs and indels from paired tumor/normal samples with Targeted Panel, WXS, and/or WGS data: Strelka2
89, Mutect2
90, Lancet
91, and VarDictJava
92.
VarDictJava
-only calls were not retained since ~ 39M calls with low VAF were uniquely called and may be potential false positives. (~1.2M calls were called by Mutect2
, Strelka2
, and Lancet
and included consensus CNV calling as described below.)
We used only Strelka2
, Mutect2
and Lancet
to analyze WXS samples from TCGA.
TCGA samples were captured using various WXS target capture kits and we downloaded the BED files from the GDC portal
.
The manufacturers provided the input interval BED files for both panel and WXS data for PBTA samples.
We padded all panel and WXS BED files were by 100 bp on each side for Strelka2
, Mutect2
, and VarDictJava
runs and by 400 bp for the Lancet
run.
For WGS calling, we utilized the non-padded BROAD Institute interval calling list wgs_calling_regions.hg38.interval_list
, comprised of the full genome minus N bases, unless otherwise noted below.
-We ran Strelka2
88 using default parameters for canonical chromosomes (chr1-22, X,Y,M), as recommended by the authors, and we filtered the final Strelka2
VCF for PASS variants.
+We ran Strelka2
89 using default parameters for canonical chromosomes (chr1-22, X,Y,M), as recommended by the authors, and we filtered the final Strelka2
VCF for PASS variants.
We ran Mutect2
from GATK
according to Broad best practices outlined from their Workflow Description Language (WDL), and we filtered the final Mutect2
VCF for PASS variants.
-To manage memory issues, we ran VarDictJava
91 using 20 Kb interval chunks of the input BED, padded by 100 bp on each side, such that if an indel occurred in between intervals, it would be captured.
+To manage memory issues, we ran VarDictJava
92 using 20 Kb interval chunks of the input BED, padded by 100 bp on each side, such that if an indel occurred in between intervals, it would be captured.
Parameters and filtering followed BCBIO standards except that variants with a variant allele frequency (VAF) >= 0.05 (instead of >= 0.10) were retained.
The 0.05 VAF increased the true positive rate for indels and decreased the false positive rate for SNVs when using VarDictJava
in consensus calling.
We filtered the final VarDictJava
VCF for PASS variants with TYPE=StronglySomatic
.
We ran Lancet
using default parameters, except for those noted below.
For input intervals to Lancet
WGS, we created a reference BED from only the UTR, exome, and start/stop codon features of the GENCODE 31 reference, augmented as recommended with PASS variant calls from Strelka2
and Mutect2
.
We then padded these intervals by 300 bp on each side during Lancet
variant calling.
-Per recommendations for WGS samples, we augmented the Lancet input intervals described above with PASS variant calls from Strelka2
and Mutect2
as validation92.
Strelka2
and Mutect2
as validation93.
We normalized INDELs with bcftools norm
on all PASS VCFs using the kfdrc_annot_vcf_sub_wf.cwl
subworkflow, release v3 (See Table S5).
-The Ensembl Variant Effect Predictor (VEP
)93, reference release 93, was used to annotate variants and bcftools was used to add population allele frequency (AF) from gnomAD94.
-We annotated SNV and INDEL hotspots from v2 of Memorial Sloan Kettering Cancer Center’s (MSKCC) database (See Key Resources Table) as well as the TERT promoter mutations C228T and C250T95.
+The Ensembl Variant Effect Predictor (VEP
)94, reference release 93, was used to annotate variants and bcftools was used to add population allele frequency (AF) from gnomAD95.
+We annotated SNV and INDEL hotspots from v2 of Memorial Sloan Kettering Cancer Center’s (MSKCC) database (See Key Resources Table) as well as the TERT promoter mutations C228T and C250T96.
We annotated SNVs by matching amino acid position (Protein_position
column in MAF file) with SNVs in the MSKCC database, we matched splice sites to HGVSp_Short
values in the MSKCC database, and we matched INDELs based on amino acid present within the range of INDEL hotspots values in the MSKCC database.
We removed non-hotspot annotated variants with a normal depth less than or equal to 7 and/or gnomAD allele frequency (AF) greater than 0.001 as potential germline variants.
-We matched TERT promoter mutations using hg38 coordinates as indicated in ref.95: C228T occurs at 5:1295113 is annotated as existing variant s1242535815
, COSM1716563
, or COSM1716558
, and is 66 bp away from the TSS; C250T occurs at Chr5:1295135, is annotated as existing variant COSM1716559
, and is 88 bp away from the TSS.
+We matched TERT promoter mutations using hg38 coordinates as indicated in ref.96: C228T occurs at 5:1295113 is annotated as existing variant s1242535815
, COSM1716563
, or COSM1716558
, and is 66 bp away from the TSS; C250T occurs at Chr5:1295135, is annotated as existing variant COSM1716559
, and is 88 bp away from the TSS.
We retained variants annotated as PASS
or HotSpotAllele=1
in the final set, and we created MAFs using MSKCC’s vcf2maf
tool.
We retained all variant calls from Strelka2
, Mutect2
, or Lancet
that overlapped with an SNV or INDEL hotspot in a hotspot-specific MAF file, which we then used for select analyses as described below.
Strelka2
, Mutect2
, and Lancet
), and we did not further consider VarDictJava
calls due to concerns it called a large number of false positives.
This decision had minimal impact on results because VarDictJava
also identified nearly every mutation that the other three callers identified, in addition to many unique mutations.
We used Control-FREEC
96,97 and CNVkit
98 for copy number variant calls.
+
We used Control-FREEC
97,98 and CNVkit
99 for copy number variant calls.
For both algorithms, the germline_sex_estimate
(described below) was used as input for sample sex and germline variant calls (above) were used as input for BAF estimation.
Control-FREEC
was run on human genome reference hg38 using the optional parameters of a 0.05 coefficient of variation, ploidy choice of 2-4, and BAF adjustment for tumor-normal pairs.
-Theta2
99 used VarDictJava
germline and somatic calls, filtered on PASS and strongly somatic, to infer tumor purity.
+Theta2
100 used VarDictJava
germline and somatic calls, filtered on PASS and strongly somatic, to infer tumor purity.
Theta2
purity was added as an optional parameter to CNVkit
to adjust copy number calls.
CNVkit
was run on human genome reference hg38 using the optional parameters of Theta2 purity and BAF adjustment for tumor-normal pairs.
-We used GISTIC
100 on the CNVkit
and the consensus CNV segmentation files to generate gene-level copy number abundance (Log R Ratio) as well as chromosomal arm copy number alterations using the parameters specified in the (run-gistic
analysis module in the OpenPBTA Analysis repository).
GISTIC
101 on the CNVkit
and the consensus CNV segmentation files to generate gene-level copy number abundance (Log R Ratio) as well as chromosomal arm copy number alterations using the parameters specified in the (run-gistic
analysis module in the OpenPBTA Analysis repository).
For each caller and sample, we called CNVs based on consensus among Control-FREEC
96,97, CNVkit
98, and Manta
101.
+
For each caller and sample, we called CNVs based on consensus among Control-FREEC
97,98, CNVkit
99, and Manta
102.
We specifically included CNVs called significant by Control-FREEC
(p-value < 0.01) and Manta
calls that passed all filters in consensus calling.
-We removed sample and consensus caller files with more than 2,500 CNVs because we expected these to be noisy and derive poor quality samples based on cutoffs used in GISTIC
100.
+We removed sample and consensus caller files with more than 2,500 CNVs because we expected these to be noisy and derive poor quality samples based on cutoffs used in GISTIC
101.
For each sample, we included the regions in the final consensus set: 1) regions with reciprocal overlap of 50% or more between at least two of the callers; 2) smaller CNV regions in which more than 90% of regions are covered by another caller.
We did not include any copy number alteration called by a single algorithm in the consensus file.
We defined copy number as NA
for any regions that had a neutral call for the samples included in the consensus file.
We merged CNV regions within 10,000 bp of each other with the same direction of gain or loss into single region.
We filtered out any CNVs that overlapped 50% or more with immunoglobulin, telomeric, centromeric, segment duplicated regions, or that were shorter than 3000 bp.
We used Manta
101 for structural variant (SV) calls, and we limited to regions used in Strelka2
.
+
We used Manta
102 for structural variant (SV) calls, and we limited to regions used in Strelka2
.
The hg38 reference for SV calling used was limited to canonical chromosome regions.
-We used AnnotSV
102 to annotate Manta
output.
+We used AnnotSV
103 to annotate Manta
output.
All associated workflows are available in the workflows GitHub repository.
We used STAR
103 to align paired-end RNA-seq reads, and we used the associated alignment for all subsequent RNA analysis.
+
We used STAR
104 to align paired-end RNA-seq reads, and we used the associated alignment for all subsequent RNA analysis.
We used Ensembl GENCODE 27 “Comprehensive gene annotation” (see Key Resources Table) as a reference.
-We used RSEM
104 for both FPKM and TPM transcript- and gene-level quantification.
RSEM
105 for both FPKM and TPM transcript- and gene-level quantification.
To enable downstream analyses, we next identified gene symbols that map to multiple Ensembl gene identifiers (in GENCODE v27, 212 gene symbols map to 1866 Ensembl gene identifiers), known as multi-mapped gene symbols, and ensured unique mappings (collapse-rnaseq
analysis module in the OpenPBTA Analysis repository).
To this end, we first removed genes with no expression from the RSEM
abundance data by requiring an FPKM > 0 in at least 1 sample across the PBTA cohort.
@@ -2333,7 +2334,7 @@
We set up Arriba
105 and STAR-Fusion
106 fusion detection tools using CWL on CAVATICA.
+
We set up Arriba
106 and STAR-Fusion
107 fusion detection tools using CWL on CAVATICA.
For both of these tools, we used aligned BAM and chimeric SAM files from STAR
as inputs and GRCh38_gencode_v27
GTF for gene annotation.
We ran STAR-Fusion
with default parameters and annotated all fusion calls with the GRCh38_v27_CTAT_lib_Feb092018.plug-n-play.tar.gz
file from the STAR-Fusion
release.
For Arriba
, we used a blacklist file blacklist_hg38_GRCh38_2018-11-04.tsv.gz
from the Arriba
release to remove recurrent fusion artifacts and transcripts present in healthy tissue.
@@ -2345,19 +2346,19 @@
All p-values are two-sided unless otherwise stated. Z-scores were calculated using the formula \(z=(x –\mu)/\sigma\) where \(x\) is the value of interest, \(\mu\) is the mean, and \(\sigma\) is the standard deviation.
tumor-purity-exploration
module)Estimating tumor fraction from RNA directly is challenging because most assume tumor cells comprise all non-immune cells107, which is not a valid assumption for many diagnoses in the PBTA cohort. +
Estimating tumor fraction from RNA directly is challenging because most assume tumor cells comprise all non-immune cells108, which is not a valid assumption for many diagnoses in the PBTA cohort.
We therefore used Theta2 (as described in the “Somatic Copy Number Variant Calling section” Methods section) to infer tumor purity from WGS samples, further assuming that co-extracted RNA and DNA samples had the same tumor purity.
We then created a set of stranded RNA-Seq data thresholded by median tumor purity of the cancer group to rerun selected transcriptomic analyses: telomerase-activity-prediction
, tp53_nf1_score
, transcriptomic-dimension-reduction
, immune-deconv
, and gene-set-enrichment-analysis
.
Note that these thresholded analyses, which only considered stranded RNA samples that also had co-extracted DNA, were performed in their respective OpenPBTA analyses modules (not within tumor-purity-exploration
).
interaction-plots
analysis module)Using the consensus SNV calls, we identified genes that were recurrently mutated in the OpenPBTA cohort, including nonsynonymous mutations with a VAF > 5% among the set of independent samples.
-We used VEP
93 annotations, including “High” and “Moderate” consequence types as defined in the R package Maftools
108, to determine the set of nonsynonymous mutations.
+We used VEP
94 annotations, including “High” and “Moderate” consequence types as defined in the R package Maftools
109, to determine the set of nonsynonymous mutations.
For each gene, we then tallied the number of samples that had at least one nonsynonymous mutation.
For genes that contained nonsynonymous mutations in multiple samples, we calculated pairwise mutation co-occurrence scores. This score was defined as \(I(-\log_{10}(P))\) where \(I\) is 1 when the odds ratio is > 1 (indicating co-occurrence), and -1 when the odds ratio is < 1 (indicating mutual exclusivity), with \(P\) defined by Fisher’s Exact Test.
focal-cn-file-preparation
analysis module)We added the ploidy inferred via Control-FREEC
to the consensus CNV segmentation file and used the ploidy and copy number values to define gain and loss values broadly at the chromosome level.
-We used bedtools coverage
109 to add cytoband status using the UCSC cytoband file110 (See Key Resources Table).
+We used bedtools coverage
110 to add cytoband status using the UCSC cytoband file111 (See Key Resources Table).
The output status call fractions, which are values of the loss, gain, and callable fractions of each cytoband region, were used to define dominant status at the cytoband-level.
We calculated the weighted means of each status call fraction using band length.
We used the weighted means to define the dominant status at the chromosome arm-level.
IRanges
package in R111 to find overlaps between the segments in the consensus CNV file and the exons in the GENCODE v27 annotation file (See Key Resources Table) .
+To obtain the gene-level data, we used the IRanges
package in R112 to find overlaps between the segments in the consensus CNV file and the exons in the GENCODE v27 annotation file (See Key Resources Table) .
If the copy number value was 0, we set the status to “deep deletion”.
For autosomes only, we set the status to “amplification” when the copy number value was greater than two times the ploidy value.
-We plotted genome-wide gains and losses in (Figure S3C) using the R package ComplexHeatmap
112.
+We plotted genome-wide gains and losses in (Figure S3C) using the R package ComplexHeatmap
113.
chromosomal-instability
analysis module)We defined breakpoint density as the number of breaks per genome or exome per sample. For Manta SV calls, we filtered to retain “PASS” variants and used breakpoints from the algorithm. @@ -2378,20 +2379,20 @@
chromothripsis
analysis module)Considering only chromosomes 1-22 and X, we identified candidate chromothripsis regions in the set of independent tumor WGS samples with ShatterSeek113, using Manta SV calls that passed all filters and consensus CNV calls. +
Considering only chromosomes 1-22 and X, we identified candidate chromothripsis regions in the set of independent tumor WGS samples with ShatterSeek114, using Manta SV calls that passed all filters and consensus CNV calls.
We modified the consensus CNV data to fit ShatterSeek
input requirements as follows: we set CNV-neutral or excluded regions as the respective sample’s ploidy value from Control-FREEC
, and we then merged consecutive segments with the same copy number value.
We classified candidate chromothripsis regions as high- or low-confidence using the statistical criteria described by the ShatterSeek
authors.
immune-deconv
analysis module)We used the R package immunedeconv
pubmed:31510660? with the method quanTIseq
114 to deconvolute various immune cell types in tumors using collapsed FPKM RNA-seq, with samples batched by library type and then combined.
+
We used the R package immunedeconv
pubmed:31510660? with the method quanTIseq
115 to deconvolute various immune cell types in tumors using collapsed FPKM RNA-seq, with samples batched by library type and then combined.
The quanTIseq
deconvolution method directly estimates absolute fractions of 10 immune cell types that represent inferred proportions of the cell types in the mixture.
Therefore, we utilized quanTIseq
for inter-sample, intra-sample, and inter-histology score comparisons.
gene-set-enrichment-analysis
analysis module)We performed Gene Set Variation Analysis (GSVA) on collapsed, log2-transformed RSEM FPKM data for stranded RNA-Seq samples using the GSVA
Bioconductor package115.
-We specified the parameter mx.diff=TRUE
to obtain Gaussian-distributed scores for each of the MSigDB hallmark gene sets116.
+
We performed Gene Set Variation Analysis (GSVA) on collapsed, log2-transformed RSEM FPKM data for stranded RNA-Seq samples using the GSVA
Bioconductor package116.
+We specified the parameter mx.diff=TRUE
to obtain Gaussian-distributed scores for each of the MSigDB hallmark gene sets117.
We compared GSVA scores among histology groups using ANOVA and subsequent Tukey tests; p-values were Bonferroni-corrected for multiple hypothesis testing.
-We plotted scores by cancer group using the ComplexHeatmap
R package (Figure 5B)112.
ComplexHeatmap
R package (Figure 5B)113.
transcriptomic-dimension-reduction
analysis module)We applied Uniform Manifold Approximation and Projection (UMAP)117 to log2-transformed FPKM data for stranded RNA-Seq samples using the umap
R package (See Key Resources Table).
+
We applied Uniform Manifold Approximation and Projection (UMAP)118 to log2-transformed FPKM data for stranded RNA-Seq samples using the umap
R package (See Key Resources Table).
We considered all stranded RNA-Seq samples for this analysis, but we removed genes whose FPKM sum across samples was less than 100.
We set the UMAP number of neighbors parameter to 15.
fusion_filtering
analysis module)oncoprint-landscape
analysis module)We used Maftools
108 to generate oncoprints depicting the frequencies of canonical somatic gene mutations, CNVs, and fusions for the top 20 genes mutated across primary tumors within broad histologies of the OpenPBTA dataset.
-We collated canonical genes from the literature for low-grade gliomas (LGGs)24, embryonal tumors25,27,28,128,129, high-grade gliomas (HGGs)14,30,31,130, and other tumors: ependymomas, craniopharyngiomas, neuronal-glial mixed tumors, histiocytic tumors, chordoma, meningioma, and choroid plexus tumors32,131–137,pubmed:28623522?,pubmed:12466115?.
We used Maftools
109 to generate oncoprints depicting the frequencies of canonical somatic gene mutations, CNVs, and fusions for the top 20 genes mutated across primary tumors within broad histologies of the OpenPBTA dataset.
+We collated canonical genes from the literature for low-grade gliomas (LGGs)24, embryonal tumors25,27,28,129,130, high-grade gliomas (HGGs)14,30,31,131, and other tumors: ependymomas, craniopharyngiomas, neuronal-glial mixed tumors, histiocytic tumors, chordoma, meningioma, and choroid plexus tumors32,132–138,pubmed:28623522?,pubmed:12466115?.
mutational-signatures
analysis module)We obtained weights (i.e., exposures) for signature sets using the deconstructSigs
R package function whichSignatures()
138 from consensus SNVs with the BSgenome.Hsapiens.UCSC.hg38 annotations (see Key Resources Table).
+
We obtained weights (i.e., exposures) for signature sets using the deconstructSigs
R package function whichSignatures()
139 from consensus SNVs with the BSgenome.Hsapiens.UCSC.hg38 annotations (see Key Resources Table).
Specifically, we estimated signature weights across samples for eight signatures previously identified in the Signal reference set of signatures (“RefSig”) as associated with adult central nervous system (CNS) tumors40.
These eight RefSig signatures are 1, 3, 8, 11, 18, 19, N6, and MMR2.
Weights for signatures fall in the range zero to one inclusive.
deconstructSigs
estimates the weights for each signature across samples and allows for a proportion of unassigned weights referred to as “Other” in the text.
-These results do not include signatures with small contributions; deconstructSigs
drops signature weights that are less than 6%138.
-We plotted mutational signatures for patients with hypermutant tumors (Figure 4E) using the R package ComplexHeatmap
112.
deconstructSigs
drops signature weights that are less than 6%139.
+We plotted mutational signatures for patients with hypermutant tumors (Figure 4E) using the R package ComplexHeatmap
113.
snv-callers
analysis module)We consider tumor mutation burden (TMB) to be the number of consensus SNVs per effectively surveyed base of the genome. We considered base pairs to be effectively surveyed if they were in the intersection of the genomic ranges considered by the callers used to generate the consensus and where appropriate, regions of interest, such as coding sequences. @@ -2440,24 +2441,24 @@
Embryonal tumors were included in non-MB and non-ATRT embryonal tumor subtyping (molecular-subtyping-embryonal
analysis module) if they met any of the following criteria:
MN1::PATZ1
since it is an entity separate of CNS HGNET-MN1 tumors139.MN1::PATZ1
since it is an entity separate of CNS HGNET-MN1 tumors140.Non-MB and non-ATRT embryonal tumors identified with the above criteria were further subtyped (molecular-subtyping-embryonal
analysis module) using the criteria below140,141,pubmed:30249036?,pubmed:26389418?.
Non-MB and non-ATRT embryonal tumors identified with the above criteria were further subtyped (molecular-subtyping-embryonal
analysis module) using the criteria below141,142,pubmed:30249036?,pubmed:26389418?.
ETMR, C19MC-altered
(Embryonal tumor with multilayer rosettes, chromosome 19 miRNA cluster altered)122,142.ETMR, NOS
(Embryonal tumor with multilayer rosettes, not otherwise specified)122,142.ETMR, C19MC-altered
(Embryonal tumor with multilayer rosettes, chromosome 19 miRNA cluster altered)123,143.ETMR, NOS
(Embryonal tumor with multilayer rosettes, not otherwise specified)123,143.CNS HGNET-MN1
[Central nervous system (CNS) high-grade neuroepithelial tumor with MN1 alteration].CNS HGNET-BCOR
(CNS high-grade neuroepithelial tumor with BCOR alteration).CNS HGNET-BCOR
(CNS high-grade neuroepithelial tumor with BCOR alteration).CNS NB-FOXR2
(CNS neuroblastoma with FOXR2 activation).CNS EFT-CIC
(CNS Ewing sarcoma family tumor with CIC alteration)121CNS EFT-CIC
(CNS Ewing sarcoma family tumor with CIC alteration)122CNS Embryonal, NOS
(CNS Embryonal tumor, not otherwise specified).Neurocytoma subtypes central neurocytoma (CNC) and extraventricular neurocytoma (EVN) were assigned (molecular-subtyping-neurocytoma
analysis module) based on the primary site of the tumor144.
+
Neurocytoma subtypes central neurocytoma (CNC) and extraventricular neurocytoma (EVN) were assigned (molecular-subtyping-neurocytoma
analysis module) based on the primary site of the tumor145.
If the tumor’s primary site was “ventricles,” we assigned the subtype as CNC; otherwise, we assigned the subtype as EVN.
Craniopharyngiomas (CRANIO) were subtyped (molecular-subtyping-CRANIO
analysis module) into adamantinomatous (CRANIO, ADAM
), papillary (CRANIO, PAP
) or undetermined (CRANIO, To be classified
) based on the following criteria145,146:
Craniopharyngiomas (CRANIO) were subtyped (molecular-subtyping-CRANIO
analysis module) into adamantinomatous (CRANIO, ADAM
), papillary (CRANIO, PAP
) or undetermined (CRANIO, To be classified
) based on the following criteria146,147:
CRANIO, PAP
.CRANIO, ADAM
.If pathology diagnosis was Subependymal Giant Cell Astrocytoma (SEGA)
, the LGG
portion of molecular subtype was recoded to SEGA
.
Lastly, for all LGG- and GNT- subtyped samples, if the tumors were glialneuronal in origin, based on pathology_free_text_diagnosis
entries of desmoplastic infantile
,desmoplastic infantile ganglioglioma
, desmoplastic infantile astrocytoma
or glioneuronal
, each was recoded as follows:
If pathology diagnosis is Low-grade glioma/astrocytoma (WHO grade I/II)
or Ganglioglioma
, the LGG
portion of the molecular subtype was recoded to GNT
.
Ependymomas (EPN) were subtyped (molecular-subtyping-EPN
analysis module) into EPN, ST RELA
, EPN, ST YAP1
, EPN, PF A
and EPN, PF B
based on evidence for these molecular subgroups as described in Pajtler et al.131.
+
Ependymomas (EPN) were subtyped (molecular-subtyping-EPN
analysis module) into EPN, ST RELA
, EPN, ST YAP1
, EPN, PF A
and EPN, PF B
based on evidence for these molecular subgroups as described in Pajtler et al.132.
Briefly, fusion, CNV and gene expression data were used to subtype EPN as follows:
RELA
as fusion partner, e.g., C11orf95::RELA
, LTBP3::RELA
, was subtyped as EPN, ST RELA
.tp53_nf1_score
analysis module)We annotated TP53 altered HGG samples as either TP53 lost
or TP53 activated
and integrated this within the molecular subtype.
To this end, we applied a TP53 inactivation classifier originally trained on TCGA pan-cancer data42 to the matched RNA expression data, with samples batched by library type.
-Along with the TP53 classifier scores, we collectively used consensus SNV and CNV, SV, and reference databases that list TP53 hotspot mutations147,148 and functional domains149 to determine TP53 alteration status for each sample.
+Along with the TP53 classifier scores, we collectively used consensus SNV and CNV, SV, and reference databases that list TP53 hotspot mutations148,149 and functional domains150 to determine TP53 alteration status for each sample.
We adopted the following rules for calling either TP53 lost
or TP53 activated
:
If a sample had either of the two well-characterized TP53 gain-of-function mutations, p.R273C or p.R248W43, we assigned TP53 activated
status.
Samples were annotated as TP53 lost
if they contained i) a TP53 hotspot mutation as defined by IARC TP53 database or the MSKCC cancer hotspots database147,148 (see also, Key Resources Table), ii) two TP53 alterations, including SNV, CNV or SV, indicative of probable bi-allelic alterations; iii) one TP53 somatic alteration, including SNV, CNV, or SV or a germline TP53 mutation indicated by the diagnosis of Li-Fraumeni syndrome (LFS)150, or iv) one germline TP53 mutation indicated by LFS and the TP53 classifier score for matched RNA-Seq was greater than 0.5.
Samples were annotated as TP53 lost
if they contained i) a TP53 hotspot mutation as defined by IARC TP53 database or the MSKCC cancer hotspots database148,149 (see also, Key Resources Table), ii) two TP53 alterations, including SNV, CNV or SV, indicative of probable bi-allelic alterations; iii) one TP53 somatic alteration, including SNV, CNV, or SV or a germline TP53 mutation indicated by the diagnosis of Li-Fraumeni syndrome (LFS)151, or iv) one germline TP53 mutation indicated by LFS and the TP53 classifier score for matched RNA-Seq was greater than 0.5.
Participant metadata included a reported gender. @@ -2552,8 +2553,8 @@
survival-analysis
analysis module)We calculated overall survival (OS) as days since initial diagnosis and performed several survival analyses on the OpenPBTA cohort using the survival
R package.
-We performed survival analysis for patients by HGG subtype using the Kaplan-Meier estimator151 and a log-rank test (Mantel-Cox test)pubmed:5910392? on the different HGG subtypes.
-Next, we used multivariate Cox (proportional hazards) regression analysis152 to model the following: a) tp53 scores + telomerase scores + extent of tumor resection + LGG group + HGG group
, in which tp53 scores
and telomerase scores
are numeric, extent of tumor resection
is categorical, and LGG group
and HGG group
are binary variables indicating whether the sample is in either broad histology grouping, b) tp53 scores + telomerase scores + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), and c) quantiseq cell type fractions + CD274 expression + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), in which quantiseq cell type fractions
and CD274 expression
are numeric.
tp53 scores + telomerase scores + extent of tumor resection + LGG group + HGG group
, in which tp53 scores
and telomerase scores
are numeric, extent of tumor resection
is categorical, and LGG group
and HGG group
are binary variables indicating whether the sample is in either broad histology grouping, b) tp53 scores + telomerase scores + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), and c) quantiseq cell type fractions + CD274 expression + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), in which quantiseq cell type fractions
and CD274 expression
are numeric.
OpenPBTA workflows repository | this project | -153 | +154 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
OpenPBTA analysis repository | this project | -154 | +155 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
OpenPBTA manuscript repository | @@ -3105,195 +3106,198 @@
Broad histology group | +OpenPBTA molecular subtype | +Patients | +Tumors | +
---|---|---|---|
Chordoma | +CHDM, conventional | +2 | +2 | +
Chordoma | +CHDM, poorly differentiated | +2 | +4 | +
Embryonal tumor | +CNS Embryonal, NOS | +13 | +13 | +
Embryonal tumor | +CNS HGNET-MN1 | +1 | +1 | +
Embryonal tumor | +CNS NB-FOXR2 | +2 | +3 | +
Embryonal tumor | +ETMR, C19MC-altered | +5 | +5 | +
Embryonal tumor | +ETMR, NOS | +1 | +1 | +
Embryonal tumor | +MB, Group3 | +14 | +14 | +
Embryonal tumor | +MB, Group4 | +48 | +49 | +
Embryonal tumor | +MB, SHH | +24 | +30 | +
Embryonal tumor | +MB, WNT | +10 | +10 | +
Ependymoma | +EPN, H3 K28 | +1 | +1 | +
Ependymoma | +EPN, ST RELA | +25 | +28 | +
Ependymoma | +EPN, ST YAP1 | +3 | +4 | +
High-grade glioma | +DMG, H3 K28 | +18 | +24 | +
High-grade glioma | +DMG, H3 K28, TP53 activated | +10 | +13 | +
High-grade glioma | +DMG, H3 K28, TP53 loss | +30 | +40 | +
High-grade glioma | +HGG, H3 G35 | +3 | +3 | +
High-grade glioma | +HGG, H3 G35, TP53 loss | +1 | +1 | +
High-grade glioma | +HGG, H3 wildtype | +26 | +31 | +
High-grade glioma | +HGG, H3 wildtype, TP53 activated | +5 | +5 | +
High-grade glioma | +HGG, H3 wildtype, TP53 loss | +14 | +21 | +
High-grade glioma | +HGG, IDH, TP53 activated | +1 | +2 | +
High-grade glioma | +HGG, IDH, TP53 loss | +1 | +1 | +
Low-grade glioma | +GNG, BRAF V600E | +13 | +13 | +
Low-grade glioma | +GNG, BRAF V600E, CDKN2A/B | +1 | +1 | +
Low-grade glioma | +GNG, FGFR | +1 | +1 | +
Low-grade glioma | +GNG, H3 | +1 | +1 | +
Low-grade glioma | +GNG, IDH | +1 | +2 | +
Low-grade glioma | +GNG, KIAA1549-BRAF | +5 | +5 | +
Low-grade glioma | +GNG, MYB/MYBL1 | +1 | +1 | +
Low-grade glioma | +GNG, NF1-germline | +1 | +1 | +
Low-grade glioma | +GNG, NF1-somatic, BRAF V600E | +1 | +1 | +
Low-grade glioma | +GNG, other MAPK | +4 | +4 | +
Low-grade glioma | +GNG, other MAPK, IDH | +1 | +1 | +
Low-grade glioma | +GNG, RTK | +2 | +3 | +
Low-grade glioma | +GNG, wildtype | +14 | +14 | +
Low-grade glioma | +LGG, BRAF V600E | +25 | +27 | +
Low-grade glioma | +LGG, BRAF V600E, CDKN2A/B | +5 | +5 | +
Low-grade glioma | +LGG, FGFR | +8 | +8 | +
Low-grade glioma | +LGG, IDH | +3 | +3 | +
Low-grade glioma | +LGG, KIAA1549-BRAF | +106 | +113 | +
Low-grade glioma | +LGG, KIAA1549-BRAF, NF1-germline | +1 | +1 | +
Low-grade glioma | +LGG, KIAA1549-BRAF, other MAPK | +1 | +1 | +
Low-grade glioma | +LGG, MYB/MYBL1 | +2 | +2 | +
Low-grade glioma | +LGG, NF1-germline | +6 | +6 | +
Low-grade glioma | +LGG, NF1-germline, CDKN2A/B | +1 | +1 | +
Low-grade glioma | +LGG, NF1-germline, FGFR | +1 | +2 | +
Low-grade glioma | +LGG, NF1-somatic | +2 | +2 | +
Low-grade glioma | +LGG, NF1-somatic, FGFR | +1 | +1 | +
Low-grade glioma | +LGG, NF1-somatic, NF1-germline, CDKN2A/B | +1 | +1 | +
Low-grade glioma | +LGG, other MAPK | +11 | +12 | +
Low-grade glioma | +LGG, RTK | +8 | +10 | +
Low-grade glioma | +LGG, RTK, CDKN2A/B | +1 | +1 | +
Low-grade glioma | +LGG, wildtype | +33 | +34 | +
Low-grade glioma | +SEGA, RTK | +1 | +1 | +
Low-grade glioma | +SEGA, wildtype | +10 | +11 | +
Mesenchymal non-meningothelial tumor | +EWS | +9 | +11 | +
Neuronal and mixed neuronal-glial tumor | +CNC | +2 | +2 | +
Neuronal and mixed neuronal-glial tumor | +EVN | +1 | +1 | +
Neuronal and mixed neuronal-glial tumor | +GNT, BRAF V600E | +1 | +1 | +
Neuronal and mixed neuronal-glial tumor | +GNT, KIAA1549-BRAF | +1 | +2 | +
Neuronal and mixed neuronal-glial tumor | +GNT, other MAPK | +1 | +1 | +
Neuronal and mixed neuronal-glial tumor | +GNT, other MAPK, FGFR | +1 | +1 | +
Neuronal and mixed neuronal-glial tumor | +GNT, RTK | +1 | +2 | +
Tumor of sellar region | +CRANIO, ADAM | +27 | +27 | +
+ | Total | +577 | +644 | +
We performed a comprehensive genomic analysis of somatic SNVs, CNVs, SVs, and fusions across all 1,074 PBTA tumors (N = 1,019 RNA-Seq, N = 918 WGS, N = 32 WXS/Panel) and 22 cell lines (N = 16 RNA-Seq, N = 22 WGS), from 943 patients, 833 with paired normal specimens (N = 801 WGS, N = 32 WXS/Panel). +Tumor purity across PBTA samples was high (median 76%), though we observed some cancer groups with lower purity, including SEGA, PXA, and teratoma (Figure S3A). +Unless otherwise noted, each analysis was performed for diagnostic tumors using one tumor per patient.
+ +SNV consensus calling (Figure S1 and Figure S2A-G) revealed, as expected, lower tumor mutation burden (TMB) (Figure S2H) in pediatric tumors compared to adult brain tumors from The Cancer Genome Atlas (TCGA) (Figure S2I), with hypermutant (> 10 Mut/Mb) and ultra-hypermutant (> 100 Mut/Mb) tumors23 only found within HGGs and embryonal tumors. +Figure 2 and Figure S3B depict oncoprints recapitulating known histology-specific driver genes in primary tumors across OpenPBTA histologies, and Table S2 summarizes all detected alterations across cancer groups.
+As expected, most (62%, 140/226) LGGs harbored a somatic alteration in BRAF, with canonical BRAF::KIAA1549 fusions as the major oncogenic driver24 (Figure 2A). +We observed additional mutations in FGFR1 (2%), PIK3CA (2%), KRAS (2%), TP53 (1%), and ATRX (1%) and fusions in NTRK2 (2%), RAF1 (2%), MYB (1%), QKI (1%), ROS1 (1%), and FGFR2 (1%), concordant with previous studies reporting near-universal upregulation of the RAS/MAPK pathway in LGGs20,24. +Indeed, gene set variant analysis (GSVA) revealed significant upregulation (ANOVA Bonferroni-corrected p < 0.01) of the KRAS signaling pathway in LGGs (Figure 5B).
++Most (N = 95) embryonal tumors were medulloblastomas from four characterized molecular subtypes (WNT, SHH, Group3, and Group 4; see Molecular Subtyping of CNS Tumors), as identified by subtype-specific canonical mutations (Figure 2B). +We detected canonical SMARCB1/SMARCA4 deletions or inactivating mutations in atypical teratoid rhabdoid tumors (ATRTs; Table S2) and C19MC amplification in ETMRs (displayed within “Other embryonal tumors” in Figure 2B)25–28.
+Across HGGs, TP53 (57%, 36/63) and H3F3A (54%, 34/63) were both most mutated and co-occurring genes (Figure 2A and C), followed by frequent mutations in ATRX (29%, 18/63) which is commonly mutated in gliomas29. +We observed recurrent amplifications and fusions in EGFR, MET, PDGFRA, and KIT, highlighting that these tumors leverage multiple oncogenic mechanisms to activate tyrosine kinases, as previously reported14,30,31. +GSVA showed upregulation (ANOVA Bonferroni-corrected p < 0.01) of DNA repair, G2M checkpoint, and MYC pathways as well as downregulation of the TP53 pathway (Figure 5B). +The two ultra-hypermutated tumors (> 100 Mutations/Mb) were from patients with mismatch repair deficiency syndrome13.
+We observed that 25% (15/60) of ependymomas were C11orf95::RELA (now, ZFTA::RELA) fusion-positive32 and 68% (21/31) of craniopharyngiomas contained CTNNB1 mutations (Figure 2D). +We observed somatic mutations or fusions in NF2 in 41% (7/17) of meningiomas, 5% (3/60) of ependymomas, and 25% (3/12) of schwannomas, as well as rare fusions in ERBB4, YAP1, and/or QKI in 10% (6/60) of ependymomas. +DNETs harbored alterations in MAPK/PI3K pathway genes, as was previously reported33, including FGFR1 (21%, 4/19), PDGFRA (10%, 2/19), and BRAF (5%, 1/19).
+We analyzed mutational co-occurrence across the OpenPBTA, using a single tumor from each patient (N = 668) with WGS. +The top 50 mutated genes (see STAR Methods for details) in primary tumors are shown in Figure 3 by tumor type (A, bar plots), with co-occurrence scores illustrated in the heatmap (B). +As expected, TP53 was the most frequently mutated gene across the OpenPBTA (8.7%, 58/668), significantly co-occurring with H3F3A (OR = 30.05, 95% CI: 14.5 - 62.3, q = 2.34e-16), ATRX (OR = 23.3, 95% CI: 9.6 - 56.3, q = 8.72e-9), NF1 (OR = 8.26, 95% CI: 3.5 - 19.4, q = 7.40e-5), and EGFR (OR = 17.5, 95% CI: 4.8 - 63.9, q = 2e-4), with all of these driven by HGGs and consistent with previous reports30,34,35.
+In embryonal tumors, CTNNB1 mutations significantly co-occurred with TP53 mutations (OR = 43.6 95% CI: 7.1 - 265.8, q = 1.52e-3) as well as with DDX3X mutations (OR = 21.4, 95% CI: 4.7 - 97.9, q = 4.15e-3), events driven by medulloblastomas as previously reported36,37. +FGFR1 and PIK3CA mutations significantly co-occurred in LGGs (OR = 77.25, 95% CI: 10.0 - 596.8, q = 3.12e-3), consistent with previous findings37,38. +Of HGG tumors with TP53 or PPM1D mutations, 53/55 (96.3%) had mutations in only one of these genes (OR = 0.17, 95% CI: 0.04 - 0.89, q = 0.056), recapitulating previous observations that these mutations are usually mutually exclusive in HGGs39.
+CNV and SV analyses revealed that HGG, DMG, and medulloblastoma tumors had the most unstable genomes, while craniopharyngiomas and schwannomas generally lacked somatic CNV (Figure S3C). +These CNV patterns largely aligned with our TMB estimates (Figure S2H). +SV and CNV breakpoint densities were significantly correlated (linear regression p = 1.05e-38; Figure 3C), and as expected, the number of chromothripsis regions called increased with breakpoint density (Figure S3D-E). +We identified chromothripsis events in 31% (N = 12/39) of DMGs and in 44% (N = 21/48) of other HGGs (Figure 3D), and found evidence of chromothripsis in over 15% of sarcomas, PXAs, metastatic secondary tumors, chordomas, glial-neuronal tumors, germinomas, meningiomas, ependymomas, medulloblastomas, ATRTs, and other embryonal tumors.
+We assessed the contributions of eight adult CNS-specific mutational signatures from the RefSig database40 across tumors (Figure 3E and Figure S4A). +Signature 1, which reflects normal spontaneous deamination of 5-methylcytosine, predominated in stage 0 and/or 1 tumors characterized by low TMBs (Figure S2H) such as pilocytic astrocytomas, gangliogliomas, other LGGs, and craniopharyngiomas (Figure S4A). +Signature 1 weights were generally higher in tumors sampled at diagnosis (pre-treatment) compared to tumors from later phases of therapy (Figure S4B). +This trend may have emerged from therapy-induced mutations that produced additional signatures (e.g., temozolomide treatment has been suggested to drive Signature 1141), subclonal expansion, and/or acquisition of additional driver mutations during tumor progression, leading to detection of additional signatures. +We observed the CNS-specific signature N6 in nearly all tumors. +Signature 18 drivers (TP53, APC, NOTCH1; found at https://signal.mutationalsignatures.com/explore/referenceCancerSignature/31/drivers) are also canonical medulloblastoma drivers, and indeed, Signature 18 had the highest signature weight in medulloblastomas. +Finally, signatures 3, 8, 18, and MMR2 were prevalent in HGGs, including DMGs.
+ +Most RNA-Seq samples in the PBTA were prepared with ribosomal RNA depletion followed by stranded sequencing (N = 977), while remaining samples were prepared with poly-A selection (N = 58). +Since batch correction was not feasible (see Limitations of the Study and Figure S7A), the following transcriptomic analyses considered only stranded samples.
+We applied a TCGA-trained classifier42 to calculate a TP53 score, a proxy for TP53 gene or pathway dysregulation, and subsequently infer tumor TP53 inactivation status. +We identified “true positive” TP53 alterations from high-confidence SNVs, CNVs, SVs, and fusions in TP53, annotating tumors as “activated” if they harbored one of p.R273C or p.R248W gain-of-function mutations43, or “lost” if 1) the patient had a Li Fraumeni Syndrome (LFS) predisposition diagnosis, 2) the tumor harbored a known hotspot mutation, or 3) the tumor contained two hits (e.g. both SNV and CNV), suggesting both alleles were affected. +If the TP53 mutation did not reside within the DNA-binding domain or no alterations in TP53 were detected, we annotated the tumor as “other,” indicating an unknown TP53 alteration status. +The classifier achieved a high accuracy (AUROC = 0.86) for rRNA-depleted, stranded tumors, but it did not perform as well on the poly-A tumors in this cohort (AUROC = 0.62; Figure S5A).
+We observed that “activated” and “lost” tumors had similar TP53 scores (Figure 4B, Wilcoxon p = 0.92), contrasting our expectation that “lost” tumors would have higher TP53 scores. +This difference suggests that classifier scores > 0.5 may actually represent an oncogenic, or altered, TP53 phenotype rather than solely TP53 inactivation, as interpreted previously42. +However, “activated” tumors showed higher TP53 expression compared to those with TP53 “loss” mutations (Wilcoxon p = 0.006, Figure 4C). +DMGs, medulloblastomas, HGGs, DNETs, ependymomas, and craniopharyngiomas, all known to harbor TP53 mutations, had the highest median TP53 scores (Figure 4D). +By contrast, gangliogliomas, LGGs, meningiomas, and schwannomas had the lowest median scores.
+We hypothesized that tumors (N = 10) from patients with LFS (N = 8) would have higher TP53 scores, which we indeed observed for 8/10 tumors (Table S3).
+Although two tumors had low TP53 scores (BS_DEHJF4C7
at 0.09 and BS_ZD5HN296
at 0.28), pathology reports confirmed that both patients were diagnosed with LFS and harbored a TP53 pathogenic germline variant.
+These two LFS tumors also had low tumor purity (16% and 37%, respectively), suggesting that accurate classification may require a certain level of tumor content.
+We suggest this classifier could be generally applied to infer TP53 function in the absence of a predicted oncogenic TP53 alteration or DNA sequencing.
We used gene expression data to predict telomerase activity using EXpression-based Telomerase ENzymatic activity Detection (EXTEND
)44 as a surrogate measure of malignant potential44,45, where higher EXTEND
scores indicate higher telomerase activity.
+Aggressive tumors such as DMGs, other HGGs, and MB had high EXTEND
scores (Figure 4D), and low-grade lesions such as schwannomas, GNGs, DNETs, and other LGGs had among the lowest scores (Table S3), supporting previous reports that aggressive tumor phenotypes have higher telomerase activity46–49.
+While EXTEND
scores were not significantly higher in tumors with TERT promoter (TERTp) mutations (N = 6; Wilcoxon p-value = 0.1196), scores were significantly correlated with TERC (R = 0.619, p < 0.01) and TERT (R = 0.491, p < 0.01) log2 FPKM expression values (Figure S5B-C).
+Since catalytically-active telomerase requires full-length TERT, TERC, and certain accessory proteins50, we expect that EXTEND
scores may not be exclusively correlated with TERT alterations and expression.
We investigated the mutational signature profiles of hypermutant (TMB > 10 Mut/Mb; N = 3) and ultra-hypermutant (TMB > 100 Mut/Mb; N = 4) tumors and/or derived cell lines from six patients in OpenPBTA (Figure 4E).
+Five tumors were HGGs and one was a brain metastasis of a MYCN non-amplified neuroblastoma tumor.
+Signature 11, which is associated with exposure to temozolomide plus MGMT promoter and/or mismatch repair deficiency51, was indeed present in tumors with previous exposure to the drug (Table 2).
+We detected the MMR2 signature in tumors of four patients (PT_0SPKM4S8, PT_3CHB9PK5, PT_JNEV57VK, and PT_VTM2STE3) diagnosed with either constitutional mismatch repair deficiency (CMMRD) or Lynch syndrome (Table 2), genetic predisposition syndromes caused by a variant in a mismatch repair gene such as PMS2, MLH1, MSH2, MSH6, or others52.
+Three of these patients harbored pathogenic germline variants in one of the aforementioned genes.
+While we did not detect a known pathogenic variant in the germline of PT_VTM2STE3, this patient’s pathology report contained a self-reported PMS2 variant, and we indeed found 19 intronic variants of unknown significance (VUS) in their PMS2.
+This is not surprising since an estimated 49% of germline PMS2 variants in patients with CMMRD and/or Lynch syndrome are VUS52.
+Interestingly, while the cell line derived from patient PT_VTM2STE3’s tumor at progression was not hypermutated (TMB = 5.7 Mut/Mb), it only contained the MMR2 signature, suggesting selective pressure to maintain a mismatch repair (MMR) phenotype in vitro.
+Only one of the two cell lines derived from patient PT_JNEV57VK
’s progressive tumor was hypermutated (TMB = 35.9 Mut/Mb).
+The hypermutated cell line was strongly weighted towards signature 11, while the non-hypermutated cell line showed several lesser signature weights (1, 11, 18, 19, MMR2; Table S2).
+This mutational process plasticity highlights the importance of careful genomic characterization and model selection for preclinical studies.
Signature 18, which has been associated with high genomic instability and can induce a hypermutator phenotype40, was uniformly represented among hypermutant solid tumors. +Additionally, all hypermutant HGG tumors or cell lines had dysfunctional TP53 (Table 2), consistent with previous findings that tumors with high genomic instability signatures require TP53 dysregulation40. +With one exception, hypermutant and ultra-hypermutant tumors had high TP53 scores (> 0.5) and telomerase activity. +Interestingly, none of the hypermutant tumors showed evidence of signature 3 (present in homologous recombination deficient tumors), signature 8 (arises from double nucleotide substitutions/unknown etiology), or signature N6 (a universal CNS tumor signature). +The mutual exclusivity of signatures 3 and MMR2 corroborates previous suggestions that tumors do not generally feature both deficient homologous repair and mismatch repair42.
+Kids First Participant ID | +Kids First Biospecimen ID | +CBTN ID | +Phase of therapy | +Composition | +Therapy post-biopsy | +Cancer predisposition | +Pathogenic germline variant | +TMB | +OpenPBTA molecular subtype | +
---|---|---|---|---|---|---|---|---|---|
PT_0SPKM4S8 | +BS_VW4XN9Y7 | +7316-2640 | +Initial CNS Tumor | +Solid Tissue | +Radiation, Temozolomide, CCNU | +None documented | +NM_000535.7(PMS2):c.137G>T (p.Ser46Ile) (LP) | +187.4 | +HGG, H3 wildtype, TP53 activated | +
PT_3CHB9PK5 | +BS_20TBZG09 | +7316-515 | +Initial CNS Tumor | +Solid Tissue | +Radiation, Temozolomide, Irinotecan, Bevacizumab | +CMMRD | +NM_000179.3(MSH6):c.3439-2A>G (LP) | +307 | +HGG, H3 wildtype, TP53 loss | +
PT_3CHB9PK5 | +BS_8AY2GM4G | +7316-2085 | +Progressive | +Solid Tissue | +Radiation, Temozolomide, Irinotecan, Bevacizumab | +CMMRD | +NM_000179.3(MSH6):c.3439-2A>G (LP) | +321.6 | +HGG, H3 wildtype, TP53 loss | +
PT_EB0D3BXG | +BS_F0GNWEJJ | +7316-3311 | +Progressive | +Solid Tissue | +Radiation, Nivolumab | +None documented | +None detected | +26.3 | +Metastatic NBL, MYCN non-amplified | +
PT_JNEV57VK | +BS_85Q5P8GF | +7316-2594 | +Initial CNS Tumor | +Solid Tissue | +Radiation, Temozolomide | +Lynch Syndrome | +NM_000251.3(MSH2):c.1906G>C (p.Ala636Pro) (P) | +4.7 | +DMG, H3 K28, TP53 loss | +
PT_JNEV57VK | +BS_HM5GFJN8 | +7316-3058 | +Progressive | +Derived Cell Line | +Radiation, Temozolomide, Nivolumab | +Lynch Syndrome | +NM_000251.3(MSH2):c.1906G>C (p.Ala636Pro) (P) | +35.9 | +DMG, H3 K28, TP53 loss | +
PT_JNEV57VK | +BS_QWM9BPDY | +7316-3058 | +Progressive | +Derived Cell Line | +Radiation, Temozolomide, Nivolumab | +Lynch Syndrome | +NM_000251.3(MSH2):c.1906G>C (p.Ala636Pro) (P) | +7.4 | +DMG, H3 K28, TP53 loss | +
PT_JNEV57VK | +BS_P0QJ1QAH | +7316-3058 | +Progressive | +Solid Tissue | +Radiation, Temozolomide, Nivolumab | +Lynch Syndrome | +NM_000251.3(MSH2):c.1906G>C (p.Ala636Pro) (P) | +6.3 | +DMG, H3 K28, TP53 activated | +
PT_S0Q27J13 | +BS_P3PF53V8 | +7316-2307 | +Initial CNS Tumor | +Solid Tissue | +Radiation, Temozolomide, Irinotecan | +None documented | +None detected | +15.5 | +HGG, H3 wildtype, TP53 activated | +
PT_VTM2STE3 | +BS_ERFMPQN3 | +7316-2189 | +Progressive | +Derived Cell Line | +Unknown | +Lynch Syndrome | +None detected | +5.7 | +HGG, H3 wildtype, TP53 loss | +
PT_VTM2STE3 | +BS_02YBZSBY | +7316-2189 | +Progressive | +Solid Tissue | +Unknown | +Lynch Syndrome | +None detected | +274.5 | +HGG, H3 wildtype, TP53 activated | +
Next, we asked whether transcriptomic classification of TP53 dysregulation and/or telomerase activity recapitulate these oncogenic biomarkers’ known prognostic influence.
+We identified several expected trends, including a significant overall survival benefit following full tumor resection (HR = 0.35, 95% CI = 0.2 - 0.62, p < 0.001) or if the tumor was an LGG (HR = 0.046, 95% CI = 0.0062 - 0.34, p = 0.003), and a significant risk if the tumor was an HGG (HR = 6.2, 95% CI = 4.0 - 9.5, p < 0.001) (Figure 4F; STAR Methods).
+High telomerase scores were associated with poor prognosis across brain tumor histologies (HR = 20, 95% CI = 6.4 - 62, p < 0.001), demonstrating that EXTEND
scores calculated from RNA-Seq are an effective rapid surrogate measure for telomerase activity.
+Higher TP53 scores were associated with significant survival risks (Table S4) within DMGs (HR = 6436, 95% CI = 2.67 - 1.55e7, p = 0.03) and ependymomas (HR = 2003, 95% CI = 9.9 - 4.05e5, p = 0.005).
+Given this result, we next assessed whether different HGG molecular subtypes carry different survival risks if stratified by TP53 status.
+We found that DMG H3 K28 tumors with TP53 loss had significantly worse prognosis (HR = 2.8, CI = 1.4-5.6, p = 0.003) than those with wildtype TP53 (Figure 4G and Figure 4H), recapitulating results from two recent restrospective analyses of DIPG tumors10,53.
UMAP visualization of gene expression variation across brain tumors (Figure 5A) showed expected histological clustering of brain tumors. +We further observed that, except for three outliers, C11orf95::RELA (ZFTA::RELA) fusion-positive ependymomas fell within distinct clusters (Figure S6A). +Medulloblastoma (MB) tumors clustered by molecular subtype, with WNT and SHH in distinct clusters and Groups 3 and 4 showing some expected overlap (Figure S6B). +Notably, two MB tumors annotated as SHH did not cluster with the other MB tumors and one clustered with Group 3/4 tumors, suggesting potential subtype misclassification or different underlying biology of these two tumors. +BRAF-driven LGGs (Figure S6C) fell into three separate clusters, suggesting additional shared biology within each cluster. +Histone H3 G35-mutant HGGs generally clustered together and away from K28-mutant tumors (Figure S6D). +Interestingly, although H3 K28-mutant and H3 wildtype tumors have different biological drivers54, they did not form distinct clusters. +This pattern suggests these subtypes may be driven by common transcriptional programs, have other much stronger biological drivers than their known distinct epigenetic drivers, or we lack power to detect transcriptional differences.
+We performed GSVA for Hallmark cancer gene sets (Figure 5B) and quantified immune cell fractions using quanTIseq (Figure 5C and Figure S6E), results from which recapitulated previously-described tumor biology. +For example, HGG, DMG, MB, and ATRT tumors are known to upregulate MYC55 which in turn activates E2F and S phasepubmed:11511364?. +Indeed, we detected significant (Bonferroni-corrected p < 0.05) upregulation of MYC and E2F targets, as well as G2M (cell cycle phase following S phase) in MBs, ATRTs, and HGGs compared to several other cancer groups. +In contrast, LGGs showed significant downregulation (Bonferroni-corrected p < 0.05, multiple cancer group comparisons) of these pathways. +Schwannomas and neurofibromas, which have an inflammatory immune microenvironment of T and B lymphocytes and tumor-associated macrophages (TAMs), are driven by upregulation of cytokines such as IFN\(\gamma\), IL-1, and IL-6, and TNF\(\alpha\)56. +GSVA revealed significant upregulation of these cytokines in hallmark pathways (Bonferroni-corrected p < 0.05, multiple cancer group comparisons) (Figure 5B), and monocytes dominated these tumors’ immune cell repertoire (Figure 5C). +We also observed significant upregulation of pro-inflammatory cytokines IFN\(\alpha\) and IFN\(\gamma\) in both LGGs and craniopharyngiomas when compared to either medulloblastoma or ependymomas (Bonferroni-corrected p < 0.05) (Figure 5B). +Together, these results support previous proteogenomic findings that aggressive medulloblastomas and ependymomas have lower immune infiltration compared to BRAF-driven LGGs and craniopharyngiomas57.
+Although CD8+ T-cell infiltration across all cancer groups was minimal (Figure 5C), we observed signal in specific cancer molecular subtypes (Groups 3 and 4 medulloblastoma) as well as outlier tumors (BRAF-driven LGG, BRAF-driven and wildtype ganglioglioma, and CNS embryonal NOS; Figure S6E) +Surprisingly, the classically immunologically-cold HGGs and DMGs58,59 contained higher overall fractions of immune cells, primarily monocytes, dendritic cells, and NK cells (Figure 5C). +Thus, quanTIseq might have actually captured microglia within these immune cell fractions.
+While we did not detect notable prognostic effects of immune cell infiltration on overall survival in HGGs or DMGs, we found that high levels of macrophage M1 and monocytes were associated with poorer overall survival (monocyte HR = 2.1e18, 95% CI = 3.80e5 - 1.2e31, p = 0.005, multivariate Cox) in medulloblastomas (Figure 5D). +We further reproduced previous findings (Figure 5E) that medulloblastomas typically have low expression of CD274 (PD-L1)60. +We also found that higher expression of CD274 was significantly associated with improved overall prognosis for medulloblastoma tumors, although marginal (HR = 0.0012, 95% CI = 7.5e−06 - 0.18, p = 0.008, multivariate Cox) (Figure 5D). +This result may be explained by the higher expression of CD274 observed in WNT subtype tumors by us and others61, as this diagnosis carries the best prognosis of all medulloblastoma subgroups (Figure 5E).
+We additionally explored the ratio of CD8+ to CD4+ T cells across tumor subtypes. +This ratio has been associated with better immunotherapy response and prognosis following PD-L1 inhibition in non-small cell lung cancer or adoptive T cell therapy in multiple stage III or IV cancers62,63. +While adamantinomatous craniopharyngiomas and Group 3 and Group 4 medulloblastomas had the highest ratios (Figure S6F), very few tumors had ratios greater than 1, highlighting an urgent need to identify novel therapeutics for pediatric brain tumors with poor prognosis.
+Finally, we explored the potential influence of tumor purity by repeating selected transcriptomic analyses restricted to only samples with high tumor purity (see STAR Methods). +Results from these analyses were broadly consistent (Figure S7D-I) with results derived from all stranded RNA-Seq samples.
+The CBTN released the PBTA raw genomic data in September 2018 without embargo, allowing researchers immediate access to begin making discoveries on behalf of children with CNS tumors everywhere. +Since this publication, the CBTN has approved over 200 data research projects4 from 69 different institutions, with 60% from non-CBTN sites. +We created OpenPBTA as an open, real-time, reproducible analysis framework to genomically characterize pediatric brain tumors, bringing together basic and translational researchers, clinicians, and data scientists. +We provide reusable code and data resources, paired with cloud-based availability of source and derived data resources, to the pediatric oncology community, encouraging interdisciplinary collaboration. +To our knowledge, this initiative represents the first large-scale, collaborative, open analysis of genomic data coupled with open manuscript writing, wherein we comprehensively analyzed the PBTA cohort. +Using available WGS, WXS, and RNA-Seq data, we generated high-confidence consensus SNV and CNV calls, prioritized putative oncogenic fusions, and established over 40 scalable and rigorously-reviewed modules to perform common downstream cancer genomics analyses. +We detected expected patterns of genomic lesions, mutational signatures, and aberrantly regulated signaling pathways across multiple pediatric brain tumor histologies.
+Assembling large, pan-histology cohorts of fresh frozen samples and associated clinical phenotypes and outcomes requires a multi-year, multi-institutional framework, like those provided by CBTN and PNOC.
+As such, uniform clinical molecular subtyping was largely not performed for this cohort at the time of sample collection.
+Since DNA methylation data for these samples were not yet available to classify molecular subtypes, we created RNA- and DNA-based subtyping modules aligned with WHO molecularly-defined diagnoses.
+We worked closely with pathologists and clinicians to assign research-grade integrated diagnoses for 60% of tumors while discovering incorrectly diagnosed or mis-identified samples in the OpenPBTA cohort.
+For example, we subtyped medulloblastoma tumors, of which only 35% (43/122) had prior subtype information from pathology reports, using MMS2
or MedulloClassifier
21,22 and subsequently applied the consensus of these methods to subtype all medulloblastomas.
We advanced the integrative analyses and cross-cohort comparison via a number of validated modules. +We used an expression classifier to determine whether tumors have dysfunctional TP5342 and the EXTEND algorithm to determine their degree of telomerase activity using a 13-gene signature44. +Interestingly, we found that hypermutant HGGs universally displayed TP53 dysregulation, unlike adult cancers like colorectal cancer and gastric adenocarcinoma where TP53 dysregulation in hypermutated tumors is less common64,65. +Furthermore, high TP53 scores were a significant prognostic marker for poor overall survival for patients with tumor types including H3 K28-mutant DMGs and ependymomas. +We also show that EXTEND scores are a robust surrogate measure for telomerase activity in pediatric brain tumors. +By assessing TP53 and telomerase activity prospectively from expression data, information usually only attainable with DNA sequencing and/or qPCR, we incorporated oncogenic biomarker and prognostic knowledge thereby expanding our biological understanding of these tumors.
+We identified enrichment of hallmark cancer pathways and characterized the immune cell landscape across pediatric brain tumors, demonstrating tumors in some histologies, such as schwannomas, craniopharyngiomas, and low-grade gliomas, may have a inflammatory tumor microenvironment. +Notably, we observed upregulation of IFN\(\gamma\), IL-1, and IL-6, and TNF\(\alpha\) in craniopharyngiomas, tumors difficult to resect due to their anatomical location and critical surrounding structures. +Neurotoxic side effects have been reported in response to IFN\(\alpha\) immunotherapy66,67, leading researchers to propose additional immune vulnerabilities, such as IL-6 inhibition and immune checkpoint blockade, as cystic adamantinomatous craniopharyngiomas therapies68–70,pubmed:34966342?,pubmed:32075140?. +Our results support this endeavor. +Finally, we reproduced the overall known poor infiltration of CD8+ T cells and general low expression of CD274 (PD-L1) in pediatric brain tumors, highlighting that we urgently need novel therapeutic strategies for tumors unlikely to respond to immune checkpoint blockade therapy.
+While large-scale collaborative efforts may take a longer time to complete, adoption an open science framework substantially mitigated this concern. +By maintaining all data, analytical code, and results in public repositories, we ensured that such logistics did not hinder progress in pediatric cancer research. +Indeed, OpenPBTA is already a foundational data analysis and processing layer for several discovery research and translational projects which will continue to add other genomic modalities and analyses, including germline, epigenomic, single-cell, splicing, imaging, and model drug response data. +For example, the OpenPBTA RNA fusion filtering module led to the development of the R package annoFuse71 and an R Shiny application shinyFuse. +Leveraging OpenPBTA’s medulloblastoma subtyping and immune deconvolution analyses, Dang and colleagues showed that SHH tumors are enriched with monocyte and microglia-derived macrophages, which may accumulate following radiation therapy9. +Expression and CNV analyses demonstrated that GPC2 is a highly expressed and copy-number gained immunotherapeutic target in ETMRs, medulloblastomas, choroid plexus carcinomas, H3 wildtype high-grade gliomas, and DMGs. +Foster and colleagues therefore developed a chimeric antigen receptor (CAR) directed against GPC2, which shows preclinical efficacy in mouse models11. +Another study harnessed OpenPBTA to integrate germline variants, discovering that pediatric HGG patients with alternative telomere lengthening are enriched for pathogenic or likely pathogenic germline variants in the MMR pathway, possess oncogenic ATRX mutations and have increased TMB12. +Moreover, OpenPBTA has enabled a framework to support real-time integration of clinical trial subjects as they enrolled on the PNOC008 high-grade glioma clinical trial72 or PNOC027 medulloblastoma clinical trial73, allowing researchers and clinicians to link tumor biology to translational impact through clinical decision support during tumor board discussions. +Finally, as part of the NCI’s CCDI, OpenPBTA was recently expanded into OpenPedCan, a pan-pediatric cancer effort (https://github.com/PediatricOpenTargets/OpenPedCan-analysis) which enabled creation of the pediatric Molecular Targets Platform (https://moleculartargets.ccdi.cancer.gov/) in support of the RACE Act. +An additional, large-scale cohort of >1,500 tumor samples and associated germline DNA is undergoing harmonization as part of CBTN CCDI-Kids First NCI and Common Fund project (https://commonfund.nih.gov/kidsfirst/2021X01projects#FY21_Resnick) and will be immediately integrated with OpenPBTA data through OpenPedCan. +OpenPBTA has paved the way for new modes of collaborative data-driven discovery using open, reproducible, and scalable analyses that will continue to grow over time. +We anticipate this foundational work will have an ongoing, long-term impact for pediatric oncology researchers, ultimately accelerating translation and leading to improved outcomes for children with cancer.
+All code and processed data are openly available through GitHub, CAVATICA, Zenodo, and PedcBioPortal (see STAR METHODS).
+Notably, PBTA brain tumor samples were collected over decades, and RNA samples were prepared using two distinct library preparations (stranded or poly-A, Figure S7A) by multiple sequencing centers. +While we noted a strong library preparation batch effect (Figure S7B) and a possible sequencing center batch effect (Figure S7C), cancer groups are highly unbalanced across library preparations (Figure S7A). +We did not perform batch correction because removing batch effects across unbalanced groups may induce false differences among groups74,75. +Instead, we circumvent batch effects by grouping only stranded RNA-Seq expression data, which comprises the vast majority of the PBTA cohort, for transcriptomic analyses presented in Figure ?? and Figure ?? . +As batch correction strategy depends highly on research goals75, we provide library preparation-specific expression matrices in the OpenPBTA data release for others to adapt to their needs. +A second potential limitation is that performing analyses with all samples, rather than samples with high tumor purity, might result in loss of information, such as subclonal variants or low-level oncogenic pathway expression. +To this end, we re-performed transcriptomic analyses using only samples with high tumor purity (see Methods for details), and indeed, results were broadly consistent with those derived from the full cohort (Figure S7D-I). +To enable more robust statistical analysis and presentation of results, we randomly selected one independent specimen from patients with duplicate sequenced samples per tumor event rather than combining the data. +This practice did not induce notable differences if the selected specimen changed over time, e.g., with a new data release. +Finally, because this initial PBTA cohort mostly contains samples collected at diagnosis from one tumor section/punch, we could not reliably perform systematic intratumoral and/or longitudinal analyses, though we expect nearly 100 paired longitudinal tumors from the (NIH X01 CA267587-01 pediatric brain tumor cohort) to be released through OpenPedCan for future exploration.
+We graciously thank the patients and families who have donated tumors to CBTN and/or PNOC, without which this research would not be possible.
+Philanthropic support has ensured the CBTN’s ability to collect, store, manage, and distribute specimen and data. +The following donors have provided leadership level support for CBTN: CBTN Executive Council members, Brain Tumor Board of Visitors, Children’s Brain Tumor Foundation, Easie Family Foundation, Kortney Rose Foundation, Lilabean Foundation, Minnick Family Charitable Fund, Perricelli Family, Psalm 103 Foundation, and Swifty Foundation.
+This work was funded through the Alex’s Lemonade Stand Foundation (ALSF) Childhood Cancer Data Lab (CSG), ALSF Young Investigator Award (JLR), ALSF Catalyst Award (JLR, ACR, PBS), ALSF Catalyst Award (SJS), ALSF CCDL Postdoctoral Training Grant (SMF), Children’s Hospital of Philadelphia Division of Neurosurgery (PBS and ACR), Australian Government, Department of Education (APH), St. Anna Kinderkrebsforschung, Austria (ARP), the Mildred Scheel Early Career Center Dresden P2, funded by the German Cancer Aid (ARP), NIH Grants 3P30 CA016520-44S5 (ACR), U2C HL138346-03 (ACR, APH), U24 CA220457-03 (ACR), K12GM081259 (SMF), R03-CA23036 (SJD), NIH Contract Nos. HHSN261200800001E (SJD) and 75N91019D00024, Task Order No. 75N91020F00003 (JLR, ACR, APH), Intramural Research Program of the Division of Cancer Epidemiology and Genetics of the National Cancer Institute +The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products or organizations imply endorsement by the U.S. Government.
+The authors thank the following collaborators who contributed or supervised analyses present in the analysis repository that were not included in the manuscript: William Amadio, Holly C. Beale, Ellen T. Kephart, A. Geoffrey Lyle, and Olena M. Vaske. +Finally, we thank Yuanchao Zhang for adding to the project codebase, Jessica B. Foster for helpful discussions while drafting the manuscript, and Gina D. Mawla for identifying and reporting OpenPBTA data issues.
+Author | +Contributions | +
---|---|
Joshua A. Shapiro | +Methodology, Software, Validation, Formal analysis, Investigation, Writing - Original draft, Writing - Review and editing, Visualization, Supervision | +
Krutika S. Gaonkar | +Data curation, Formal analysis, Investigation, Methodology, Software, Writing – Original draft, Writing - Review and editing | +
Stephanie J. Spielman | +Validation, Formal analysis, Writing - Original draft, Writing - Review and editing, Investigation, Software, Visualization, Supervision, Funding acquisition | +
Candace L. Savonen | +Methodology, Software, Validation, Formal analysis, Investigation, Writing - Original draft, Writing - Review and editing, Visualization | +
Chante J. Bethell | +Methodology, Validation, Formal analysis, Investigation, Writing - Original draft, Visualization | +
Run Jin | +Data curation, Formal analysis, Visualization, Writing – Original draft, Writing - Review and editing | +
Komal S. Rathi | +Formal analysis, Investigation, Methodology, Writing – Original draft | +
Yuankun Zhu | +Data curation, Formal analysis, Investigation, Methodology, Supervision | +
Laura E. Egolf | +Formal analysis, Writing - Original draft | +
Bailey K. Farrow | +Data curation, Software | +
Daniel P. Miller | +Formal analysis | +
Yang Yang | +Formal analysis, Software | +
Tejaswi Koganti | +Formal analysis, Investigation | +
Nighat Noureen | +Formal analysis, Visualization, Writing - Original draft | +
Mateusz P. Koptyra | +Formal analysis, Writing – Original draft | +
Nhat Duong | +Formal analysis, Investigation, Methodology | +
Mariarita Santi | +Investigation, Validation, Writing - Review and editing | +
Jung Kim | +Investigation, Writing - Review and editing | +
Shannon Robins | +Data curation | +
Phillip B. Storm | +Conceptualization, Funding acquisition, Resources | +
Stephen C. Mack | +Writing - Review and editing | +
Jena V. Lilly | +Conceptualization, Funding acquisition, Project administration | +
Hongbo M. Xie | +Methodology, Supervision | +
Payal Jain | +Data curation, Investigation, Validation | +
Pichai Raman | +Conceptualization, Formal analysis, Methodology | +
Brian R. Rood | +Conceptualization | +
Rishi R. Lulla | +Conceptualization | +
Javad Nazarian | +Conceptualization | +
Adam A. Kraya | +Methodology | +
Zalman Vaksman | +Formal analysis, Investigation | +
Allison P. Heath | +Project administration, Funding acquisition | +
Cassie Kline | +Supervision, Investigation, Writing - Review and editing | +
Laura Scolaro | +Data curation | +
Angela N. Viaene | +Investigation, Validation | +
Xiaoyan Huang | +Formal analysis | +
Gregory P. Way | +Investigation, Writing - Review and editing | +
Steven M. Foltz | +Validation, Funding acquisition | +
Bo Zhang | +Data curation, Formal analysis | +
Anna R. Poetsch | +Formal analysis, Funding acquisition, Writing – Review and editing | +
Sabine Mueller | +Conceptualization | +
Brian M. Ennis | +Data curation, Formal analysis | +
Michael Prados | +Conceptualization | +
Sharon J. Diskin | +Investigation, Supervision, Validation, Funding acquisition, Writing - Review and editing | +
Siyuan Zheng | +Formal analysis, Visualization, Writing - Original draft, Supervision, Writing - Review and editing | +
Yiran Guo | +Formal analysis, Writing - Review and editing | +
Shrivats Kannan | +Formal analysis, Methodology, Writing – Original draft | +
Angela J. Waanders | +Supervision, Conceptualization | +
Ashley S. Margol | +Writing - Review and editing | +
Meen Chul Kim | +Data curation | +
Derek Hanson | +Validation | +
Nicholas Van Kuren | +Data curation, Software | +
Jessica Wong | +Writing – Original draft | +
Rebecca S. Kaufman | +Formal analysis, Investigation, Validation | +
Noel Coleman | +Data curation | +
Christopher Blackden | +Resources | +
Kristina A. Cole | +Writing - Review and editing | +
Jennifer L. Mason | +Supervision | +
Peter J. Madsen | +Writing – Review & editing | +
Carl J. Koschmann | +Conceptualization | +
Douglas R. Stewart | +Supervision, Writing - Review and editing | +
Eric Wafula | +Formal analysis, Software | +
Miguel A. Brown | +Data curation, Methodology, Formal analysis | +
Adam C. Resnick | +Conceptualization, Funding acquisition, Resources, Supervision | +
Casey S. Greene | +Conceptualization, Funding acquisition, Methodology, Project administration, Software, Supervision, Writing – Review & editing | +
Jo Lynne Rokita | +Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Software, Supervision, Writing – Original draft, Writing - Review and editing | +
Jaclyn N. Taroni | +Methodology, Software, Validation, Formal analysis, Investigation, Data curation, Writing - Review and editing, Visualization, Supervision, Project administration | +
Children’s Brain Tumor Network | +Conceptualization | +
Pacific Pediatric Neuro-Oncology Consortium | +Conceptualization | +
Except for the first and last four authors, authorship order was determined as follows: Authors who contributed to the OpenPBTA code base are listed based on number of modules included in the manuscript to which that individual contributed and, in the case of ties, a random order is used. All remaining authors are then listed in a random order.
+Code for determining authorship order can be found in the count-contributions
module of the OpenPBTA analysis repository.
CSG’s spouse was an employee of Alex’s Lemonade Stand Foundation, which was a sponsor of this research. +JAS, CLS, CJB, SJS, and JNT are or were employees of Alex’s Lemonade Stand Foundation, a sponsor of this research. +AJW is a member of the Scientific Advisory boards for Alexion and DayOne Biopharmaceuticals.
+Figure 1. Overview of the OpenPBTA Project. +A, CBTN and PNOC collected tumors from 943 patients. 22 tumor cell lines were created, and over 2000 specimens were sequenced (N = 1035 RNA-Seq, N = 940 WGS, and N = 32 WXS or targeted panel). The Kids First Data Resource Center Data harmonized the data using Amazon S3 through CAVATICA. Panel created with BioRender.com. B, Number of biospecimens across phases of therapy, with one broad histology per panel. Each bar denotes a cancer group. (Abbreviations: GNG = ganglioglioma, Other LGG = other low-grade glioma, PA = pilocytic astrocytoma, PXA = pleomorphic xanthoastrocytoma, SEGA = subependymal giant cell astrocytoma, DIPG = diffuse intrinsic pontine glioma, DMG = diffuse midline glioma, Other HGG = other high-grade glioma, ATRT = atypical teratoid rhabdoid tumor, MB = medulloblastoma, Other ET = other embryonal tumor, EPN = ependymoma, PNF = plexiform neurofibroma, DNET = dysembryoplastic neuroepithelial tumor, CRANIO = craniopharyngioma, EWS = Ewing sarcoma, CPP = choroid plexus papilloma). C, Overview of the open analysis and manuscript contribution models. Contributors proposed analyses, implemented it in their fork, and filed a pull request (PR) with proposed changes. PRs underwent review for scientific rigor and accuracy. Container and continuous integration technologies ensured that all software dependencies were included and code was not sensitive to underlying data changes. Finally, a contributor filed a PR documenting their methods and results to the Manubot-powered manuscript repository for review. D, A potential path for an analytical PR. Arrows indicate revisions.
+Figure 2. Mutational landscape of PBTA tumors.
+Frequencies of canonical somatic gene mutations, CNVs, fusions, and TMB (top bar plot) for the top mutated genes across primary tumors within the OpenPBTA dataset. A, LGGs (N = 226): pilocytic astrocytoma (N = 104), other LGG (N = 68), ganglioglioma (N = 35), pleomorphic xanthoastrocytoma (N = 9), subependymal giant cell astrocytoma (N = 10). B, Embryonal tumors (N = 129): medulloblastoma (N = 95), atypical teratoid rhabdoid tumor (N = 24), other embryonal tumor (N = 10). C, HGGs (N = 63): diffuse midline glioma (N = 36) and other HGG (N = 27). D, Other CNS tumors (N = 153): ependymoma (N = 60), craniopharyngioma (N = 31), meningioma (N = 17), dysembryoplastic neuroepithelial tumor (N = 19), Ewing sarcoma (N = 7), schwannoma (N = 12), and neurofibroma plexiform (N = 7). Rare CNS tumors are displayed in Figure S3B. Histology (Cancer Group
) and sex (Germline sex estimate
) annotations are displayed under each plot. Only tumors with mutations in the listed genes are shown. Multiple CNVs are denoted as a complex event. N denotes the number of unique tumors (one tumor per patient).
Figure 3. Mutational co-occurrence and signatures highlight key oncogenic drivers. +A, Nonsynonymous mutations for 50 most commonly-mutated genes across all histologies. “Other” denotes a histology with <10 tumors. B, Co-occurrence and mutual exclusivity of mutated genes. The co-occurrence score is defined as \(I(-\log_{10}(P))\) where \(P\) is Fisher’s exact test and \(I\) is 1 when mutations co-occur more often than expected or -1 when exclusivity is more common. C, Number of SV and CNV breaks are significantly correlated (Adjusted R = 0.443, p = 1.05e-38). D, Chromothripsis frequency across cancer groups with N >= 3 tumors. E, Sina plots of RefSig signature weights for signatures 1, 11, 18, 19, 3, 8, N6, MMR2, and Other across cancer groups. Boxplot represents 5% (lower whisker), 25% (lower box), 50% (median), 75% (upper box), and 95% (upper whisker) quantiles.
+Figure 4. TP53 and telomerase activity +A, Receiver Operating Characteristic for TP53 classifier run on stranded FPKM RNA-Seq. B, Violin and strip plots of TP53 scores plotted by TP53 alteration type (Nactivated = 11, Nlost = 100, Nother = 866). C, Violin and strip plots of TP53 RNA expression plotted by TP53 activation status (Nactivated = 11, Nlost = 100, Nother = 866). D, Boxplots of TP53 and telomerase (EXTEND) scores across cancer groups. TMB status is highlighted in orange (hypermutant) or red (ultra-hypermutant). E, Heatmap of RefSig mutational signatures for patients with at least one hypermutant tumor or cell line. F, Forest plot depicting prognostic effects of TP53 and telomerase scores on overall survival (OS), controlling for extent of tumor resection, LGG group, and HGG group. G, Forest plot depicting the effect of molecular subtype on HGG OS. Hazard ratios (HR) with 95% confidence intervals and p-values (multivariate Cox) are given in F and G. Black diamonds denote significant p-values, and gray diamonds denote reference groups. H, Kaplan-Meier curve of HGGs by molecular subtype. Boxplot represents 5% (lower whisker), 25% (lower box), 50% (median), 75% (upper box), and 95% (upper whisker) quantiles.
+Figure 5. Transcriptomic and immune landscape of pediatric brain tumors +A, First two dimensions of transcriptome data UMAP, with points colored by broad histology. B, Heatmap of GSVA scores for Hallmark gene sets with tumors ordered by cancer group. C, Boxplots of quanTIseq estimates of immune cell proportions in cancer groups with N > 15 tumors. Note: other HGGs and other LGGs have immune cell proportions similar to DMG and pilocytic astrocytoma, respectively, and are not shown. D, Forest plot depicting additive effects of CD274 expression, immune cell proportion, and extent of tumor resection on OS of medulloblastoma patients. HRs with 95% confidence intervals and p-values (multivariate Cox) are listed. Black diamonds denote significant p-values, and gray diamonds denote reference groups. Note: the Macrophage M1 HR was 0 (coefficient = -9.90e+4) with infinite upper and lower CIs, and thus was not included in the figure. E, Boxplot of CD274 expression (log2 FPKM) for medulloblastomas grouped by subtype. Bonferroni-corrected p-values from Wilcoxon tests are shown. Boxplot represents 5% (lower whisker), 25% (lower box), 50% (median), 75% (upper box), and 95% (upper whisker) quantiles. Only stranded RNA-Seq data is plotted.
+Table 1. Molecular subtypes generated through the OpenPBTA project. +Broad tumor histologies, molecular subtypes generated, and number of patients and tumors subtyped within OpenPBTA.
+Table 2. Patients with hypermutant tumors. +Patients with at least one hypermutant or ultra-hypermutant tumor or cell line. Pathogenic (P) or likely pathogenic (LP) germline variants, coding region TMB, phase of therapy, therapeutic interventions, cancer predisposition (CMMRD = Constitutional mismatch repair deficiency), and molecular subtypes are included.
+Requests for access to OpenPBTA raw data and/or specimens may be directed to, and will be fulfilled by Jo Lynne Rokita (rokita@chop.edu).
+This study did not create new, unique reagents.
+Raw and harmonized WGS, WXS, and RNA-Seq data derived from human samples are available within the KidsFirst Portal76 upon access request to the CBTN (https://cbtn.org/) as of the date of the publication. +In addition, merged summary files are openly accessible at https://cavatica.sbgenomics.com/u/cavatica/openpbta or via download script in the https://github.com/AlexsLemonade/OpenPBTA-analysis repository. +Summary data are visible within PedcBioPortal at https://pedcbioportal.kidsfirstdrc.org/study/summary?id=openpbta. +Associated DOIs are listed in the Key Resources Table. +Data underlying manuscript figures are available on Zenodo77. +
+All original code was developed within the following repositories and is publicly available as follows. +Primary data analyses can be found at https://github.com/d3b-center/OpenPBTA-workflows. +Downstream data analyses can be found at https://github.com/AlexsLemonade/OpenPBTA-analysis. +Manuscript code can be found at https://github.com/AlexsLemonade/OpenPBTA-manuscript. +Associated DOIs are listed in the Key Resources Table. +Software versions are documented in Table S5 as an appendix to the Key Resources Table.
+Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request.
+We maintained a data release folder on Amazon S3, downloadable directly from S3 or our open-access CAVATICA project, with merged files for each analysis (See Data and code availability section). +As we produced new results (e.g., tumor mutation burden calculations) that we expected to be used across multiple analyses, or identified data issues, we created new data releases in a versioned manner. +We reran all manuscript-specific analysis modules with the latest data release (v23) prior to submission and subsequently created a GitHub repository-tagged release to ensure reproducibility.
+The Pediatric Brain Tumor Atlas specimens are comprised of samples from Children’s Brain Tumor Network (CBTN) and the Pediatric Pacific Neuro-Oncology Consortium (PNOC). +The CBTN is a collaborative, multi-institutional (32 institutions worldwide) research program dedicated to the study of childhood brain tumors. +PNOC is an international consortium dedicated to bringing new therapies to children and young adults with brain tumors. +We also include blood and tumor biospecimens from newly-diagnosed diffuse intrinsic pontine glioma (DIPG) patients as part of the PNOC003 clinical trial PNOC003/NCT0227498714.
+The CBTN-generated cell lines were derived from either fresh tumor tissue directly obtained from surgery performed at Children’s Hospital of Philadelphia (CHOP) or from prospectively collected tumor specimens stored in Recover Cell Culture Freezing medium (cat# 12648010, Gibco). +We dissociated tumor tissue using enzymatic method with papain as described13. +Briefly, we washed tissue with HBSS (cat# 14175095, Gibco), and we minced and incubated the tissue with activated papain solution (cat# LS003124, SciQuest) for up to 45 minutes. +We used ovomucoid solution (cat# 542000, SciQuest) to inactivate the papain, briefly treated tissue with DNase (cat# 10104159001, Roche), passed it through the 100μm cell strainer (cat# 542000, Greiner Bio-One). +We initiated two cell culture conditions based on the number of cells available. +For cultures utilizing the fetal bovine serum (FBS), we plated a minimum density of 3×105 cells/mL in DMEM/F-12 medium (cat# D8062, Sigma) supplemented with 20% FBS (cat# SH30910.03, Hyclone), 1% GlutaMAX (cat# 35050061, Gibco), Penicillin/Streptomycin-Amphotericin B Mixture (cat# 17-745E, Lonza), and 0.2% Normocin (cat# ant-nr-2, Invivogen). +For serum-free media conditions, we plated cells at minimum density of 1×106 cells/mL in DMEM/F12 medium supplemented with 1% GlutaMAX, 1X B-27 supplement minus vitamin A (cat# 12587-010, Gibco), 1x N-2 supplement (cat# 17502001, Gibco), 20 ng/ml epidermal growth factor (cat# PHG0311L, Gibco), 20 ng/mL basic fibroblast growth factor (cat# 100-18B, PeproTech), 2.5μg/mL heparin (cat# H3149, Sigma), Penicillin/Streptomycin-Amphotericin B Mixture, and 0.2% Normocin.
+The Translational Genomic Research Institute (TGEN; Phoenix, AZ) performed DNA and RNA extractions on tumor biopsies using a DNA/RNA AllPrep Kit (Qiagen, #80204). +All RNA used for library prep had a minimum RIN of seven, but no QC thresholds were implemented for the DNA. +For library preparation, 500 ng of nucleic acids were used as input for RNA-Seq, WXS, and targeted DNA panel (panel) sequencing. +RNA library preparation was performed using the TruSeq RNA Sample Prep Kit (Illumina, #FC-122-1001) with poly-A selection, and the exome prep was performed using KAPA Library Preparation Kit (Roche, #KK8201) using Agilent’s SureSelect Human All Exon V5 backbone with custom probes. +The targeted DNA panel developed by Ashion Analytics (formerly known as the GEM Cancer panel) consisted of exonic probes against 541 cancer genes. +Both panel and WXS assays contained 44,000 probes across evenly spaced genomic loci used for genome-wide copy number analysis. +For the panel, additional probes tiled across intronic regions of 22 known tumor suppressor genes and 22 genes involved in common cancer translocations for structural analysis. +All extractions and library preparations were performed according to manufacturer’s instructions.
+Blood, tissue, and cell line DNA/RNA extractions were performed at the Biorepository Core at CHOP. +Briefly, 10-20 mg frozen tissue, 0.4-1ml of blood, or 2e6 cells pellet was used for extractions. +Tissues were lysed using a Qiagen TissueLyser II (Qiagen) with 2×30 sec at 18Hz settings using 5 mm steel beads (cat# 69989, Qiagen). +Both tissue and cell pellets processes included a CHCl3 extraction and were run on the QIACube automated platform (Qiagen) using the AllPrep DNA/RNA/miRNA Universal kit (cat# 80224, Qiagen). +Blood was thawed and treated with RNase A (cat#, 19101, Qiagen); 0.4-1ml was processed using the Qiagen QIAsymphony automated platform (Qiagen) using the QIAsymphony DSP DNA Midi Kit (cat# 937255, Qiagen). +DNA and RNA quantity and quality was assessed by PerkinElmer DropletQuant UV-VIS spectrophotometer (PerkinElmer) and an Agilent 4200 TapeStation (Agilent, USA) for RIN and DIN (RNA Integrity Number and DNA Integrity Number, respectively). +The NantHealth Sequencing Center, BGI at CHOP, or the Genomic Clinical Core at Sidra Medical and Research Center performed library preparation and sequencing. +BGI at CHOP and Sidra Medical and Research Center used in house, center-specific workflows for sample preparation. +At NantHealth Sequencing Center, DNA sequencing libraries were prepared for tumor and matched-normal DNA using the KAPA HyperPrep kit (cat# 08098107702, Roche), and tumor RNA-Seq libraries were prepared using KAPA Stranded RNA-Seq with RiboErase kit (cat# 07962304001, Roche).
+NantHealth and Sidra performed 2x150 bp WGS on paired tumor (~60X) and constitutive DNA (~30X) samples on an Illumina X/400. +BGI at CHOP performed 2x100 bp WGS sequenced at 60X depth for both tumor and normal samples. +NantHealth performed ribosomal-depleted whole transcriptome stranded RNA-Seq to an average depth of 200M. +BGI at CHOP performed poly-A or ribosomal-depleted whole transcriptome stranded RNA-Seq to an average depth of 100M. +The Translational Genomic Research Institute (TGEN; Phoenix, AZ) performed paired tumor (~200X) and constitutive whole exome sequencing (WXS) or targeted DNA panel (panel) and poly-A selected RNA-Seq (~200M reads) for PNOC tumor samples. +The panel tumor sample was sequenced to 470X, and the normal panel sample was sequenced to 308X. +PNOC 2x100 bp WXS and RNA-Seq libraries were sequenced on an Illumina HiSeq 2500.
+We used BWA-MEM
78 to align paired-end DNA-seq reads to the version 38 patch release 12 of the Homo sapiens genome reference, obtained as a FASTA file from UCSC (see Key Resources Table).
+Next, we used the Broad Institute’s Best Practices79 to process Binary Alignment/Map files (BAMs) in preparation for variant discovery.
+We marked duplicates using SAMBLASTER
80, and we merged and sorted BAMs using Sambamba
81
+We used the BaseRecalibrator
submodule of the Broad’s Genome Analysis Tool Kit GATK
82 to process BAM files.
+Lastly, for normal/germline input, we used the GATK
HaplotypeCaller
83 submodule on the recalibrated BAM to generate a genomic variant call format (GVCF) file.
+This file is used as the basis for germline calling, described in the SNV calling for B-allele Frequency (BAF) generation section.
We obtained references from the Broad Genome References on AWS bucket with a general description of references at https://s3.amazonaws.com/broad-references/broad-references-readme.html.
+To confirm sample matches and remove mis-matched samples from the dataset, we performed NGSCheckMate
84 on matched tumor/normal CRAM files.
+Briefly, we processed CRAMs using BCFtools
to filter and call 20k common single nucleotide polymorphisms (SNPs) using default parameters.
+We used the resulting VCFs to run NGSCheckMate
.
+Per NGSCheckMate
author recommendations, we used <= 0.61 as a correlation coefficient cutoff at sequencing depths > 10 to predict mis-matched samples.
+We determined RNA-Seq read strandedness by running the infer_experiment.py
script from RNA-SeQC
85 on the first 200k mapped reads.
+We removed any samples whose calculated strandedness did not match strandedness information provided by the sequencing center.
+We required that at least 60% of RNA-Seq reads mapped to the human reference for samples to be included in analysis.
+During OpenPBTA analysis, we identified some samples which were mis-identified or potentially swapped.
+Through collaborative analyses and pathology review, these samples were removed from our data releases and from the Kids First portal.
+Sample removal and associated justifications were documented in the OpenPBTA data release notes.
We performed germline haplotype calls using the GATK
Joint Genotyping Workflow on individual GVCFs from the normal sample alignment workflow.
+Using only SNPs, we applied the GATK
generic hard filter suggestions to the VCF, with an additional requirement of 10 reads minimum depth per SNP.
+We used the filtered VCF as input to Control-FREEC
and CNVkit
(below) to generate B-allele frequency (BAF) files.
+This single-sample workflow is available in the D3b GitHub repository.
+References can be obtained from the Broad Genome References on AWS bucket, and a general description of references can be found at https://s3.amazonaws.com/broad-references/broad-references-readme.html.
For patients with hypermutant samples, we first added population frequency of germline variants using ANNOVAR
86 and pathogenicity scoring from ClinVar87 using SnpSift
88.
+We then filtered for variants with read depth >= 15, variant allele fraction >= 0.20, and which were observed at < 0.1% allele frequency across each population in the Genome Aggregation Database (see Key Resources Table).
+Finally, we retained variants in genes included in the KEGG MMR gene set (see Key Resources Table), POLE, and/or TP53 which were ClinVar-annotated as pathogenic (P) or likely pathogenic (LP) with review status of >= 2 stars.
+All P/LP variants were manually reviewed by an interdisciplinary team of scientists, clinicians, and genetic counselors.
+This workflow is available in the D3b GitHub repository.
We used four variant callers to call SNVs and indels from paired tumor/normal samples with Targeted Panel, WXS, and/or WGS data: Strelka2
89, Mutect2
90, Lancet
91, and VarDictJava
92.
+VarDictJava
-only calls were not retained since ~ 39M calls with low VAF were uniquely called and may be potential false positives. (~1.2M calls were called by Mutect2
, Strelka2
, and Lancet
and included consensus CNV calling as described below.)
+We used only Strelka2
, Mutect2
and Lancet
to analyze WXS samples from TCGA.
+TCGA samples were captured using various WXS target capture kits and we downloaded the BED files from the GDC portal
.
+The manufacturers provided the input interval BED files for both panel and WXS data for PBTA samples.
+We padded all panel and WXS BED files were by 100 bp on each side for Strelka2
, Mutect2
, and VarDictJava
runs and by 400 bp for the Lancet
run.
+For WGS calling, we utilized the non-padded BROAD Institute interval calling list wgs_calling_regions.hg38.interval_list
, comprised of the full genome minus N bases, unless otherwise noted below.
+We ran Strelka2
89 using default parameters for canonical chromosomes (chr1-22, X,Y,M), as recommended by the authors, and we filtered the final Strelka2
VCF for PASS variants.
+We ran Mutect2
from GATK
according to Broad best practices outlined from their Workflow Description Language (WDL), and we filtered the final Mutect2
VCF for PASS variants.
+To manage memory issues, we ran VarDictJava
92 using 20 Kb interval chunks of the input BED, padded by 100 bp on each side, such that if an indel occurred in between intervals, it would be captured.
+Parameters and filtering followed BCBIO standards except that variants with a variant allele frequency (VAF) >= 0.05 (instead of >= 0.10) were retained.
+The 0.05 VAF increased the true positive rate for indels and decreased the false positive rate for SNVs when using VarDictJava
in consensus calling.
+We filtered the final VarDictJava
VCF for PASS variants with TYPE=StronglySomatic
.
+We ran Lancet
using default parameters, except for those noted below.
+For input intervals to Lancet
WGS, we created a reference BED from only the UTR, exome, and start/stop codon features of the GENCODE 31 reference, augmented as recommended with PASS variant calls from Strelka2
and Mutect2
.
+We then padded these intervals by 300 bp on each side during Lancet
variant calling.
+Per recommendations for WGS samples, we augmented the Lancet input intervals described above with PASS variant calls from Strelka2
and Mutect2
as validation93.
We normalized INDELs with bcftools norm
on all PASS VCFs using the kfdrc_annot_vcf_sub_wf.cwl
subworkflow, release v3 (See Table S5).
+The Ensembl Variant Effect Predictor (VEP
)94, reference release 93, was used to annotate variants and bcftools was used to add population allele frequency (AF) from gnomAD95.
+We annotated SNV and INDEL hotspots from v2 of Memorial Sloan Kettering Cancer Center’s (MSKCC) database (See Key Resources Table) as well as the TERT promoter mutations C228T and C250T96.
+We annotated SNVs by matching amino acid position (Protein_position
column in MAF file) with SNVs in the MSKCC database, we matched splice sites to HGVSp_Short
values in the MSKCC database, and we matched INDELs based on amino acid present within the range of INDEL hotspots values in the MSKCC database.
+We removed non-hotspot annotated variants with a normal depth less than or equal to 7 and/or gnomAD allele frequency (AF) greater than 0.001 as potential germline variants.
+We matched TERT promoter mutations using hg38 coordinates as indicated in ref.96: C228T occurs at 5:1295113 is annotated as existing variant s1242535815
, COSM1716563
, or COSM1716558
, and is 66 bp away from the TSS; C250T occurs at Chr5:1295135, is annotated as existing variant COSM1716559
, and is 88 bp away from the TSS.
+We retained variants annotated as PASS
or HotSpotAllele=1
in the final set, and we created MAFs using MSKCC’s vcf2maf
tool.
We retained all variant calls from Strelka2
, Mutect2
, or Lancet
that overlapped with an SNV or INDEL hotspot in a hotspot-specific MAF file, which we then used for select analyses as described below.
Our SNV calling process led to separate sets of predicted mutations for each caller.
+We considered mutations to describe the same change if they were identical for the following MAF fields: Chromosome
, Start_Position
, Reference_Allele
, Allele
, and Tumor_Sample_Barcode
.
+Strelka2
does not call multinucleotide variants (MNV), but instead calls each component SNV as a separate mutation, so we separated MNV calls from Mutect2
and Lancet
into consecutive SNVs before comparing them to Strelka2
calls.
+We examined VAFs produced by each caller and compared their overlap with each other (Figure S2).
+VarDictJava
calls included many variants that were not identified by other callers (Figure S2C), while the other callers produced results that were relatively consistent with one another.
+Many of these VarDictJava
-specific calls were variants with low allele frequency (Figure S2B).
+We therefore derived consensus mutation calls as those shared among the other three callers (Strelka2
, Mutect2
, and Lancet
), and we did not further consider VarDictJava
calls due to concerns it called a large number of false positives.
+This decision had minimal impact on results because VarDictJava
also identified nearly every mutation that the other three callers identified, in addition to many unique mutations.
We used Control-FREEC
97,98 and CNVkit
99 for copy number variant calls.
+For both algorithms, the germline_sex_estimate
(described below) was used as input for sample sex and germline variant calls (above) were used as input for BAF estimation.
+Control-FREEC
was run on human genome reference hg38 using the optional parameters of a 0.05 coefficient of variation, ploidy choice of 2-4, and BAF adjustment for tumor-normal pairs.
+Theta2
100 used VarDictJava
germline and somatic calls, filtered on PASS and strongly somatic, to infer tumor purity.
+Theta2
purity was added as an optional parameter to CNVkit
to adjust copy number calls.
+CNVkit
was run on human genome reference hg38 using the optional parameters of Theta2 purity and BAF adjustment for tumor-normal pairs.
+We used GISTIC
101 on the CNVkit
and the consensus CNV segmentation files to generate gene-level copy number abundance (Log R Ratio) as well as chromosomal arm copy number alterations using the parameters specified in the (run-gistic
analysis module in the OpenPBTA Analysis repository).
For each caller and sample, we called CNVs based on consensus among Control-FREEC
97,98, CNVkit
99, and Manta
102.
+We specifically included CNVs called significant by Control-FREEC
(p-value < 0.01) and Manta
calls that passed all filters in consensus calling.
+We removed sample and consensus caller files with more than 2,500 CNVs because we expected these to be noisy and derive poor quality samples based on cutoffs used in GISTIC
101.
+For each sample, we included the regions in the final consensus set: 1) regions with reciprocal overlap of 50% or more between at least two of the callers; 2) smaller CNV regions in which more than 90% of regions are covered by another caller.
+We did not include any copy number alteration called by a single algorithm in the consensus file.
+We defined copy number as NA
for any regions that had a neutral call for the samples included in the consensus file.
+We merged CNV regions within 10,000 bp of each other with the same direction of gain or loss into single region.
+We filtered out any CNVs that overlapped 50% or more with immunoglobulin, telomeric, centromeric, segment duplicated regions, or that were shorter than 3000 bp.
We used Manta
102 for structural variant (SV) calls, and we limited to regions used in Strelka2
.
+The hg38 reference for SV calling used was limited to canonical chromosome regions.
+We used AnnotSV
103 to annotate Manta
output.
+All associated workflows are available in the workflows GitHub repository.
We used STAR
104 to align paired-end RNA-seq reads, and we used the associated alignment for all subsequent RNA analysis.
+We used Ensembl GENCODE 27 “Comprehensive gene annotation” (see Key Resources Table) as a reference.
+We used RSEM
105 for both FPKM and TPM transcript- and gene-level quantification.
To enable downstream analyses, we next identified gene symbols that map to multiple Ensembl gene identifiers (in GENCODE v27, 212 gene symbols map to 1866 Ensembl gene identifiers), known as multi-mapped gene symbols, and ensured unique mappings (collapse-rnaseq
analysis module in the OpenPBTA Analysis repository).
+To this end, we first removed genes with no expression from the RSEM
abundance data by requiring an FPKM > 0 in at least 1 sample across the PBTA cohort.
+We computed the mean FPKM across all samples per gene.
+For each multi-mapped gene symbol, we chose the Ensembl identifier corresponding to the maximum mean FPKM, using the assumption that the gene identifier with the highest expression best represented the expression of the gene.
+After collapsing gene identifiers, 46,400 uniquely-expressed genes remained in the poly-A dataset, and 53,011 uniquely-expressed genes remained in the stranded dataset.
We set up Arriba
106 and STAR-Fusion
107 fusion detection tools using CWL on CAVATICA.
+For both of these tools, we used aligned BAM and chimeric SAM files from STAR
as inputs and GRCh38_gencode_v27
GTF for gene annotation.
+We ran STAR-Fusion
with default parameters and annotated all fusion calls with the GRCh38_v27_CTAT_lib_Feb092018.plug-n-play.tar.gz
file from the STAR-Fusion
release.
+For Arriba
, we used a blacklist file blacklist_hg38_GRCh38_2018-11-04.tsv.gz
from the Arriba
release to remove recurrent fusion artifacts and transcripts present in healthy tissue.
+We provided Arriba
with strandedness information for stranded samples, or we set it to auto-detection for poly-A samples.
+We used FusionAnnotator
on Arriba
fusion calls to harmonize annotations with those of STAR-Fusion
.
+The RNA expression and fusion workflows can be found in the D3b GitHub repository.
+The FusionAnnotator
workflow we used for this analysis can be found in the D3b GitHub repository.
All p-values are two-sided unless otherwise stated. +Z-scores were calculated using the formula \(z=(x –\mu)/\sigma\) where \(x\) is the value of interest, \(\mu\) is the mean, and \(\sigma\) is the standard deviation.
+tumor-purity-exploration
module)Estimating tumor fraction from RNA directly is challenging because most assume tumor cells comprise all non-immune cells108, which is not a valid assumption for many diagnoses in the PBTA cohort.
+We therefore used Theta2 (as described in the “Somatic Copy Number Variant Calling section” Methods section) to infer tumor purity from WGS samples, further assuming that co-extracted RNA and DNA samples had the same tumor purity.
+We then created a set of stranded RNA-Seq data thresholded by median tumor purity of the cancer group to rerun selected transcriptomic analyses: telomerase-activity-prediction
, tp53_nf1_score
, transcriptomic-dimension-reduction
, immune-deconv
, and gene-set-enrichment-analysis
.
+Note that these thresholded analyses, which only considered stranded RNA samples that also had co-extracted DNA, were performed in their respective OpenPBTA analyses modules (not within tumor-purity-exploration
).
interaction-plots
analysis module)Using the consensus SNV calls, we identified genes that were recurrently mutated in the OpenPBTA cohort, including nonsynonymous mutations with a VAF > 5% among the set of independent samples.
+We used VEP
94 annotations, including “High” and “Moderate” consequence types as defined in the R package Maftools
109, to determine the set of nonsynonymous mutations.
+For each gene, we then tallied the number of samples that had at least one nonsynonymous mutation.
For genes that contained nonsynonymous mutations in multiple samples, we calculated pairwise mutation co-occurrence scores. +This score was defined as \(I(-\log_{10}(P))\) where \(I\) is 1 when the odds ratio is > 1 (indicating co-occurrence), and -1 when the odds ratio is < 1 (indicating mutual exclusivity), with \(P\) defined by Fisher’s Exact Test.
+focal-cn-file-preparation
analysis module)We added the ploidy inferred via Control-FREEC
to the consensus CNV segmentation file and used the ploidy and copy number values to define gain and loss values broadly at the chromosome level.
+We used bedtools coverage
110 to add cytoband status using the UCSC cytoband file111 (See Key Resources Table).
+The output status call fractions, which are values of the loss, gain, and callable fractions of each cytoband region, were used to define dominant status at the cytoband-level.
+We calculated the weighted means of each status call fraction using band length.
+We used the weighted means to define the dominant status at the chromosome arm-level.
A status was considered dominant if more than half of the region was callable and the status call fraction was greater than 0.9 for that region. +We adopted this 0.9 threshold to ensure that the dominant status fraction call was greater than the remaining status fraction calls in a region.
+We aimed to define focal copy number units to avoid calling adjacent genes in the same cytoband or arm as copy number losses or gains where it would be more appropriate to call the broader region a loss or gain.
+To determine the most focal units, we first considered the dominant status calls at the chromosome arm-level.
+If the chromosome arm dominant status was callable but not clearly defined as a gain or loss, we instead included the cytoband-level status call.
+Similarly, if a cytoband dominant status call was callable but not clearly defined as a gain or loss, we instead included gene-level status call.
+To obtain the gene-level data, we used the IRanges
package in R112 to find overlaps between the segments in the consensus CNV file and the exons in the GENCODE v27 annotation file (See Key Resources Table) .
+If the copy number value was 0, we set the status to “deep deletion”.
+For autosomes only, we set the status to “amplification” when the copy number value was greater than two times the ploidy value.
+We plotted genome-wide gains and losses in (Figure S3C) using the R package ComplexHeatmap
113.
chromosomal-instability
analysis module)We defined breakpoint density as the number of breaks per genome or exome per sample. +For Manta SV calls, we filtered to retain “PASS” variants and used breakpoints from the algorithm. +For consensus CNV calls, if |log2 ratio| > log2(1), we annotated the segment as a break. +We then calculated breakpoint density as:
+\[\textrm{breakpoint density} = \frac{\textrm{N breaks}}{\textrm{Size in Mb of }\textit{effectively surveyed} \textrm{ genome}}\]
+chromothripsis
analysis module)Considering only chromosomes 1-22 and X, we identified candidate chromothripsis regions in the set of independent tumor WGS samples with ShatterSeek114, using Manta SV calls that passed all filters and consensus CNV calls.
+We modified the consensus CNV data to fit ShatterSeek
input requirements as follows: we set CNV-neutral or excluded regions as the respective sample’s ploidy value from Control-FREEC
, and we then merged consecutive segments with the same copy number value.
+We classified candidate chromothripsis regions as high- or low-confidence using the statistical criteria described by the ShatterSeek
authors.
immune-deconv
analysis module)We used the R package immunedeconv
pubmed:31510660? with the method quanTIseq
115 to deconvolute various immune cell types in tumors using collapsed FPKM RNA-seq, with samples batched by library type and then combined.
+The quanTIseq
deconvolution method directly estimates absolute fractions of 10 immune cell types that represent inferred proportions of the cell types in the mixture.
+Therefore, we utilized quanTIseq
for inter-sample, intra-sample, and inter-histology score comparisons.
gene-set-enrichment-analysis
analysis module)We performed Gene Set Variation Analysis (GSVA) on collapsed, log2-transformed RSEM FPKM data for stranded RNA-Seq samples using the GSVA
Bioconductor package116.
+We specified the parameter mx.diff=TRUE
to obtain Gaussian-distributed scores for each of the MSigDB hallmark gene sets117.
+We compared GSVA scores among histology groups using ANOVA and subsequent Tukey tests; p-values were Bonferroni-corrected for multiple hypothesis testing.
+We plotted scores by cancer group using the ComplexHeatmap
R package (Figure 5B)113.
transcriptomic-dimension-reduction
analysis module)We applied Uniform Manifold Approximation and Projection (UMAP)118 to log2-transformed FPKM data for stranded RNA-Seq samples using the umap
R package (See Key Resources Table).
+We considered all stranded RNA-Seq samples for this analysis, but we removed genes whose FPKM sum across samples was less than 100.
+We set the UMAP number of neighbors parameter to 15.
fusion_filtering
analysis module)We performed artifact filtering and additional annotation on fusion calls to prioritize putative oncogenic fusions. +Briefly, we considered all in-frame and frameshift fusion calls with at least one junction read and at least one gene partner expressed (TPM > 1) to be true calls. +If a fusion call had a large number of spanning fragment reads compared to junction reads (spanning fragment minus junction read greater than ten), we removed these calls as potential false positives. +We prioritized a union of fusion calls as true calls if the fused genes were detected by both callers, the same fusion was recurrent within a broad histology grouping (> 2 samples), or the fusion was specific to the given broad histology. +If either 5’ or 3’ genes fused to more than five different genes within a sample, we removed these calls as potential false positives. +We annotated putative driver fusions and prioritized fusions based on partners containing known kinases, oncogenes, tumor suppressors, curated transcription factors119, COSMIC genes, and/or known TCGA fusions from curated references. +Based on pediatric cancer literature review, we added MYBL1120, SNCAIP121, FOXR2122, TTYH1123, and TERT124–127 to the oncogene list, and we added BCOR122 and QKI128 to the tumor suppressor gene list.
+oncoprint-landscape
analysis module)We used Maftools
109 to generate oncoprints depicting the frequencies of canonical somatic gene mutations, CNVs, and fusions for the top 20 genes mutated across primary tumors within broad histologies of the OpenPBTA dataset.
+We collated canonical genes from the literature for low-grade gliomas (LGGs)24, embryonal tumors25,27,28,129,130, high-grade gliomas (HGGs)14,30,31,131, and other tumors: ependymomas, craniopharyngiomas, neuronal-glial mixed tumors, histiocytic tumors, chordoma, meningioma, and choroid plexus tumors32,132–138,pubmed:28623522?,pubmed:12466115?.
mutational-signatures
analysis module)We obtained weights (i.e., exposures) for signature sets using the deconstructSigs
R package function whichSignatures()
139 from consensus SNVs with the BSgenome.Hsapiens.UCSC.hg38 annotations (see Key Resources Table).
+Specifically, we estimated signature weights across samples for eight signatures previously identified in the Signal reference set of signatures (“RefSig”) as associated with adult central nervous system (CNS) tumors40.
+These eight RefSig signatures are 1, 3, 8, 11, 18, 19, N6, and MMR2.
+Weights for signatures fall in the range zero to one inclusive.
+deconstructSigs
estimates the weights for each signature across samples and allows for a proportion of unassigned weights referred to as “Other” in the text.
+These results do not include signatures with small contributions; deconstructSigs
drops signature weights that are less than 6%139.
+We plotted mutational signatures for patients with hypermutant tumors (Figure 4E) using the R package ComplexHeatmap
113.
snv-callers
analysis module)We consider tumor mutation burden (TMB) to be the number of consensus SNVs per effectively surveyed base of the genome. +We considered base pairs to be effectively surveyed if they were in the intersection of the genomic ranges considered by the callers used to generate the consensus and where appropriate, regions of interest, such as coding sequences. +We calculated TMB as:
+\[\textrm{TMB} = \frac{\textrm{# of coding sequence SNVs}}{\textrm{Size in Mb of }\textit{effectively surveyed} \textrm{ genome} }\]
+We used the total number coding sequence consensus SNVs for the numerator and the size of the intersection of the regions considered by Strelka2
and Mutect2
with coding regions (CDS from GENCODE v27 annotation, see Key Resources Table) as the denominator.
Table S1 contains a README, along with sample technical, clinical, and additional metadata used for this study.
+We performed molecular subtyping on tumors in the OpenPBTA to the extent possible.
+The molecular_subtype
field in pbta-histologies.tsv
contains molecular subtypes for tumor types selected from pathology_diagnosis
and pathology_free_text_diagnosis
fields as described below, following World Health Organization 2016 classification criteria19.
+We further categorized broad tumor histologies into smaller groupings we denote “cancer groups.”
Medulloblastoma (MB) subtypes SHH, WNT, Group 3, and Group 4 were predicted using the consensus of two RNA expression classifiers: MedulloClassifier
21 and MM2S
22 on the RSEM FPKM data (molecular-subtyping-MB
analysis module).
+The 43 “true positive” subtypes were manually curated from pathology reports by two independent reviewers.
High-grade glioma (HGG) subtypes were derived (molecular-subtyping-HGG
analysis module) using the following criteria:
DMG, H3 K28
.DMG, H3 K28, BRAF V600E
.HGG, H3 G35
.HGG, IDH
.BRAF V600E
.HGG, H3 wildtype
.Embryonal tumors were included in non-MB and non-ATRT embryonal tumor subtyping (molecular-subtyping-embryonal
analysis module) if they met any of the following criteria:
MN1::PATZ1
since it is an entity separate of CNS HGNET-MN1 tumors140.Non-MB and non-ATRT embryonal tumors identified with the above criteria were further subtyped (molecular-subtyping-embryonal
analysis module) using the criteria below141,142,pubmed:30249036?,pubmed:26389418?.
ETMR, C19MC-altered
(Embryonal tumor with multilayer rosettes, chromosome 19 miRNA cluster altered)123,143.ETMR, NOS
(Embryonal tumor with multilayer rosettes, not otherwise specified)123,143.CNS HGNET-MN1
[Central nervous system (CNS) high-grade neuroepithelial tumor with MN1 alteration].CNS HGNET-BCOR
(CNS high-grade neuroepithelial tumor with BCOR alteration).CNS NB-FOXR2
(CNS neuroblastoma with FOXR2 activation).CNS EFT-CIC
(CNS Ewing sarcoma family tumor with CIC alteration)122CNS Embryonal, NOS
(CNS Embryonal tumor, not otherwise specified).Neurocytoma subtypes central neurocytoma (CNC) and extraventricular neurocytoma (EVN) were assigned (molecular-subtyping-neurocytoma
analysis module) based on the primary site of the tumor145.
+If the tumor’s primary site was “ventricles,” we assigned the subtype as CNC; otherwise, we assigned the subtype as EVN.
Craniopharyngiomas (CRANIO) were subtyped (molecular-subtyping-CRANIO
analysis module) into adamantinomatous (CRANIO, ADAM
), papillary (CRANIO, PAP
) or undetermined (CRANIO, To be classified
) based on the following criteria146,147:
CRANIO, PAP
.CRANIO, ADAM
.CRANIO, To be classified
.A molecular subtype of EWS
was assigned to any tumor with a EWSR1 fusion or with a pathology_diagnosis
of Ewings Sarcoma
(molecular-subtyping-EWS
analysis module).
LGG or glialneuronal tumors (GNT) were subtyped (molecular-subtyping-LGAT
analysis module) based on SNV, fusion, and CNV status based on20 and as described below.
LGG, NF1-somatic
.LGG, NF1-germline
.LGG, IDH
.LGG, H3
.LGG, BRAF V600E
.KIAA1549::BRAF
fusion, it was subtyped as LGG, KIAA1549::BRAF
.KIAA1549::BRAF
, it was subtyped as LGG, other MAPK
.LGG, RTK
.LGG, FGFR
.LGG, MYB/MYBL1
.LGG, CDKN2A/B
.For LGG tumors that did not have any of the above molecular alterations, if both RNA and DNA samples were available, it was subtyped as LGG, wildtype
.
+Otherwise, if either RNA or DNA sample was unavailable, it was subtyped as LGG, To be classified
.
If pathology diagnosis was Subependymal Giant Cell Astrocytoma (SEGA)
, the LGG
portion of molecular subtype was recoded to SEGA
.
Lastly, for all LGG- and GNT- subtyped samples, if the tumors were glialneuronal in origin, based on pathology_free_text_diagnosis
entries of desmoplastic infantile
,desmoplastic infantile ganglioglioma
, desmoplastic infantile astrocytoma
or glioneuronal
, each was recoded as follows:
+If pathology diagnosis is Low-grade glioma/astrocytoma (WHO grade I/II)
or Ganglioglioma
, the LGG
portion of the molecular subtype was recoded to GNT
.
Ependymomas (EPN) were subtyped (molecular-subtyping-EPN
analysis module) into EPN, ST RELA
, EPN, ST YAP1
, EPN, PF A
and EPN, PF B
based on evidence for these molecular subgroups as described in Pajtler et al.132.
+Briefly, fusion, CNV and gene expression data were used to subtype EPN as follows:
RELA
as fusion partner, e.g., C11orf95::RELA
, LTBP3::RELA
, was subtyped as EPN, ST RELA
.YAP1
as fusion partner, such as C11orf95::YAP1
, YAP1::MAMLD1
and YAP1::FAM118B
, was subtyped as EPN, ST YAP1
.EPN, PF A
:EPN, PF B
:Any tumor with the above molecular characteristics would be exclusively subtyped to the designated group.
+For all other remaining EPN tumors without above molecular characteristics, they would be subtyped to EPN, ST RELA
and EPN, ST YAP1
in a non-exclusive way (e.g., a tumor could have both EPN, ST RELA
and EPN, ST YAP1
subtypes) if any of the following alterations were present.
EPN, ST RELA
:PTEN::TAS2R1
fusionEPN, ST YAP1
:C11orf95::MAML2
fusionAfter all relevant tumor samples were subtyped by the above molecular subtyping modules, the results from these modules, along with other clinical information (such as pathology diagnosis free text), were compiled in the molecular-subtyping-pathology
module and integrated into the OpenPBTA data in the molecular-subtyping-integrate
module.
tp53_nf1_score
analysis module)We annotated TP53 altered HGG samples as either TP53 lost
or TP53 activated
and integrated this within the molecular subtype.
+To this end, we applied a TP53 inactivation classifier originally trained on TCGA pan-cancer data42 to the matched RNA expression data, with samples batched by library type.
+Along with the TP53 classifier scores, we collectively used consensus SNV and CNV, SV, and reference databases that list TP53 hotspot mutations148,149 and functional domains150 to determine TP53 alteration status for each sample.
+We adopted the following rules for calling either TP53 lost
or TP53 activated
:
If a sample had either of the two well-characterized TP53 gain-of-function mutations, p.R273C or p.R248W43, we assigned TP53 activated
status.
Samples were annotated as TP53 lost
if they contained i) a TP53 hotspot mutation as defined by IARC TP53 database or the MSKCC cancer hotspots database148,149 (see also, Key Resources Table), ii) two TP53 alterations, including SNV, CNV or SV, indicative of probable bi-allelic alterations; iii) one TP53 somatic alteration, including SNV, CNV, or SV or a germline TP53 mutation indicated by the diagnosis of Li-Fraumeni syndrome (LFS)151, or iv) one germline TP53 mutation indicated by LFS and the TP53 classifier score for matched RNA-Seq was greater than 0.5.
Participant metadata included a reported gender.
+We used WGS germline data, in concert with the reported gender, to predict participant genetic sex so that we could identify sexually dimorphic outcomes.
+This analysis may also indicate samples that may have been contaminated.
+We used the idxstats
utility from SAMtools
pubmed:19505943? to calculate read lengths, the number of mapped reads, and the corresponding chromosomal location for reads to the X and Y chromosomes.
+We used the fraction of total normalized X and Y chromosome reads that were attributed to the Y chromosome as a summary statistic.
+We manually reviewed this statistic in the context of reported gender and determined that a threshold of less than 0.2 clearly delineated female samples.
+We marked fractions greater than 0.4 as predicted males, and we marked samples with values in the inclusive range 0.2-0.4 as unknown.
+We performed this analysis through CWL on CAVATICA.
+We added resulting calls to the histologies file under the column header germline_sex_estimate
.
independent-samples
analysis module)Certain analyses required that we select only a single representative specimen for each individual. +In these cases, we identified a single specimen by prioritizing primary tumors and those with whole-genome sequencing available. +If this filtering still resulted in multiple specimens, we randomly selected a single specimen from the remaining set.
+telomerase-activity-prediction
analysis module)We predicted telomerase activity of tumor samples using the recently developed EXTEND
method44, with samples batched by library type.
+Briefly, EXTEND
estimates telomerase activity based on the expression of a 13-gene signature.
+We derived this signature by comparing telomerase-positive tumors and tumors with activated alternative lengthening of telomeres pathway, a group presumably negative of telomerase activity.
survival-analysis
analysis module)We calculated overall survival (OS) as days since initial diagnosis and performed several survival analyses on the OpenPBTA cohort using the survival
R package.
+We performed survival analysis for patients by HGG subtype using the Kaplan-Meier estimator152 and a log-rank test (Mantel-Cox test)pubmed:5910392? on the different HGG subtypes.
+Next, we used multivariate Cox (proportional hazards) regression analysis153 to model the following: a) tp53 scores + telomerase scores + extent of tumor resection + LGG group + HGG group
, in which tp53 scores
and telomerase scores
are numeric, extent of tumor resection
is categorical, and LGG group
and HGG group
are binary variables indicating whether the sample is in either broad histology grouping, b) tp53 scores + telomerase scores + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), and c) quantiseq cell type fractions + CD274 expression + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), in which quantiseq cell type fractions
and CD274 expression
are numeric.
REAGENT or RESOURCE | +SOURCE | +IDENTIFIER | +
---|---|---|
Critical commercial assays | ++ | + |
Recover Cell Culture Freezing media | +Gibco | +12648010 | +
Hank’s Balanced Salt Solution (HBSS) | +Gibco | +14175095 | +
Papain | +SciQuest | +LS003124 | +
Ovomucoid | +SciQuest | +542000 | +
DNase | +Roche | +10104159001 | +
100μm cell strainer | +Greiner Bio-One | +542000 | +
DMEM/F-12 medium | +Sigma | +D8062 | +
Fetal Bovine Serum (FBS) | +Hyclone | +SH30910.03 | +
GlutaMAX | +Gibco | +35050061 | +
Penicillin/Streptomycin-Amphotericin B | +Lonza | +17-745E | +
Normocin | +Invivogen | +ant-nr-2 | +
B-27 supplement minus vitamin A | +Gibco | +12587-010 | +
N-2 supplement | +Gibco | +17502001 | +
Epidermal growth factor | +Gibco | +PHG0311L | +
Basic fibroblast growth factor | +PeproTech | +100-18B | +
Heparin | +Sigma | +H3149 | +
DNA/RNA AllPrep Kit | +Qiagen | +80204 | +
TruSeq RNA Sample Prep Kit | +Illumina | +FC-122-1001 | +
KAPA Library Preparation Kit | +Roche | +KK8201 | +
AllPrep DNA/RNA/miRNA Universal kit | +Qiagen | +80224 | +
RNase A | +Qiagen | +19101 | +
QIAsymphony DSP DNA Midi Kit | +Qiagen | +937255 | +
KAPA HyperPrep kit | +Roche | +08098107702 | +
RiboErase kit | +Roche | +07962304001 | +
Raw and harmonized WGS, WXS, Panel, RNA-Seq | +KidsFirst Data Resource Center, this project | +76 | +
Merged summary files | +this project | +https://cavatica.sbgenomics.com/u/cavatica/openpbta | +
Merged summary files and downstream analyses | +this project | +https://github.com/AlexsLemonade/OpenPBTA-analysis/ | +
Processed data | +this project | +https://pedcbioportal.kidsfirstdrc.org/study/summary?id=openpbta | +
Experimental models: Cell lines | ++ | + |
CBTN pediatric brain tumor-derived cell lines | +13 | +See Table S1 for identifiers | +
Software and algorithms | ++ | + |
Data processing and analysis software | +Multiple | +See Table S5 for identifiers | +
OpenPBTA workflows repository | +this project | +154 | +
OpenPBTA analysis repository | +this project | +155 | +
OpenPBTA manuscript repository | +this project | ++ |
+ | + | + |
Other | ++ | + |
TCGA WXS dataset | +National Institutes of Health The Cancer Genome Atlas (TCGA) | +dbGAP phs000178.v11.p8 | +
Cancer hotspots | +MSKCC | +https://www.cancerhotspots.org/#/download (v2) | +
Reference genomes | +Broad | +https://s3.console.aws.amazon.com/s3/buckets/broad-references/hg38/v0/ | +
Reference genome hg38, patch release 12 | +UCSC | +http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/ | +
Human Cytoband file | +UCSC | +http://hgdownload.cse.ucsc.edu/goldenpath/hg38/database/cytoBand.txt.gz | +
CDS from GENCODE v27 annotation | +GENCODE | +https://www.gencodegenes.org/human/release_27.html | +
PFAM domains and locations | +UCSC | +http://hgdownload.soe.ucsc.edu/goldenPath/hg38/database/pfamDesc.txt.gz; https://pfam.xfam.org/family/PF07714 | +
BSgenome.Hsapiens.UCSC.hg38 annotations | +Bioconductor | +https://bioconductor.org/packages/release/data/annotation/html/BSgenome.Hsapiens.UCSC.hg38.html | +
gnomAD v2.1.1 (exome and genome) | +Genome Aggregation Database | +https://gnomad.broadinstitute.org/downloads#v2-liftover-variants | +
KEGG MMR gene set v7.5.1 | +BROAD Institute | +https://www.gsea-msigdb.org/gsea/msigdb/download_geneset.jsp?geneSetName=KEGG_MISMATCH_REPAIR | +
ClinVar Database (2022-05-07) | +NCBI | +https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/archive_2.0/2022/clinvar_20220507.vcf.gz | +
Table S1. Related to Figure 1. +Table of specimens and associated metadata, clinical data, and histological data utilized in the OpenPBTA project.
+Table S2. Related to Figures 2 and 3. +Excel file with four sheets, where the first three represent tables of TMB, eight CNS mutational signatures, and chromothripsis events per sample, respectively, and the fourth sheet shows summarized genomic alterations across cancer groups.
+Table S3. Related to Figures 4 and 5. +Excel file with three sheets representing tables of TP53 scores, telomerase EXTEND scores, and quanTIseq immune scores, respectively.
+Table S4. Related to Figures 4 and 5. +Excel file with six sheets representing the survival analyses performed for this manuscript. +See Star Methods for details.
+Table S5. Related to Figure 1. +Excel file with four sheets representing of all software and their respective versions used for the OpenPBTA project, including the R packages in the OpenPBTA Docker image, Python packages i the OpenPBTA Docker image, other command line tools in the OpenPBTA Docker image, and all software used in the OpenPBTA workflows, respectively. +Note that all software in the OpenPBTA Docker image was utilized within the analysis repository, but not all software was used for the final manuscript.
+The past and present members of the Children’s Brain Tumor Network who contributed to the generation of specimens and data are Adam C. Resnick, Alexa Plisiewicz, Allison M. Morgan, Allison P. Heath, Alyssa Paul, Amanda Saratsis, Amy Smith, Ana Aguilar, Ana Guerreiro Stücklin, Anastasia Arynchyna, Andrea Franson, Angela J. Waanders, Angela N. Viaene, Anita Nirenberg, Anna Maria Buccoliero, Anna Yaffe, Anny Shai, Anthony Bet, Antoinette Price, Arlene Luther, Ashley Plant, Augustine Eze, Bailey K. Farrow, Baoli Hu, Beth Frenkel, Bo Zhang, Bobby Moulder, Bonnie Cole, Brian M. Ennis, Brian R. Rood, Brittany Lebert, Carina A. Leonard, Carl Koschmann, Caroline Caudill, Caroline Drinkwater, Cassie N. Kline, Catherine Sullivan, Chanel Keoni, Chiara Caporalini, Christine Bobick-Butcher, Christopher Mason, Chunde Li, Claire Carter, Claudia MaduroCoronado, Clayton Wiley, Cynthia Wong, David E. Kram, David Haussler, David Kram, David Pisapia, David Ziegler, Denise Morinigo, Derek Hanson, Donald W. Parsons, Elizabeth Appert, Emily Drake, Emily Golbeck, Ena Agbodza, Eric H. Raabe, Eric M. Jackson, Erin Alexander, Esteban Uceda, Eugene Hwang, Fausto Rodriquez, Gabrielle S. Stone, Gary Kohanbash, Gavriella Silverman, George Rafidi, Gerald Grant, Gerri Trooskin, Gilad Evrony, Graham Keyes, Hagop Boyajian, Holly B. Lindsay, Holly C. Beale, Ian F. Pollack, James Johnston, James Palmer, Jane Minturn, Jared Pisapia, Jason E. Cain, Jason R. Fangusaro, Javad Nazarian, Jeanette Haugh, Jeff Stevens, Jeffrey P. Greenfield, Jeffrey Rubens, Jena V. Lilly, Jennifer L. Mason, Jessica B. Foster, Jim Olson, Jo Lynne Rokita, Joanna J. Phillips, Jonathan Waller, Josh Rubin, Judy E. Palma, Justin McCroskey, Justine Rizzo, Kaitlin Lehmann, Kamnaa Arya, Karlene Hall, Katherine Pehlivan, Kenneth Seidl, Kimberly Diamond, Kristen Harnett, Kristina A. Cole, Krutika S. Gaonkar, Lamiya Tauhid, Laura Prolo, Leah Holloway, Leslie Brosig, Lina Lopez, Lionel Chow, Madhuri Kambhampati, Mahdi Sarmady, Margaret Nevins, Mari Groves, Mariarita Santi-Vicini, Marilyn M. Li, Marion Mateos, Mateusz Koptyra, Matija Snuderl, Matthew Miller, Matthew Sklar, Matthew Wood, Meghan Connors, Melissa Williams, Meredith Egan, Michael Fisher, Michael Koldobskiy, Michelle Monje, Migdalia Martinez, Miguel A. Brown, Mike Prados, Miriam Bornhorst, Mirko Scagnet, Mohamed AbdelBaki, Monique Carrero-Tagle, Nadia Dahmane, Nalin Gupta, Nathan Young, Nicholas A. Vitanza, Nicholas Tassone, Nicholas Van Kuren, Nicolas Gerber, Nithin D. Adappa, Nitin Wadhwani, Noel Coleman, Obi Obayashi, Olena M. Vaske, Olivier Elemento, Oren Becher, Philbert Oliveros, Phillip B. Storm, Pichai Raman, Prajwal Rajappa, Rintaro Hashizume, Rishi R. Lulla, Robert Keating, Robert M. Lober, Ron Firestein, Sabine Mueller, Sameer Agnihotri, Samuel G. Winebrake, Samuel Rivero-Hinojosa, Sarah Diane Black, Sarah Leary, Schuyler Stoller, Shannon Robins, Sharon Gardner, Shelly Wang, Sherri Mayans, Sherry Tutson, Shida Zhu, Sofie R. Salama, Sonia Partap, Sonika Dahiya, Sriram Venneti, Stacie Stapleton, Stephani Campion, Stephanie Stefankiewicz, Stewart Goldman, Swetha Thambireddy, Tatiana S. Patton, Teresa Hidalgo, Theo Nicolaides, Thinh Q. Nguyen, Thomas W. McLean, Tiffany Walker, Toba Niazi, Tobey MacDonald, Valeria Lopez-Gil, Valerie Baubet, Whitney Rife, Xiao-Nan Li, Ximena Cuellar, Yiran Guo, Yuankun Zhu, and Zeinab Helil.
+The past and present members of the Pacific Pediatric Neuro-Oncology Consortium who contributed to the generation of specimens and data are Adam C. Resnick, Alicia Lenzen, Alyssa Reddy, Amar Gajjar, Ana Guerreiro Stucklin, Anat Epstein, Andrea Franson, Angela Waanders, Anne Bendel, Anu Banerjee, Ashley Margol, Ashley Plant, Brian Rood, Carl Koschmann, Carol Bruggers, Caroline Hastings, Cassie N. Kline, Christina Coleman Abadi, Christopher Tinkle, Corey Raffel, Dan Runco, Daniel Landi, Daphne Adele Haas-Kogan, David Ashley, David Ziegler, Derek Hanson, Dong Anh Khuong Quang, Duane Mitchell, Elias Sayour, Eric Jackson, Eric Raabe, Eugene Hwang, Fatema Malbari, Geoffrey McCowage, Girish Dhall, Gregory Friedman, Hideho Okada, Ibrahim Qaddoumi, Iris Fried, Jae Cho, Jane Minturn, Jason Blatt, Javad Nazarian, Jeffrey Rubens, Jena V. Lilly, Jennifer Elster, Jennifer L. Mason, Jessica Schulte, Jonathan Schoenfeld, Josh Rubin, Karen Gauvain, Karen Wright, Katharine Offer, Katie Metrock, Kellie Haworth, Ken Cohen, Kristina A. Cole, Lance Governale, Linda Stork, Lindsay Kilburn, Lissa Baird, Maggie Skrypek, Marcia Leonard, Margaret Shatara, Margot Lazow, Mariella Filbin, Maryam Fouladi, Matthew Miller, Megan Paul, Michael Fisher, Michael Koldobskiy, Michael Prados, Michal Yalon Oren, Mimi Bandopadhayay, Miriam Bornhorst, Mohamed AbdelBaki, Nalin Gupta, Nathan Robison, Nicholas Whipple, Nick Gottardo, Nicholas A. Vitanza, Nicolas Gerber, Patricia Robertson, Payal Jain, Peter Sun, Priya Chan, Richard S Lemons, Robert Wechsler-Reya, Roger Packer, Russ Geyer, Ryan Velasco, Sabine Mueller, Sahaja Acharya, Sam Cheshier, Sarah Leary, Scott Coven, Sebastian M. Waszak, Sharon Gardner, Sri Gururangan, Stewart Goldman, Susan Chi, Tab Cooney, Tatiana S. Patton, Theodore Nicolaides, and Tom Belle Davidson.
+This manuscript -(permalink) +(permalink) was automatically generated -from AlexsLemonade/OpenPBTA-manuscript@cfc4c7d +from AlexsLemonade/OpenPBTA-manuscript@b2753bb on April 6, 2023.
All original code was developed within the following repositories and is publicly available as follows. Primary data analyses can be found at https://github.com/d3b-center/OpenPBTA-workflows. @@ -2230,19 +2231,19 @@
We used BWA-MEM
77 to align paired-end DNA-seq reads to the version 38 patch release 12 of the Homo sapiens genome reference, obtained as a FASTA file from UCSC (see Key Resources Table).
-Next, we used the Broad Institute’s Best Practices78 to process Binary Alignment/Map files (BAMs) in preparation for variant discovery.
-We marked duplicates using SAMBLASTER
79, and we merged and sorted BAMs using Sambamba
80
-We used the BaseRecalibrator
submodule of the Broad’s Genome Analysis Tool Kit GATK
81 to process BAM files.
-Lastly, for normal/germline input, we used the GATK
HaplotypeCaller
82 submodule on the recalibrated BAM to generate a genomic variant call format (GVCF) file.
+
We used BWA-MEM
78 to align paired-end DNA-seq reads to the version 38 patch release 12 of the Homo sapiens genome reference, obtained as a FASTA file from UCSC (see Key Resources Table).
+Next, we used the Broad Institute’s Best Practices79 to process Binary Alignment/Map files (BAMs) in preparation for variant discovery.
+We marked duplicates using SAMBLASTER
80, and we merged and sorted BAMs using Sambamba
81
+We used the BaseRecalibrator
submodule of the Broad’s Genome Analysis Tool Kit GATK
82 to process BAM files.
+Lastly, for normal/germline input, we used the GATK
HaplotypeCaller
83 submodule on the recalibrated BAM to generate a genomic variant call format (GVCF) file.
This file is used as the basis for germline calling, described in the SNV calling for B-allele Frequency (BAF) generation section.
We obtained references from the Broad Genome References on AWS bucket with a general description of references at https://s3.amazonaws.com/broad-references/broad-references-readme.html.
To confirm sample matches and remove mis-matched samples from the dataset, we performed NGSCheckMate
83 on matched tumor/normal CRAM files.
+
To confirm sample matches and remove mis-matched samples from the dataset, we performed NGSCheckMate
84 on matched tumor/normal CRAM files.
Briefly, we processed CRAMs using BCFtools
to filter and call 20k common single nucleotide polymorphisms (SNPs) using default parameters.
We used the resulting VCFs to run NGSCheckMate
.
Per NGSCheckMate
author recommendations, we used <= 0.61 as a correlation coefficient cutoff at sequencing depths > 10 to predict mis-matched samples.
-We determined RNA-Seq read strandedness by running the infer_experiment.py
script from RNA-SeQC
84 on the first 200k mapped reads.
+We determined RNA-Seq read strandedness by running the infer_experiment.py
script from RNA-SeQC
85 on the first 200k mapped reads.
We removed any samples whose calculated strandedness did not match strandedness information provided by the sequencing center.
We required that at least 60% of RNA-Seq reads mapped to the human reference for samples to be included in analysis.
During OpenPBTA analysis, we identified some samples which were mis-identified or potentially swapped.
@@ -2256,37 +2257,37 @@
For patients with hypermutant samples, we first added population frequency of germline variants using ANNOVAR
85 and pathogenicity scoring from ClinVar86 using SnpSift
87.
+
For patients with hypermutant samples, we first added population frequency of germline variants using ANNOVAR
86 and pathogenicity scoring from ClinVar87 using SnpSift
88.
We then filtered for variants with read depth >= 15, variant allele fraction >= 0.20, and which were observed at < 0.1% allele frequency across each population in the Genome Aggregation Database (see Key Resources Table).
Finally, we retained variants in genes included in the KEGG MMR gene set (see Key Resources Table), POLE, and/or TP53 which were ClinVar-annotated as pathogenic (P) or likely pathogenic (LP) with review status of >= 2 stars.
All P/LP variants were manually reviewed by an interdisciplinary team of scientists, clinicians, and genetic counselors.
This workflow is available in the D3b GitHub repository.
We used four variant callers to call SNVs and indels from paired tumor/normal samples with Targeted Panel, WXS, and/or WGS data: Strelka2
88, Mutect2
89, Lancet
90, and VarDictJava
91.
+
We used four variant callers to call SNVs and indels from paired tumor/normal samples with Targeted Panel, WXS, and/or WGS data: Strelka2
89, Mutect2
90, Lancet
91, and VarDictJava
92.
VarDictJava
-only calls were not retained since ~ 39M calls with low VAF were uniquely called and may be potential false positives. (~1.2M calls were called by Mutect2
, Strelka2
, and Lancet
and included consensus CNV calling as described below.)
We used only Strelka2
, Mutect2
and Lancet
to analyze WXS samples from TCGA.
TCGA samples were captured using various WXS target capture kits and we downloaded the BED files from the GDC portal
.
The manufacturers provided the input interval BED files for both panel and WXS data for PBTA samples.
We padded all panel and WXS BED files were by 100 bp on each side for Strelka2
, Mutect2
, and VarDictJava
runs and by 400 bp for the Lancet
run.
For WGS calling, we utilized the non-padded BROAD Institute interval calling list wgs_calling_regions.hg38.interval_list
, comprised of the full genome minus N bases, unless otherwise noted below.
-We ran Strelka2
88 using default parameters for canonical chromosomes (chr1-22, X,Y,M), as recommended by the authors, and we filtered the final Strelka2
VCF for PASS variants.
+We ran Strelka2
89 using default parameters for canonical chromosomes (chr1-22, X,Y,M), as recommended by the authors, and we filtered the final Strelka2
VCF for PASS variants.
We ran Mutect2
from GATK
according to Broad best practices outlined from their Workflow Description Language (WDL), and we filtered the final Mutect2
VCF for PASS variants.
-To manage memory issues, we ran VarDictJava
91 using 20 Kb interval chunks of the input BED, padded by 100 bp on each side, such that if an indel occurred in between intervals, it would be captured.
+To manage memory issues, we ran VarDictJava
92 using 20 Kb interval chunks of the input BED, padded by 100 bp on each side, such that if an indel occurred in between intervals, it would be captured.
Parameters and filtering followed BCBIO standards except that variants with a variant allele frequency (VAF) >= 0.05 (instead of >= 0.10) were retained.
The 0.05 VAF increased the true positive rate for indels and decreased the false positive rate for SNVs when using VarDictJava
in consensus calling.
We filtered the final VarDictJava
VCF for PASS variants with TYPE=StronglySomatic
.
We ran Lancet
using default parameters, except for those noted below.
For input intervals to Lancet
WGS, we created a reference BED from only the UTR, exome, and start/stop codon features of the GENCODE 31 reference, augmented as recommended with PASS variant calls from Strelka2
and Mutect2
.
We then padded these intervals by 300 bp on each side during Lancet
variant calling.
-Per recommendations for WGS samples, we augmented the Lancet input intervals described above with PASS variant calls from Strelka2
and Mutect2
as validation92.
Strelka2
and Mutect2
as validation93.
We normalized INDELs with bcftools norm
on all PASS VCFs using the kfdrc_annot_vcf_sub_wf.cwl
subworkflow, release v3 (See Table S5).
-The Ensembl Variant Effect Predictor (VEP
)93, reference release 93, was used to annotate variants and bcftools was used to add population allele frequency (AF) from gnomAD94.
-We annotated SNV and INDEL hotspots from v2 of Memorial Sloan Kettering Cancer Center’s (MSKCC) database (See Key Resources Table) as well as the TERT promoter mutations C228T and C250T95.
+The Ensembl Variant Effect Predictor (VEP
)94, reference release 93, was used to annotate variants and bcftools was used to add population allele frequency (AF) from gnomAD95.
+We annotated SNV and INDEL hotspots from v2 of Memorial Sloan Kettering Cancer Center’s (MSKCC) database (See Key Resources Table) as well as the TERT promoter mutations C228T and C250T96.
We annotated SNVs by matching amino acid position (Protein_position
column in MAF file) with SNVs in the MSKCC database, we matched splice sites to HGVSp_Short
values in the MSKCC database, and we matched INDELs based on amino acid present within the range of INDEL hotspots values in the MSKCC database.
We removed non-hotspot annotated variants with a normal depth less than or equal to 7 and/or gnomAD allele frequency (AF) greater than 0.001 as potential germline variants.
-We matched TERT promoter mutations using hg38 coordinates as indicated in ref.95: C228T occurs at 5:1295113 is annotated as existing variant s1242535815
, COSM1716563
, or COSM1716558
, and is 66 bp away from the TSS; C250T occurs at Chr5:1295135, is annotated as existing variant COSM1716559
, and is 88 bp away from the TSS.
+We matched TERT promoter mutations using hg38 coordinates as indicated in ref.96: C228T occurs at 5:1295113 is annotated as existing variant s1242535815
, COSM1716563
, or COSM1716558
, and is 66 bp away from the TSS; C250T occurs at Chr5:1295135, is annotated as existing variant COSM1716559
, and is 88 bp away from the TSS.
We retained variants annotated as PASS
or HotSpotAllele=1
in the final set, and we created MAFs using MSKCC’s vcf2maf
tool.
We retained all variant calls from Strelka2
, Mutect2
, or Lancet
that overlapped with an SNV or INDEL hotspot in a hotspot-specific MAF file, which we then used for select analyses as described below.
Strelka2
, Mutect2
, and Lancet
), and we did not further consider VarDictJava
calls due to concerns it called a large number of false positives.
This decision had minimal impact on results because VarDictJava
also identified nearly every mutation that the other three callers identified, in addition to many unique mutations.
We used Control-FREEC
96,97 and CNVkit
98 for copy number variant calls.
+
We used Control-FREEC
97,98 and CNVkit
99 for copy number variant calls.
For both algorithms, the germline_sex_estimate
(described below) was used as input for sample sex and germline variant calls (above) were used as input for BAF estimation.
Control-FREEC
was run on human genome reference hg38 using the optional parameters of a 0.05 coefficient of variation, ploidy choice of 2-4, and BAF adjustment for tumor-normal pairs.
-Theta2
99 used VarDictJava
germline and somatic calls, filtered on PASS and strongly somatic, to infer tumor purity.
+Theta2
100 used VarDictJava
germline and somatic calls, filtered on PASS and strongly somatic, to infer tumor purity.
Theta2
purity was added as an optional parameter to CNVkit
to adjust copy number calls.
CNVkit
was run on human genome reference hg38 using the optional parameters of Theta2 purity and BAF adjustment for tumor-normal pairs.
-We used GISTIC
100 on the CNVkit
and the consensus CNV segmentation files to generate gene-level copy number abundance (Log R Ratio) as well as chromosomal arm copy number alterations using the parameters specified in the (run-gistic
analysis module in the OpenPBTA Analysis repository).
GISTIC
101 on the CNVkit
and the consensus CNV segmentation files to generate gene-level copy number abundance (Log R Ratio) as well as chromosomal arm copy number alterations using the parameters specified in the (run-gistic
analysis module in the OpenPBTA Analysis repository).
For each caller and sample, we called CNVs based on consensus among Control-FREEC
96,97, CNVkit
98, and Manta
101.
+
For each caller and sample, we called CNVs based on consensus among Control-FREEC
97,98, CNVkit
99, and Manta
102.
We specifically included CNVs called significant by Control-FREEC
(p-value < 0.01) and Manta
calls that passed all filters in consensus calling.
-We removed sample and consensus caller files with more than 2,500 CNVs because we expected these to be noisy and derive poor quality samples based on cutoffs used in GISTIC
100.
+We removed sample and consensus caller files with more than 2,500 CNVs because we expected these to be noisy and derive poor quality samples based on cutoffs used in GISTIC
101.
For each sample, we included the regions in the final consensus set: 1) regions with reciprocal overlap of 50% or more between at least two of the callers; 2) smaller CNV regions in which more than 90% of regions are covered by another caller.
We did not include any copy number alteration called by a single algorithm in the consensus file.
We defined copy number as NA
for any regions that had a neutral call for the samples included in the consensus file.
We merged CNV regions within 10,000 bp of each other with the same direction of gain or loss into single region.
We filtered out any CNVs that overlapped 50% or more with immunoglobulin, telomeric, centromeric, segment duplicated regions, or that were shorter than 3000 bp.
We used Manta
101 for structural variant (SV) calls, and we limited to regions used in Strelka2
.
+
We used Manta
102 for structural variant (SV) calls, and we limited to regions used in Strelka2
.
The hg38 reference for SV calling used was limited to canonical chromosome regions.
-We used AnnotSV
102 to annotate Manta
output.
+We used AnnotSV
103 to annotate Manta
output.
All associated workflows are available in the workflows GitHub repository.
We used STAR
103 to align paired-end RNA-seq reads, and we used the associated alignment for all subsequent RNA analysis.
+
We used STAR
104 to align paired-end RNA-seq reads, and we used the associated alignment for all subsequent RNA analysis.
We used Ensembl GENCODE 27 “Comprehensive gene annotation” (see Key Resources Table) as a reference.
-We used RSEM
104 for both FPKM and TPM transcript- and gene-level quantification.
RSEM
105 for both FPKM and TPM transcript- and gene-level quantification.
To enable downstream analyses, we next identified gene symbols that map to multiple Ensembl gene identifiers (in GENCODE v27, 212 gene symbols map to 1866 Ensembl gene identifiers), known as multi-mapped gene symbols, and ensured unique mappings (collapse-rnaseq
analysis module in the OpenPBTA Analysis repository).
To this end, we first removed genes with no expression from the RSEM
abundance data by requiring an FPKM > 0 in at least 1 sample across the PBTA cohort.
@@ -2333,7 +2334,7 @@
We set up Arriba
105 and STAR-Fusion
106 fusion detection tools using CWL on CAVATICA.
+
We set up Arriba
106 and STAR-Fusion
107 fusion detection tools using CWL on CAVATICA.
For both of these tools, we used aligned BAM and chimeric SAM files from STAR
as inputs and GRCh38_gencode_v27
GTF for gene annotation.
We ran STAR-Fusion
with default parameters and annotated all fusion calls with the GRCh38_v27_CTAT_lib_Feb092018.plug-n-play.tar.gz
file from the STAR-Fusion
release.
For Arriba
, we used a blacklist file blacklist_hg38_GRCh38_2018-11-04.tsv.gz
from the Arriba
release to remove recurrent fusion artifacts and transcripts present in healthy tissue.
@@ -2345,19 +2346,19 @@
All p-values are two-sided unless otherwise stated. Z-scores were calculated using the formula \(z=(x –\mu)/\sigma\) where \(x\) is the value of interest, \(\mu\) is the mean, and \(\sigma\) is the standard deviation.
tumor-purity-exploration
module)Estimating tumor fraction from RNA directly is challenging because most assume tumor cells comprise all non-immune cells107, which is not a valid assumption for many diagnoses in the PBTA cohort. +
Estimating tumor fraction from RNA directly is challenging because most assume tumor cells comprise all non-immune cells108, which is not a valid assumption for many diagnoses in the PBTA cohort.
We therefore used Theta2 (as described in the “Somatic Copy Number Variant Calling section” Methods section) to infer tumor purity from WGS samples, further assuming that co-extracted RNA and DNA samples had the same tumor purity.
We then created a set of stranded RNA-Seq data thresholded by median tumor purity of the cancer group to rerun selected transcriptomic analyses: telomerase-activity-prediction
, tp53_nf1_score
, transcriptomic-dimension-reduction
, immune-deconv
, and gene-set-enrichment-analysis
.
Note that these thresholded analyses, which only considered stranded RNA samples that also had co-extracted DNA, were performed in their respective OpenPBTA analyses modules (not within tumor-purity-exploration
).
interaction-plots
analysis module)Using the consensus SNV calls, we identified genes that were recurrently mutated in the OpenPBTA cohort, including nonsynonymous mutations with a VAF > 5% among the set of independent samples.
-We used VEP
93 annotations, including “High” and “Moderate” consequence types as defined in the R package Maftools
108, to determine the set of nonsynonymous mutations.
+We used VEP
94 annotations, including “High” and “Moderate” consequence types as defined in the R package Maftools
109, to determine the set of nonsynonymous mutations.
For each gene, we then tallied the number of samples that had at least one nonsynonymous mutation.
For genes that contained nonsynonymous mutations in multiple samples, we calculated pairwise mutation co-occurrence scores. This score was defined as \(I(-\log_{10}(P))\) where \(I\) is 1 when the odds ratio is > 1 (indicating co-occurrence), and -1 when the odds ratio is < 1 (indicating mutual exclusivity), with \(P\) defined by Fisher’s Exact Test.
focal-cn-file-preparation
analysis module)We added the ploidy inferred via Control-FREEC
to the consensus CNV segmentation file and used the ploidy and copy number values to define gain and loss values broadly at the chromosome level.
-We used bedtools coverage
109 to add cytoband status using the UCSC cytoband file110 (See Key Resources Table).
+We used bedtools coverage
110 to add cytoband status using the UCSC cytoband file111 (See Key Resources Table).
The output status call fractions, which are values of the loss, gain, and callable fractions of each cytoband region, were used to define dominant status at the cytoband-level.
We calculated the weighted means of each status call fraction using band length.
We used the weighted means to define the dominant status at the chromosome arm-level.
IRanges
package in R111 to find overlaps between the segments in the consensus CNV file and the exons in the GENCODE v27 annotation file (See Key Resources Table) .
+To obtain the gene-level data, we used the IRanges
package in R112 to find overlaps between the segments in the consensus CNV file and the exons in the GENCODE v27 annotation file (See Key Resources Table) .
If the copy number value was 0, we set the status to “deep deletion”.
For autosomes only, we set the status to “amplification” when the copy number value was greater than two times the ploidy value.
-We plotted genome-wide gains and losses in (Figure S3C) using the R package ComplexHeatmap
112.
+We plotted genome-wide gains and losses in (Figure S3C) using the R package ComplexHeatmap
113.
chromosomal-instability
analysis module)We defined breakpoint density as the number of breaks per genome or exome per sample. For Manta SV calls, we filtered to retain “PASS” variants and used breakpoints from the algorithm. @@ -2378,20 +2379,20 @@
chromothripsis
analysis module)Considering only chromosomes 1-22 and X, we identified candidate chromothripsis regions in the set of independent tumor WGS samples with ShatterSeek113, using Manta SV calls that passed all filters and consensus CNV calls. +
Considering only chromosomes 1-22 and X, we identified candidate chromothripsis regions in the set of independent tumor WGS samples with ShatterSeek114, using Manta SV calls that passed all filters and consensus CNV calls.
We modified the consensus CNV data to fit ShatterSeek
input requirements as follows: we set CNV-neutral or excluded regions as the respective sample’s ploidy value from Control-FREEC
, and we then merged consecutive segments with the same copy number value.
We classified candidate chromothripsis regions as high- or low-confidence using the statistical criteria described by the ShatterSeek
authors.
immune-deconv
analysis module)We used the R package immunedeconv
pubmed:31510660? with the method quanTIseq
114 to deconvolute various immune cell types in tumors using collapsed FPKM RNA-seq, with samples batched by library type and then combined.
+
We used the R package immunedeconv
pubmed:31510660? with the method quanTIseq
115 to deconvolute various immune cell types in tumors using collapsed FPKM RNA-seq, with samples batched by library type and then combined.
The quanTIseq
deconvolution method directly estimates absolute fractions of 10 immune cell types that represent inferred proportions of the cell types in the mixture.
Therefore, we utilized quanTIseq
for inter-sample, intra-sample, and inter-histology score comparisons.
gene-set-enrichment-analysis
analysis module)We performed Gene Set Variation Analysis (GSVA) on collapsed, log2-transformed RSEM FPKM data for stranded RNA-Seq samples using the GSVA
Bioconductor package115.
-We specified the parameter mx.diff=TRUE
to obtain Gaussian-distributed scores for each of the MSigDB hallmark gene sets116.
+
We performed Gene Set Variation Analysis (GSVA) on collapsed, log2-transformed RSEM FPKM data for stranded RNA-Seq samples using the GSVA
Bioconductor package116.
+We specified the parameter mx.diff=TRUE
to obtain Gaussian-distributed scores for each of the MSigDB hallmark gene sets117.
We compared GSVA scores among histology groups using ANOVA and subsequent Tukey tests; p-values were Bonferroni-corrected for multiple hypothesis testing.
-We plotted scores by cancer group using the ComplexHeatmap
R package (Figure 5B)112.
ComplexHeatmap
R package (Figure 5B)113.
transcriptomic-dimension-reduction
analysis module)We applied Uniform Manifold Approximation and Projection (UMAP)117 to log2-transformed FPKM data for stranded RNA-Seq samples using the umap
R package (See Key Resources Table).
+
We applied Uniform Manifold Approximation and Projection (UMAP)118 to log2-transformed FPKM data for stranded RNA-Seq samples using the umap
R package (See Key Resources Table).
We considered all stranded RNA-Seq samples for this analysis, but we removed genes whose FPKM sum across samples was less than 100.
We set the UMAP number of neighbors parameter to 15.
fusion_filtering
analysis module)oncoprint-landscape
analysis module)We used Maftools
108 to generate oncoprints depicting the frequencies of canonical somatic gene mutations, CNVs, and fusions for the top 20 genes mutated across primary tumors within broad histologies of the OpenPBTA dataset.
-We collated canonical genes from the literature for low-grade gliomas (LGGs)24, embryonal tumors25,27,28,128,129, high-grade gliomas (HGGs)14,30,31,130, and other tumors: ependymomas, craniopharyngiomas, neuronal-glial mixed tumors, histiocytic tumors, chordoma, meningioma, and choroid plexus tumors32,131–137,pubmed:28623522?,pubmed:12466115?.
We used Maftools
109 to generate oncoprints depicting the frequencies of canonical somatic gene mutations, CNVs, and fusions for the top 20 genes mutated across primary tumors within broad histologies of the OpenPBTA dataset.
+We collated canonical genes from the literature for low-grade gliomas (LGGs)24, embryonal tumors25,27,28,129,130, high-grade gliomas (HGGs)14,30,31,131, and other tumors: ependymomas, craniopharyngiomas, neuronal-glial mixed tumors, histiocytic tumors, chordoma, meningioma, and choroid plexus tumors32,132–138,pubmed:28623522?,pubmed:12466115?.
mutational-signatures
analysis module)We obtained weights (i.e., exposures) for signature sets using the deconstructSigs
R package function whichSignatures()
138 from consensus SNVs with the BSgenome.Hsapiens.UCSC.hg38 annotations (see Key Resources Table).
+
We obtained weights (i.e., exposures) for signature sets using the deconstructSigs
R package function whichSignatures()
139 from consensus SNVs with the BSgenome.Hsapiens.UCSC.hg38 annotations (see Key Resources Table).
Specifically, we estimated signature weights across samples for eight signatures previously identified in the Signal reference set of signatures (“RefSig”) as associated with adult central nervous system (CNS) tumors40.
These eight RefSig signatures are 1, 3, 8, 11, 18, 19, N6, and MMR2.
Weights for signatures fall in the range zero to one inclusive.
deconstructSigs
estimates the weights for each signature across samples and allows for a proportion of unassigned weights referred to as “Other” in the text.
-These results do not include signatures with small contributions; deconstructSigs
drops signature weights that are less than 6%138.
-We plotted mutational signatures for patients with hypermutant tumors (Figure 4E) using the R package ComplexHeatmap
112.
deconstructSigs
drops signature weights that are less than 6%139.
+We plotted mutational signatures for patients with hypermutant tumors (Figure 4E) using the R package ComplexHeatmap
113.
snv-callers
analysis module)We consider tumor mutation burden (TMB) to be the number of consensus SNVs per effectively surveyed base of the genome. We considered base pairs to be effectively surveyed if they were in the intersection of the genomic ranges considered by the callers used to generate the consensus and where appropriate, regions of interest, such as coding sequences. @@ -2440,24 +2441,24 @@
Embryonal tumors were included in non-MB and non-ATRT embryonal tumor subtyping (molecular-subtyping-embryonal
analysis module) if they met any of the following criteria:
MN1::PATZ1
since it is an entity separate of CNS HGNET-MN1 tumors139.MN1::PATZ1
since it is an entity separate of CNS HGNET-MN1 tumors140.Non-MB and non-ATRT embryonal tumors identified with the above criteria were further subtyped (molecular-subtyping-embryonal
analysis module) using the criteria below140,141,pubmed:30249036?,pubmed:26389418?.
Non-MB and non-ATRT embryonal tumors identified with the above criteria were further subtyped (molecular-subtyping-embryonal
analysis module) using the criteria below141,142,pubmed:30249036?,pubmed:26389418?.
ETMR, C19MC-altered
(Embryonal tumor with multilayer rosettes, chromosome 19 miRNA cluster altered)122,142.ETMR, NOS
(Embryonal tumor with multilayer rosettes, not otherwise specified)122,142.ETMR, C19MC-altered
(Embryonal tumor with multilayer rosettes, chromosome 19 miRNA cluster altered)123,143.ETMR, NOS
(Embryonal tumor with multilayer rosettes, not otherwise specified)123,143.CNS HGNET-MN1
[Central nervous system (CNS) high-grade neuroepithelial tumor with MN1 alteration].CNS HGNET-BCOR
(CNS high-grade neuroepithelial tumor with BCOR alteration).CNS HGNET-BCOR
(CNS high-grade neuroepithelial tumor with BCOR alteration).CNS NB-FOXR2
(CNS neuroblastoma with FOXR2 activation).CNS EFT-CIC
(CNS Ewing sarcoma family tumor with CIC alteration)121CNS EFT-CIC
(CNS Ewing sarcoma family tumor with CIC alteration)122CNS Embryonal, NOS
(CNS Embryonal tumor, not otherwise specified).Neurocytoma subtypes central neurocytoma (CNC) and extraventricular neurocytoma (EVN) were assigned (molecular-subtyping-neurocytoma
analysis module) based on the primary site of the tumor144.
+
Neurocytoma subtypes central neurocytoma (CNC) and extraventricular neurocytoma (EVN) were assigned (molecular-subtyping-neurocytoma
analysis module) based on the primary site of the tumor145.
If the tumor’s primary site was “ventricles,” we assigned the subtype as CNC; otherwise, we assigned the subtype as EVN.
Craniopharyngiomas (CRANIO) were subtyped (molecular-subtyping-CRANIO
analysis module) into adamantinomatous (CRANIO, ADAM
), papillary (CRANIO, PAP
) or undetermined (CRANIO, To be classified
) based on the following criteria145,146:
Craniopharyngiomas (CRANIO) were subtyped (molecular-subtyping-CRANIO
analysis module) into adamantinomatous (CRANIO, ADAM
), papillary (CRANIO, PAP
) or undetermined (CRANIO, To be classified
) based on the following criteria146,147:
CRANIO, PAP
.CRANIO, ADAM
.If pathology diagnosis was Subependymal Giant Cell Astrocytoma (SEGA)
, the LGG
portion of molecular subtype was recoded to SEGA
.
Lastly, for all LGG- and GNT- subtyped samples, if the tumors were glialneuronal in origin, based on pathology_free_text_diagnosis
entries of desmoplastic infantile
,desmoplastic infantile ganglioglioma
, desmoplastic infantile astrocytoma
or glioneuronal
, each was recoded as follows:
If pathology diagnosis is Low-grade glioma/astrocytoma (WHO grade I/II)
or Ganglioglioma
, the LGG
portion of the molecular subtype was recoded to GNT
.
Ependymomas (EPN) were subtyped (molecular-subtyping-EPN
analysis module) into EPN, ST RELA
, EPN, ST YAP1
, EPN, PF A
and EPN, PF B
based on evidence for these molecular subgroups as described in Pajtler et al.131.
+
Ependymomas (EPN) were subtyped (molecular-subtyping-EPN
analysis module) into EPN, ST RELA
, EPN, ST YAP1
, EPN, PF A
and EPN, PF B
based on evidence for these molecular subgroups as described in Pajtler et al.132.
Briefly, fusion, CNV and gene expression data were used to subtype EPN as follows:
RELA
as fusion partner, e.g., C11orf95::RELA
, LTBP3::RELA
, was subtyped as EPN, ST RELA
.tp53_nf1_score
analysis module)We annotated TP53 altered HGG samples as either TP53 lost
or TP53 activated
and integrated this within the molecular subtype.
To this end, we applied a TP53 inactivation classifier originally trained on TCGA pan-cancer data42 to the matched RNA expression data, with samples batched by library type.
-Along with the TP53 classifier scores, we collectively used consensus SNV and CNV, SV, and reference databases that list TP53 hotspot mutations147,148 and functional domains149 to determine TP53 alteration status for each sample.
+Along with the TP53 classifier scores, we collectively used consensus SNV and CNV, SV, and reference databases that list TP53 hotspot mutations148,149 and functional domains150 to determine TP53 alteration status for each sample.
We adopted the following rules for calling either TP53 lost
or TP53 activated
:
If a sample had either of the two well-characterized TP53 gain-of-function mutations, p.R273C or p.R248W43, we assigned TP53 activated
status.
Samples were annotated as TP53 lost
if they contained i) a TP53 hotspot mutation as defined by IARC TP53 database or the MSKCC cancer hotspots database147,148 (see also, Key Resources Table), ii) two TP53 alterations, including SNV, CNV or SV, indicative of probable bi-allelic alterations; iii) one TP53 somatic alteration, including SNV, CNV, or SV or a germline TP53 mutation indicated by the diagnosis of Li-Fraumeni syndrome (LFS)150, or iv) one germline TP53 mutation indicated by LFS and the TP53 classifier score for matched RNA-Seq was greater than 0.5.
Samples were annotated as TP53 lost
if they contained i) a TP53 hotspot mutation as defined by IARC TP53 database or the MSKCC cancer hotspots database148,149 (see also, Key Resources Table), ii) two TP53 alterations, including SNV, CNV or SV, indicative of probable bi-allelic alterations; iii) one TP53 somatic alteration, including SNV, CNV, or SV or a germline TP53 mutation indicated by the diagnosis of Li-Fraumeni syndrome (LFS)151, or iv) one germline TP53 mutation indicated by LFS and the TP53 classifier score for matched RNA-Seq was greater than 0.5.
Participant metadata included a reported gender. @@ -2552,8 +2553,8 @@
survival-analysis
analysis module)We calculated overall survival (OS) as days since initial diagnosis and performed several survival analyses on the OpenPBTA cohort using the survival
R package.
-We performed survival analysis for patients by HGG subtype using the Kaplan-Meier estimator151 and a log-rank test (Mantel-Cox test)pubmed:5910392? on the different HGG subtypes.
-Next, we used multivariate Cox (proportional hazards) regression analysis152 to model the following: a) tp53 scores + telomerase scores + extent of tumor resection + LGG group + HGG group
, in which tp53 scores
and telomerase scores
are numeric, extent of tumor resection
is categorical, and LGG group
and HGG group
are binary variables indicating whether the sample is in either broad histology grouping, b) tp53 scores + telomerase scores + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), and c) quantiseq cell type fractions + CD274 expression + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), in which quantiseq cell type fractions
and CD274 expression
are numeric.
tp53 scores + telomerase scores + extent of tumor resection + LGG group + HGG group
, in which tp53 scores
and telomerase scores
are numeric, extent of tumor resection
is categorical, and LGG group
and HGG group
are binary variables indicating whether the sample is in either broad histology grouping, b) tp53 scores + telomerase scores + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), and c) quantiseq cell type fractions + CD274 expression + extent of tumor resection
for each cancer_group
with an N>=3 deceased patients (DIPG, DMG, HGG, MB, and EPN), in which quantiseq cell type fractions
and CD274 expression
are numeric.
OpenPBTA workflows repository | this project | -153 | +154 |
OpenPBTA analysis repository | this project | -154 | +155 |
OpenPBTA manuscript repository | @@ -3105,195 +3106,198 @@