Phenome-Wide Association Study of SPTLC3

Introduction

About PheWAS

In this Phenome-Wide Association Study, logistic regression is used to evaluate associations between single-nucleotide polymorphisms (SNPs) in the SPTLC3 gene and phenotypes derived from clinical ICD10 diagnosis billing codes. For each logistic regression model, we estimate whether the genotype is associated with a given diagnosis.

About Phecodes

Phecodes are standardized phenotype groupings that aggregate thousands of individual diagnosis codes from electronic health records into broader, clinically interpretable categories. For example, multiple ICD codes corresponding to related diabetes diagnoses can be collapsed into a single phecode category. This approach reduces sparsity and coding heterogeneity while improving the detection of clinically meaningful genotype–phenotype associations. Phecodes can be mapped back to ICD codes using the Phecode catalog.

About All of Us

All of Us (AoU) is a large-scale, longitudinal research program launched by the National Institutes of Health (NIH) in 2015. The initiative enrolled its first participants in 2017 and has since expanded to include more than 880,000 individuals across the United States. AoU is particularly valuable for genetic and phenomic research because it integrates multiple data modalities, including participant surveys, physical measurements, wearable device data, electronic health records (EHRs), and genomic sequencing data.

• About the Researcher Hub
• Browse aggregate-level data

About SPTLC3

Sphingolipids are bioactive lipid molecules that serve as structural components of cell membranes and participate in cellular signaling pathways. Their de novo biosynthesis is initiated by serine palmitoyltransferase (SPT), an enzyme complex that catalyzes the first committed step in sphingolipid production. The composition of the SPT complex influences the sphingoid base products generated:

• SPTLC1 + SPTLC2 produces typical 18-carbon sphingolipids.
• SPTLC1 + SPTLC3 produces atypical 16-carbon sphingolipids.

These atypical sphingolipids are present in many tissues but their biological function is not fully understood. In 2024, two independent studies found that reducing SPTLC3 activity had a protective effect on liver and cardiac cells. In cardiac tissue, Kovilakath et al. reported that SPTLC3 contributes to ischemic cardiomyopathy by regulating mitochondrial complex I activity, with cardiomyocyte-specific SPTLC3 depletion showing protective effects after myocardial infarction (Kovilakath, 2024). In hepatic tissue, Montefusco et al. found that SPTLC3 regulates plasma membrane sphingolipid composition and supports hepatic gluconeogenesis through adenylate cyclase/cAMP signaling (Montefusco, 2024). Together, these findings suggest that genetic variation in SPTLC3 may have physiologically relevant consequences for metabolic and cardiovascular related phenotypes.

---

Methods

SNP Selection

Candidate variants were identified using the All of Us Data Browser. Variants in exonic regions that had an allele frequency greater than 0.015 were selected for analysis. Two additional intronic variants of interest, rs117004417 and rs77523068, were also included, resulting in a final set of eight SNPs.

• Allele Count: The total number of times a specific allele (variant) is observed across all chromosomes in a dataset.
• Allele Number: The total number of valid alleles sampled for a given variant, the combined pool of all tested chromosomes.
• Allele Frequency: The proportion of all alleles in the population that are the specific variant. (AC/AN)
• Homozygous Count: The total number of individuals who possess two identical copies of a specific variant.

variant_metadata |> 
  select(rsid, 
         exon,
         consequence,
         allele_count,
         allele_num,
         homozygous_count,
         aou_link
  ) |> 
  mutate(aou_link = paste0("<a target='_blank', href='", aou_link, "'>AoU</a>"),
         allele_count = scales::comma(allele_count),
         allele_num = scales::comma(allele_num),
         homozygous_count = scales::comma(homozygous_count)
  ) |> 
  datatable(escape = FALSE, 
            options = list(dom = "t")
  )

SNP Location

Scroll right on the image to see the full gene. Click the image to view an interactive version on NCBI.

Cohort Construction

A separate cohort was generated for each variant. Participants carrying at least one copy of the alternate allele were classified as carriers, whereas participants with no copies were classified as non-carriers.

Covariates were incorporated into the statistical model to account for participant-level factors that may confound genotype–phenotype associations. The adjusted covariates included current age, sex assigned at birth, electronic health record length, and the first 10 genetic principal components. Participants with missing covariate values were excluded prior to analysis.

After covariate filtering, each variant-specific cohort was stratified by genetic ancestry group. This produced up to six ancestry-specific subgroups per variant:

• African (AFR)
• Americas (AMR)
• East Asian (EAS)
• European (EUR)
• Middle Eastern (MID)
• South Asian (SAS)

PheWAS

Phenome-wide association analyses were performed using PheTK, an open-source toolkit developed for PheWAS workflows within the All of Us Researcher Workbench. Phecodes were generated from ICD diagnosis codes using Phecode X with U.S. ICD mappings. Diagnosis codes were mapped to phecodes and aggregated at the participant–phecode level.

For each SPTLC3 variant, participants were classified using a carrier-based dominant genotype model. Individuals with genotype 0/1 or 1/1 were classified as variant carriers, and individuals with genotype 0/0 were classified as non-carriers. Covariates were added to each variant-specific cohort, including current age, sex assigned at birth, EHR length, genetically inferred ancestry, and the first 10 genetic principal components. Genetic ancestry was used to stratify analyses, whereas current age, sex assigned at birth, EHR length, and the first 10 genetic principal components were included as regression covariates. Participants with missing covariate data were excluded.

Analyses were stratified by genetically inferred ancestry group. Within each variant–ancestry cohort, logistic regression was used to test the association between variant carrier status and each phecode-defined phenotype. Carrier status was specified as the independent variable of interest, and models were adjusted for available covariates. Phecodes were tested only if at least 50 cases were present in the cohort. Participants were classified as phecode cases only if they had at least two recorded instances of the corresponding diagnosis-derived phecode.

---

Results

This section summarizes ancestry-stratified PheWAS results for each variant. Model convergence is shown first as a quality-control metric. Results should be interpreted cautiously when model convergence is low or case counts are sparse. Manhattan plots show the strength of association between the variant and each tested phecode. Larger points are statistically significant.

rs243887

Convergence

eur

amr

afr

eas

sas

mid

rs243888

Convergence

eur

amr

afr

eas

sas

mid

rs117004417

Convergence

eur

amr

afr

eas

sas

mid

rs77523068

Convergence

eur

amr

afr

eas

No converged models for this ancestry.

sas

mid

No converged models for this ancestry.

rs6109692

Convergence

eur

amr

afr

eas

sas

mid

rs77696068

Convergence

eur

amr

afr

eas

sas

mid

rs6078938

Convergence

eur

amr

afr

eas

sas

mid

rs61738161

Convergence

eur

amr

afr

eas

sas

mid

---

Meta Analysis: Manhattan Plots

This section presents meta-analysis Manhattan plots combining PheWAS results across all ancestry groups for each variant. Pooled estimates were derived using a random-effects model (REML). Results with fewer contributing ancestry groups should be interpreted with additional caution.

rs243887

rs243888

rs117004417

rs77523068

rs6109692

rs77696068

rs6078938

rs61738161

---

Meta Analysis: Category Heatmap

This section displays meta-analysis results organized by phecode category, with one tab per category. Each cell represents a variant–phecode pair, colored by the pooled beta estimate (log odds ratio).

• Blue indicates a negative association (decreased risk).
• Red indicates a positive association (increased risk).

Blood/Immune

Cardiovascular

Congenital

Dermatological

Endocrine/Metab

Gastrointestinal

Genetic

Genitourinary

Infections

Mental

Muscloskeletal

Neoplasms

Neurological

Pregnancy

Respiratory

Sense organs

Symptoms

---

Meta Analysis: QC

QQ Plot

For each variant, we rank all phecode p-values from most to least significant and plot the observed -log10(p) against what we expect if every phecode were a true null (i.e., no association anywhere). The bulk of phecodes should fall along the diagonal because most phecodes shouldn’t associate with any variant. Departure above the line at the upper right suggests a signal (or inflation). Departure below the line means deflation.

qq_data <- meta_results |>
  filter(!is.na(pval)) |>
  group_by(rsid) |>
  arrange(pval, .by_group = TRUE) |>
  mutate(
    expected = -log10(seq_along(pval) / (n() + 1)),
    observed = -log10(pval),
    lambda = median(qchisq(1 - pval, df = 1), na.rm = TRUE) / qchisq(0.5, df = 1)
  ) |>
  ungroup()

# diagonal reference line spanning the full range
ref_max <- max(c(qq_data$expected, qq_data$observed), na.rm = TRUE)

plot_ly(qq_data,
        x = ~expected,
        y = ~observed,
        color = ~rsid,
        type  = "scatter",
        mode  = "markers",
        text  = ~paste0("<b>rsid</b>: ", rsid, "<br>",
                        "<b>Phecode</b>: ", phecode, "<br>",
                        "<b>Category</b>: ", phecode_category, "<br>",
                        "<b>Description</b>: ", phecode_string, "<br>",
                        "<b>Cases</b>: ", scales::comma(meta_total_cases), "<br>",
                        "<b>Controls</b>: ", scales::comma(meta_total_controls)),
        hoverinfo = "text"
  ) |>
  add_segments(
    x = 0, 
    xend = ref_max,
    y = 0, 
    yend = ref_max,
    line      = list(color = "firebrick", width = 1.5, dash = "dash"),
    inherit   = FALSE,
    showlegend = FALSE,
    hoverinfo  = "skip"
  ) |>
  layout(
    title  = "QQ Plot by Variant",
    xaxis  = list(title = "Expected -log10(p)"),
    yaxis  = list(title = "Observed -log10(p)"),
    legend = list(orientation = "h", x = 0.5, y = -0.2, xanchor = "center")
  )

Ancestries Per Meta-Analysis

In this plot, we show how often we actually had multi-ancestry data. For each variant, we counted how many phecode-variant pairs were meta-analyzed with ancestry groups (k).

Heterogeneity (I²)

I² shows the proportion of variance in the effect estimate attributable to true between-ancestry differences vs. sampling error. This shows the distribution of I² values for phecode-variant pairs with k > 2 (I² is unreliable at k=2 (wide CI, bimodal distribution near 0% or 100%).

• High I² suggests ancestries are telling different stories.
• Low I² (near 0) is what we expect in a trans-ancestry PheWAS because most phecodes should show no association in any ancestry.

plot_ly(
  meta_results |> filter(n_groups > 2),
  x     = ~i2,
  color = ~rsid,
  type  = "histogram",
  opacity   = 0.7,
  hoverinfo = "x+y"
) |>
  layout(
    barmode = "overlay",
    title   = "Heterogeneity Distribution by Variant",
    xaxis   = list(title = "I² (%)"),
    yaxis   = list(title = "Count"),
    legend = list(orientation = "h", x = 0.5, y = -0.2, xanchor = "center")
  )

REML Boundary

REML (Random Effects ) estimates τ² = 0 when the data doesn’t show between-ancestry variance. This plot shows the percentage of phecode-variant meta-analyses where the random effects model collapsed to a fixed-effects model after filtering out data where the number of ancestry groups is < 2.

• High percentage indicate low heterogeneity (the ancestries are effectively agreeing).
• Lower percentage indicates greater between-ancestry variation (the ancestries are effectively disagreeing).

boundary_df <- meta_results |>
  filter(n_groups > 2) |>
  group_by(rsid) |>
  summarise(
    n_total    = n(),
    n_boundary = sum(tau2_boundary, na.rm = TRUE),
    pct        = round(100 * n_boundary / n_total, 1),
    .groups    = "drop"
  )

plot_ly(
  boundary_df,
  x    = ~reorder(rsid, pct),
  y    = ~pct,
  type = "bar",
  text        = ~paste0(pct, "% (", n_boundary, "/", n_total, ")"),
  textposition = "none"
) |>
  layout(
    title  = list(text = "REML τ² = 0 Boundary Rate", 
                  y = .99),
    xaxis  = list(title = ""),
    yaxis  = list(title = "% of Meta-Analyses at Boundary", range = c(0, 100))
  )