Thyroid Cancer Multi-Omics: Clinical Covariates, Transcriptomic Programs, Somatic Drivers, and Precision-Oncology Integration
Expanded technical portfolio case study of non-medullary thyroid cancer integrating clinical phenotype engineering, covariate-adjusted RNA-seq, pathway/signaling analysis, inflammatory programs, somatic driver architecture, demographic associations, and staged DNA–RNA integration.
This project is not a single differential-expression analysis. It was designed around four linked precision-oncology questions:
- Sociodemographic genomic disparities — do somatic drivers, molecular groups, pathway alterations, and mutation burden differ by race, age, sex, or disease status?
- Immune and inflammasome architecture — do inflammatory/immune programs differ by AWD/NED, autoimmune thyroiditis, or demographics?
- Prognostic pathway stratification — are MAPK, PI3K/AKT/mTOR, stress-response, and other tumor programs associated with AWD versus NED?
- Multi-omics integration / validation — do DNA driver states map to RNA phenotypes, and can major findings be compared with TCGA-THCA after appropriate harmonization?
| Layer | Aggregate real project value |
|---|---|
| Final non-medullary RNA AWD/NED cohort | 130 |
| NED | 91 |
| AWD | 39 |
| RNA samples with reported race | 107 |
| Genes tested | 24,364 |
| DNA MAF samples | 139 |
| DNA samples with clinical metadata | 121 |
| DNA AWD/NED subset | 101 |
| Low-depth DNA samples flagged | 20 |
Clinical metadata
AWD/NED + race + age + sex + histology + thyroiditis evidence
|
cohort / endpoint harmonization
|
+-------------+-------------+
| |
RNA-seq DNA
Salmon gene counts Funcotator MAFs
| |
count QC / filtering variant QC / reconciliation
| |
DESeq2 adjusted model curated driver layer
age + sex + disease status |
| patient alteration matrix
VST / PCA |
| MAPK / PI3K / repair
gene-level effects |
| clinical exact models
Hallmark / GO / KEGG / Reactome |
| |
thyroid + stress + immune programs |
+-------------+-------------+
|
matched DNA–RNA
integration
The RNA analysis does not treat the clinical table as passive metadata.
NED is the reference and AWD is the disease group of interest. Endpoint mapping is kept auditable because biochemical/structural incomplete response, recurrence, persistence, or similar labels can materially change the phenotype.
Histology is refined using structured fields and pathology text. Medullary thyroid cancer is excluded. A study-level covariate collapses histology to:
- Papillary / Follicular / FVPTC
- All other non-medullary
The analysis retains:
- age as a continuous covariate;
- age groups
<45,45–54,55–64,65+; - a
<55versus55+indicator; - sex;
- harmonized race group.
A separate clinical model evaluates AWD in relation to available demographics and histology:
glm(
AWD ~ age_at_biopsy + sex + race_group + study_histology_group,
family = binomial()
)See docs/COVARIATE_MODELING.md.
Salmon-derived counts are validated, aligned to clinical metadata, and filtered before DESeq2.
The primary source model is:
design = ~ age_at_biopsy + sex + disease_groupSamples missing required model covariates are excluded explicitly rather than silently retained.
DESeq2 outputs include:
- all genes;
- BH-FDR significant genes;
- strong-effect genes (
FDR < 0.05,|log2FC| >= 1); - genes higher in AWD;
- genes higher in NED;
- normalized counts;
- VST expression;
- PCA coordinates.
| Result | N |
|---|---|
| Genes tested | 24,364 |
| FDR < 0.05 | 354 |
| Higher in AWD, FDR < 0.05 and |log2FC| ≥ 1 | 110 |
| Higher in NED, FDR < 0.05 and |log2FC| ≥ 1 | 41 |
| Race-adjusted FDR < 0.05 sensitivity result | 230 |
The race-adjusted result is kept conceptually separate from the primary age/sex-adjusted source model.
The project evaluated effect size together with descriptive discrimination rather than ranking genes on P value alone.
| Gene | log2FC | AUC | 95% CI |
|---|---|---|---|
| IL1A | 4.59 | 0.67 | 0.570–0.777 |
| SLC2A1 | 2.30 | 0.64 | 0.528–0.749 |
| GJB6 | 3.18 | 0.61 | 0.505–0.721 |
| GJB2 | 3.37 | 0.60 | 0.495–0.714 |
| ZPLD1 | 2.55 | 0.59 | 0.483–0.697 |
ROC analyses are discovery-level and are not presented as clinically validated classifiers.
DESeq2 statistics are used for preranked FGSEA across:
- Hallmark
- GO Biological Process
- KEGG
- Reactome
The analysis tracks normalized enrichment score, FDR, and whether a pathway is enriched toward AWD or NED.
The project also evaluates coordinated sample-level programs:
| Program | Genes detected | NED median | AWD median | FDR |
|---|---|---|---|---|
| Mitochondrial metabolism | 10 | -0.0009 | 0.148 | 0.590 |
| Oxidative stress and ferroptosis | 10 | -0.0343 | 0.112 | 0.590 |
| MAPK signaling | 12 | -0.212 | -0.0415 | 0.749 |
| PI3K/AKT signaling | 10 | -0.350 | -0.322 | 0.749 |
| Thyroid differentiation / iodine metabolism | 11 | 0.052 | -0.0481 | 0.846 |
| EMT / invasion | 10 | 0.0387 | 0.0033 | 0.896 |
| Inflammation / interferon response | 10 | 0.0052 | -0.130 | 0.915 |
No prespecified program remained FDR-significant in the unadjusted AWD/NED program-level comparison. Mitochondrial and oxidative-stress distributions nevertheless motivated covariate-adjusted follow-up.
See docs/PATHWAY_AND_SIGNALING.md.
The prespecified inflammasome panel contains:
NLRP3, NLRP1, NLRC4, AIM2, PYCARD, CASP1, CASP4, CASP5, IL18, GSDMD, NFKB1, RELA, TXNIP.
| Gene | NED median | AWD median | Difference | P | FDR |
|---|---|---|---|---|---|
| CASP4 | 7.95 | 8.58 | +0.627 | 0.0154 | 0.200 |
| GSDMD | 7.72 | 7.86 | +0.132 | 0.121 | 0.514 |
| NLRP3 | 7.29 | 7.02 | -0.270 | 0.206 | 0.535 |
| CASP1 | 7.04 | 7.17 | +0.126 | 0.382 | 0.828 |
CASP4 is the strongest individual-gene signal, but it does not remain FDR-significant.
The composite inflammasome score is constructed from gene-wise standardized VST expression.
Aggregate project results include:
- race association: Kruskal
P=0.0155, FDR0.031; - Black vs White pairwise BH
P=0.019; - adjusted Black vs White coefficient
β=-0.439,P=0.0285; - adjusted AWD coefficient
β=0.367,P=0.0317.
These results illustrate why covariate-aware modeling changes the interpretation of pathway-level signals.
Autoimmune-thyroiditis evidence is defined conservatively using pathology and/or clearly positive thyroid-autoantibody evidence. When immune-deconvolution outputs are available, cell-type estimates can be evaluated by thyroiditis, disease status, race, age group, and sex.
See docs/INFLAMMASOME_AND_IMMUNE.md.
The DNA workflow reconciles provisional and curated calls into an analysis-ready driver layer and collapses them to patient-level alteration matrices.
| Gene | Mutated | Prevalence |
|---|---|---|
| BRAF | 68 | 48.9% |
| HRAS | 10 | 7.2% |
| NRAS | 9 | 6.5% |
| PIK3CA | 5 | 3.6% |
| PTEN | 5 | 3.6% |
| KRAS | 4 | 2.9% |
| RET | 4 | 2.9% |
| TP53 | 4 | 2.9% |
Rare alterations are analyzed preferentially at the pathway level rather than overinterpreted individually.
The project defines three patient-level states:
| Group | N | % |
|---|---|---|
| BRAF V600E | 64 | 46.0% |
| Other / no canonical MAPK hotspot | 56 | 40.3% |
| RAS hotspot | 19 | 13.7% |
This grouping is used for clinical interpretation and future DNA–RNA pathway integration.
The public configuration explicitly records pathway membership.
BRAF, NRAS, HRAS, KRAS, NF1, RET
PIK3CA, AKT1, PTEN, MTOR, TSC1, TSC2
TP53, RB1, APC
ATM, ATR, CHEK2, BRCA1, BRCA2
DICER1, MEN1
This makes the biological decision logic reviewable instead of burying it in code.
In the current mutation-only analysis, canonical MAPK molecular-group distribution differed between Black and White patients:
Fisher exact P = 0.000583.
The observed pattern suggested relatively greater RAS-hotspot and lower BRAF V600E representation among Black/African American patients in this cohort.
This is a discovery-level association, not a causal ancestry claim. The project explicitly calls for adjustment/sensitivity analysis incorporating histology, age, sex, sequencing depth, and sample size.
Current mutation-only observations indicate that:
- BRAF prevalence is similar in AWD and NED;
- NRAS is concentrated more heavily in AWD in the plotted subset;
- RET appears more frequent in NED in the plotted subset;
- sparse drivers require exact/pathway-level modeling.
The prognosis question cannot be fully resolved from mutation-only data because CNV and fusion events are not yet fully integrated.
The matched multi-omics layer is designed to compare DNA states such as BRAF/RAS/PI3K with RNA phenotypes including:
- MAPK signaling
- mTORC1 / PI3K activity
- E2F / G2M / MYC
- glycolysis
- oxidative phosphorylation
- EMT
- interferon / inflammation
- oxidative stress
- thyroid differentiation
See docs/MULTIOMICS_INTEGRATION.md.
A major principle of this public repository is not to present planned work as completed.
- cohort / metadata harmonization
- histology refinement and medullary exclusion
- clinical Table 1 / missingness summaries
- adjusted AWD demographic model
- DESeq2 + VST/PCA
- race-adjusted RNA sensitivity result
- FGSEA
- thyroid/stress program scoring
- inflammasome gene / composite analyses
- MAF/driver reconciliation
- primary driver prevalence
- canonical MAPK molecular groups
- discovery-level race-stratified MAPK comparison
- thyroiditis annotation
- immune deconvolution
- expanded robust demographic sensitivity modeling
- full CNV integration
- full fusion/rearrangement integration
- final mutational-signature layer
- fully integrated mutation+CNV+fusion outcome models
- TCGA-THCA validation
See docs/COMPLETION_STATUS.md.
config/
analysis_parameters.yml
gene_sets.json
data/synthetic/
docs/
scripts/R/
results/
aggregate_real_results/
example_tables/
scripts/python/
.github/workflows/
The synthetic public layer is designed to exercise the workflow without redistributing restricted data.
python scripts/python/validate_inputs.py
python scripts/python/privacy_check.py
Rscript scripts/R/01_build_clinical_cohort.R
Rscript scripts/R/02_table1_and_missingness.R
Rscript scripts/R/03_adjusted_awd_clinical_model.R
Rscript scripts/R/04_prepare_counts.R
Rscript scripts/R/05_deseq2_primary_awd_ned.R
Rscript scripts/R/06_vst_pca.R
Rscript scripts/R/09_thyroid_program_scores.R
Rscript scripts/R/13_driver_reconciliation.R
Rscript scripts/R/14_mapk_molecular_groups.R
Rscript scripts/R/15_driver_prevalence.RClinical bioinformatics: endpoint harmonization, histology curation, missing-data review, covariate construction, logistic regression.
RNA-seq: Salmon counts, DESeq2, LFC shrinkage, VST/PCA, adjusted differential expression, effect-size prioritization.
Pathway biology: Hallmark, GO, KEGG, Reactome, FGSEA, sample-level pathway scoring, MAPK, PI3K/AKT, DNA repair, oxidative stress, hypoxia, EMT, interferon, thyroid differentiation.
Tumor immunology: thyroiditis phenotype engineering, immune deconvolution scaffolding, inflammasome gene and composite-score analysis.
Cancer genomics: MAF processing, curated driver reconciliation, hotspot grouping, mutation prevalence, pathway-level alteration matrices, exact testing.
Multi-omics: DNA–RNA matched analysis design, molecular-group to transcriptional-program integration, external-cohort harmonization strategy.
Statistical rigor: covariate adjustment, sparse-data exact tests, BH-FDR, sensitivity analyses, discovery-versus-validation distinction.
No real patient-level identifiers, sequencing matrices, institutional paths, protected clinical records, or restricted CARIS data are distributed. Real results are limited to aggregate non-patient-level summaries; synthetic data are used for public reproducibility.