TCGA 通用生存分析（突变样本 + 表达分组 + OS）

When to use

• User asks: TCGA 某癌种、某基因突变样本里、某表达基因高低组的生存分析（原需求：TCGA-LUAD + EGFR 突变 + COL1A1/STAB1 表达与预后）。
• Or: "把 TCGA_analysis.R 泛化" / "换个癌种/基因跑生存分析"。

Canonical script

• /workspace/yanyali/other_work/PM-202396-hongxuan-288/TCGA-LUAD_analysis/TCGA_analysis_general.R.txt
  (copy in this skill: scripts/TCGA_analysis_general.R.txt; original TCGA_analysis.R is the LUAD/EGFR/COL1A1-STAB1 hardcoded version — NEVER modify it, only use/generalize the general version)
• R environment: Rscript is NOT in PATH — use /workspace/yanyali/miniconda3/envs/seuratv4/bin/Rscript (has TCGAbiolinks, survminer, ggplot2, survival, dplyr; no optparse → base-R arg parsing).

Usage

R=/workspace/yanyali/miniconda3/envs/seuratv4/bin/Rscript
$R TCGA_analysis_general.R.txt --project TCGA-LUAD --mut_gene EGFR \
    --exp_genes COL1A1,STAB1 --quantiles 0.25,0.5,0.75 \
    --outdir TCGA-LUAD_EGFR_analysis
# 换癌种/基因:
$R TCGA_analysis_general.R.txt --project TCGA-BRCA --mut_gene TP53 \
    --exp_genes ESR1,ERBB2 --outdir TCGA-BRCA_TP53_analysis

Params: --project (any TCGA-*), --mut_gene (Hugo_Symbol, or ALL = any-mutated samples), --exp_genes (comma list), --quantiles (comma list, default 0.5), --outdir, --clinical_url (override Xena URL), --sample_types (default 1 = primary tumor), --exclude (comma barcodes), --overwrite (redownload).

Data flow (same as original)

1. Mutation: GDCquery "Simple Nucleotide Variation"/"Masked Somatic Mutation"/"Aliquot Ensemble Somatic Variant Merging and Masking" → filter Hugo_Symbol == mut_gene; sample = substr(Tumor_Sample_Barcode,1,15).
2. Expression: GDCquery "Transcriptome Profiling"/"Gene Expression Quantification"/"STAR - Counts" → assay tpm_unstrand → log2(TPM+1); rowData gene_name as rownames.
3. Clinical: UCSC Xena clinicalMatrix, auto URL https://tcga-xena-hub.s3.us-east-1.amazonaws.com/download/TCGA.<CANCER>.sampleMap/<CANCER>_clinicalMatrix (verified for LUAD, BRCA; download.file method="libcurl"). Keep sample_type_id in --sample_types; event from vital_status; days: death→days_to_death, censor→days_to_last_followup; OS_time = days/30.4375 (months).
4. Merge meta + expression by 15-char barcode (inner_join sampleID=Sample) → per gene × quantile: split High/Low at quantile, survfit(Surv(OS_time,event)~Expression), ggsurvplot PDF, logrank p via survdiff → results CSV.

Outputs

<project>_<mut_gene>_mutations.txt, <project>_TPM_exp.rds (cache), <project>_mutations.rds (cache), <cancer>_clinicalMatrix.txt (cache), <project>_<mut_gene>.samples_OS_meta.rds (merged), <gene>_quantile<q>_expression_OS_plot.pdf (KM + risk table), <project>_<mut_gene>_survival_results.csv (gene, quantile, n_samples, n_high, n_low, logrank_p).

Pitfalls (learned the hard way)

• survminer named palette → "No shared levels found between names(values)..." warnings even when names match strata ("Expression=High" etc.). Fix: use an UNNAMED 2-color vector in level order c("#2166AC", "#B2182B") (Low=blue, High=red) — zero warnings. (Verified: named → 22 warnings, unnamed → 0.)
• GDC downloads are huge (whole-project MAF + expression). Script caches rds/txt in outdir; reuse unless --overwrite. Real GDC run takes a long time — don't run it as part of a quick test.
• Barcode lengths differ: mutation MAF = 16-char TCGA-XX-XXXX-01A, expression/Xena = 15-char TCGA-XX-XXXX-01 → always substr(...,1,15) before matching.
• vital_status encoding varies by cohort (LIVING/DECEASED vs ALIVE/DEAD vs 0/1) — build_meta normalizes all.
• LUAD had duplicate-sample problems (original script hardcoded 4 exclusions: TCGA-44-6147-01A/01B, TCGA-44-5645-01A/01B) — the general script exposes --exclude AND auto-dedups: the STAR expression matrix can contain multiple aliquots (01A/01B) of the same 15-char sample → duplicated Sample values inflate inner_join results. exp_df[!duplicated(exp_df$Sample), ] removes them with a printed count. Verified: 79 EGFR-mut samples with expression → 4 duplicates removed → 75 → 73 with survival info.
• Reusing the original run's caches: the old TCGA-LUAD_TPM_exp.rds format is list(exp = raw TPM matrix (Ensembl rownames), gene = gene_names) — the new script's cache-hit path consumes it directly (positional gene match). Values are raw TPM (not log2) but quantile grouping is monotonic → identical High/Low groups and p-values.
• mut_gene not found → stop() with a hint to check Hugo_Symbol spelling; exp gene not found → warning + skip (error if all missing).
• Testing without GDC downloads: source("TCGA_analysis_general.R.txt") only loads functions (main is gated by sys.nframe()==0), then test with the real local LUAD_clinicalMatrix.txt + mock expression data — verifies build_meta, survival_plot, PDF output, and catches palette/strata issues. Real LUAD clinical: 511 primary tumor samples (184 deceased / 327 living), OS_time 0–238 months.
• Fontconfig warning from cairo_pdf on this cluster is harmless.

Real-data context (PM-202396)

Original run outputs live in /workspace/yanyali/other_work/PM-202396-hongxuan-288/TCGA-LUAD_analysis/ (GDCdata/ [2.5GB: 618 mutation MAFs + 601 STAR files], TCGA_LUAD_EGFR.txt, TCGA-LUAD_TPM_exp.rds, TCGA-LUAD_EGFR.samples_OS_meta.rds, LUAD_clinicalMatrix.txt, 6 KM PDFs). Reuse these cached files (copy into outdir or point --outdir there) to skip re-downloading.
Test run of the general script (validated end-to-end): TCGA-LUAD_analysis/TCGA-LUAD_EGFR_test/ — caches: TCGA-LUAD_mutations.rds (built once via GDCquery+GDCprepare reading local GDCdata, 194729x140), TCGA-LUAD_TPM_exp.rds, LUAD_clinicalMatrix.txt. Results: 79 EGFR-mut samples with expression → 4 dups removed → 75 → 73 with survival (30 events); COL1A1 q0.25 p=0.098 / q0.5 p=0.405 / q0.75 p=0.877; STAB1 q0.25 p=0.023 / q0.5 p=0.078 / q0.75 p=0.631; 6 KM PDFs + TCGA-LUAD_EGFR_survival_results.csv. Sample set identical to the original run's 73.

> Packaging note. Scripts and tables in this skill carry a .txt suffix on their
> original extension (e.g. .R.txt, .tsv.txt) because the SkillHub upload allowlist
> does not accept .R or .tsv. This is cosmetic: Rscript is extension-agnostic, so
> Rscript scripts/<name>.R.txt ... behaves exactly like .R, and tab-separated data
> parses identically from .tsv.txt. Keep the suffix when repackaging.