Classify Head and Neck Squamous Cell Carcinoma Samples into TCGA Subtypes

This function takes a gene expression matrix, performs quality checks, optional log2 transformation, single-sample gene set enrichment (ssGSEA) at the pathway level, and then predicts the TCGA-based molecular subtype of each sample using a pre-trained random forest model.

Usage

classifyHNSC(
  input_expr = NULL,
  idType = "SYMBOL",
  outputType = c("class", "prob")
)

Arguments

input_expr

A numeric gene expression matrix (or data frame) with genes in rows and samples in columns. Row and column names are mandatory. The expression values must be non-negative (e.g., TPM, FPKM, or normalized counts). If the data are not on a log2 scale, they will be automatically log2-transformed via log2(x + 1) internally.

idType

Character string specifying the gene identifier type used in the rownames of input_expr. Acceptable values are "SYMBOL" (default), "ENSEMBL", "ENTREZID", or "REFSEQ". If any type other than "SYMBOL" is provided, the function will convert the rownames to gene symbols using the org.Hs.eg.db annotation package. This requires org.Hs.eg.db and AnnotationDbi to be installed (they are listed under Suggests of the package). For multiple identifiers mapping to the same gene symbol, expression values are averaged.

outputType

Character indicating the type of prediction result to return. Use "class" (default) to obtain the assigned subtype label for each sample, or "prob" to get a matrix of class probabilities (one row per sample, one column per subtype).

Value

• If outputType = "class", a named character vector of predicted subtypes (e.g., "Atypical", "Basal", "Classical", "Mesenchymal").

• If outputType = "prob", a numeric matrix with samples in rows, subtypes in columns, and values ranging from 0 to 1 representing the predicted probability of each subtype.

Details

The function proceeds through five main steps:

1. Input validation: Checks for missing values (NA), non-numeric columns, negative expression values, missing row/column names, and correct sample number.

2. Gene ID conversion (if idType is not "SYMBOL"): Converts row identifiers to gene symbols using org.Hs.eg.db. If multiple identifiers map to the same gene symbol, the maximum expression value is retained for each sample.

3. Log2 transformation: If the data are not already on a log2 scale (determined by an internal heuristic), they are transformed as log2(x + 1).

4. Pathway-level representation: Single-sample gene set enrichment analysis (ssGSEA) is performed using a curated set of pathway gene sets stored internally as required.sets. This step harmonizes the input data with the training cohort at the biological pathway level.

5. Random forest classification: The pathway scores are centered, scaled, and fed into a pre-trained random forest model (finalModel) to generate the final subtype prediction.

Both required.sets and finalModel are internal package data and are loaded automatically. The function requires the GSVA package.

Note

If using idType values other than "SYMBOL", the packages org.Hs.eg.db and AnnotationDbi must be installed. You can install them with:

BiocManager::install("org.Hs.eg.db")
BiocManager::install("AnnotationDbi")

The conversion step may remove genes that lack a unique symbol mapping; a warning will indicate the number of genes retained.

See also

gsva for the underlying pathway scoring engine.

Examples

if (FALSE) { # \dontrun{
# Example with gene symbols (default)
data(TCGA_LUSC)
subtypes <- classifyHNSC(TCGA_LUSC, outputType = "class")

# Example with Ensembl IDs
# Ensure org.Hs.eg.db is installed
subtypes_ens <- classifyHNSC(TCGA_LUSC_ENSEMBL,
                             idType = "ENSEMBL",
                             outputType = "prob")
} # }