TCGA Transcriptomics End-to-End Pipeline (Download → Gene × Sample Expression Matrix)
When to Use
- You need to turn TCGA STAR-Counts expression data of any cancer project (COAD/READ/BRCA/LUAD, etc.) into a gene × sample matrix for downstream analysis such as differential expression or survival analysis.
- You already have a downloaded
gdc_datafolder (one TSV per sample) or need to download from scratch; you need to handle duplicate aliquots, normalize column names, and output multiple expression matrices (counts/TPM/FPKM).
Overview
End-to-end workflow: GDC query & download (TCGAbiolinks) → directory inspection → per-sample STAR gene counts TSV parsing → gene_id-aligned merge → file_id → TCGA barcode mapping → output counts/TPM/FPKM matrices → duplicate-aliquot deduplication (keep latest batch) → column-name normalization → validation.
Validated on TCGA-COAD (524 samples) + TCGA-READ (177 samples); all core parameters are configurable so any TCGA project can reuse it directly.
Inputs (Parameters)
| Parameter | Description | Example |
|---|---|---|
| projects | TCGA project list | c("TCGA-COAD","TCGA-READ") or "TCGA-BRCA" |
| gdc_data_dir | Root of downloaded data (skip download if exists) | <root>/gdc_data |
| query_tsv_dir | Directory of GDCquery manifests (gdc_query_<PROJ>_STAR_Counts.tsv) | <root>/Stat |
| expr_types | Output types, multi-select | c("counts","tpm","fpkm"); counts=unstranded, tpm=tpm_unstranded, fpkm=fpkm_uq_unstranded (use fpkm_unstranded if plain FPKM preferred) |
| sample_type_filter | Sample scope | NULL=all; "Primary Tumor"=primary only; c("Primary Tumor","Solid Tissue Normal") |
| keep_dup_latest | Deduplicate same-aliquot files keeping the latest batch | TRUE |
| output_format | csv / tsv | "csv" |
Steps
1. Query & Download (skip if gdc_data already exists)
- Query:
GDCquery(project, data.category="Transcriptome Profiling", data.type="Gene Expression Quantification", workflow.type="STAR - Counts"); clinical:GDCquery_clinic(project, type="clinical"). - Output: query RDS +
gdc_query_<PROJ>_STAR_Counts.tsvmanifest (file_id, cases.submitter_id, sample.submitter_id, sample_type, created_datetime, updated_datetime, analysis_id, file_size). - Download:
GDCdownload(query, directory=gdc_data_dir, files.per.chunk=100, method="api"). - Note: TCGAbiolinks >= 2.30 removed the legacy/file.type arguments; GDCdownload method is only "api"/"client"; install pitfalls in Pitfalls. Files are auto-extracted to
gdc_data/<PROJ>/Transcriptome_Profiling/Gene_Expression_Quantification/<file_uuid>/<file_uuid>.rna_seq.augmented_star_gene_counts.tsv.
2. File Format Inspection (always inspect before merging)
One TSV per sample: line 1 is a comment # gene-model: GENCODE v36 (skip); line 2 header gene_id, gene_name, gene_type, unstranded, stranded_first, stranded_second, tpm_unstranded, fpkm_unstranded, fpkm_uq_unstranded; lines 3-6 are non-gene summary rows N_unmapped/N_multimapping/N_noFeature/N_ambiguous (must skip; optionally keep for QC); from line 7 on are genes with versioned ENSEMBL IDs (ENSG00000000003.15). All files within a project share the same gene set (GENCODE v36 / hg38).
3. Merge into Gene × Sample Matrix (core)
Input: per-sample TSVs + query manifest. Logic: per file data.table::fread(comment.char="#") → filter grepl("^ENSG", gene_id) → align to the first file's gene order with match() → extract needed columns → bind into a matrix.
- Sample column names = TCGA barcode: the UUID folder is the GDC file_id; map file_id → cases.submitter_id (case-level barcode, e.g. TCGA-A6-2670) or sample.submitter_id (aliquot-level) via the manifest.
- Key pitfall: matrix column names contain hyphens; never build the matrix via as.data.frame() (make.names converts - to ., breaking clinical-table matching); use do.call(cbind, list) + explicit colnames()<-, or as.data.frame(check.names=FALSE).
- Output: three matrices (gene_id + gene_name annotation columns + sample columns) + sample_mapping.csv (file_uuid, barcode, aliquot, project, sample_type, created_datetime, updated_datetime, total_counts).
4. Duplicate-Aliquot Identification & Verification
- Group sample_mapping by aliquot (sample.submitter_id): N>1 means the same sample was sequenced/analyzed more than once. In COAD/READ there were 10 groups (2 files each, independent GDC analyses of the same aliquot, different created_datetime batches).
- Do NOT confuse with same-case different-aliquot files (e.g. tumor + adjacent normal; 52 cases in COAD/READ) — those are different biological samples, not duplicate sequencing.
- Optional consistency check: per pair compute Pearson (raw counts, log2(counts+1), TPM) and "exactly identical" fraction. Observed: log2 correlation 0.92-0.96 (largely concordant) but raw/TPM correlation low (0.2-0.6 / 0.02-0.15) because a few technical-artifact genes dominate — RN7SL1 (structural RNA; multimapping assignment differs between STAR versions, up to millions of counts), mitochondrial genes (MT-ND4/CO1/CO2/CO3/CYB), and a few highly expressed genes (EEF1A1, RPS6, EEF2). Judge consistency by log2-counts correlation.
5. Deduplication (keep the latest batch per aliquot)
Priority: later created_datetime (analysis batch, UTC-normalized) → later updated_datetime → higher total counts (deeper library). All 10 pairs in the real run were decided by the time criterion. Rename the kept column to the standard aliquot barcode (e.g. TCGA-A6-2677-01A); record dropped columns in dedup_mapping.csv (keep/drop columns, timestamps, total counts, reason).
6. Column-Name Normalization & Validation
- Merge-phase collision names of the form
TCGA-XX-XXXX_TCGA-XX-XXXX-XXA(non-duplicate samples, label-only issue) are normalized to the pure aliquot barcode; log in colname_normalization.csv. - Final validation: dimensions (n_samples = total files − dropped columns), unique column names, all matching the standard barcode regex
^TCGA-[A-Z0-9]{2}-[A-Z0-9]{4}(-[0-9]{2}[A-Z])?$, no NA.
Parameters to Customize
output_dir <- "<OUTPUT_ROOT>" # required
gdc_data_dir <- "<GDC_DATA_DIR>" # required
query_tsv_dir <- "<QUERY_TSV_DIR>" # required
projects <- c("TCGA-COAD","TCGA-READ")
expr_types <- c("counts","tpm","fpkm")
sample_type_filter <- NULL
keep_dup_latest <- TRUE
output_format <- "csv"
Reference scripts: scripts/tcga_expression_pipeline.R.txt (main workflow; edit the config block only); scripts/verify_duplicates_template.py (optional verification).
Evidence Criteria / Expected Output
- Query non-empty; downloaded file count equals manifest file count;
du -shtotal matches manifest file_size sum. - Matrix dimensions = gene count (~60k) × sample count; no NA; unique standard TCGA barcode column names.
- Spot check: a gene value (e.g. TSPAN6/ENSG00000000003.15) in one sample matches the source TSV.
- Post-dedup sample count = total files − extra files from duplicate aliquots; dedup_mapping fully documents each decision.
Pitfalls
- biomaRt install failure: under Bioc 3.22, dplyr < 1.2 lacks the filter_out export causing lazy-loading failure; run
install.packages("dplyr")before installing TCGAbiolinks. - TCGAbiolinks >= 2.30 API: GDCquery has no legacy argument; clinical via GDCquery_clinic; GDCdownload method only api/client; result column is analysis_workflow_type (underscore).
- fread needs
comment.char="#"to skip the# gene-modelline; filter N_* summary rows with^ENSG. - Hyphenated column names get mangled by make.names → use do.call(cbind)/check.names=FALSE.
- Duplicate-aliquot grouping is by aliquot, not by case; same-case different-aliquot files are different samples.
- RN7SL1/mitochondrial genes are technical artifacts; do not judge replicate quality from them; prefer time over depth for dedup (unless timestamps are unavailable).
- GDC download progress uses \r carriage returns; monitor with
find <dir> -name "*.tsv" | wc -l+du -sh. - Memory: read files one at a time keeping only needed columns; three matrices of ~60k × 701 samples take ~1 GB, safely within a 50 GB limit.
Extensions
- Sample filtering: by sample_type (primary tumor only), tumor purity (ESTIMATE/ABSOLUTE), exclude FFPE.
- Clinical integration: merge clinical_<PROJ>.tsv (GDCquery_clinic) on case barcode prefix for grouping and survival analysis.
- Downstream: feed counts to DESeq2/edgeR (strip version suffixes from IDs or map symbols via gene_name); TPM/FPKM for cross-sample comparisons and visualization.
- Pan-cancer merge: run per project then merge on common genes (mind batch effects; correct them).
- File-level QC: extract N_* summary rows for alignment/multimapping-rate QC and sample filtering.
- Output format extensions: Excel via writexl, TSV via fwrite(sep="\t"), RDS archives for R workflows.
> 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.
