C9 penetrance vs. Family History Burden

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.

library(readxl)
library(ggplot2)
df <- read_excel("~/Documents/family_data.xlsx", sheet = 1)
names(df)

df <- df %>%
  mutate(
    family_history = if_else(Number_Affected >= 2, "High", "Low"),
    family_history = factor(family_history, levels = c("Low", "High"))
  )

ggplot(df, aes(x = Number_Affected, y = Penetrance_score, color = Total_Family)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se = TRUE) +
  scale_color_viridis_c(name = "Total Family\nSize") +
  labs(
    x     = "Number of Affected Relatives",
    y     = "C9 Penetrance Score",
    title = "Penetrance vs. Family Burden\n(colored by Total Family Size)"
  ) +
  theme_minimal(base_size = 14)

lm1 <- lm(Penetrance_score ~ Number_Affected, data = df)
summary(lm1)

lm2 <- lm(Penetrance_score ~ Number_Affected + Total_Family, data = df)
summary(lm2)

ggplot(df, aes(x = family_history, y = Penetrance_score, fill = family_history)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  labs(
    x = "Family History Group",
    y = "C9 Penetrance Score",
    title = "Low vs. High Family History"
  ) +
  theme_classic(base_size = 13) +
  scale_fill_viridis_d()

library(readxl)
library(dplyr)

G3;
Attaching package: ‘dplyr’

gG3;The following objects are masked from ‘package:stats’:

    filter, lag

gG3;The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

g

library(ggplot2)
library(viridis)  # for the color scale

G3;Loading required package: viridisLite
g

# 1. Read in your data
df <- read_excel("~/Documents/family_data.xlsx", sheet = 1)

# 2. (Re)create the Low/High factor if you need it later
df <- df %>%
  mutate(
    family_history = if_else(Number_Affected >= 2, "High", "Low"),
    family_history = factor(family_history, levels = c("Low","High"))
  )

# 3. Continuous scatter + regression, coloring by Total_family
ggplot(df, aes(
    x     = Number_Affected,
    y     = Penetrance_score,
    color = Total_family
  )) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se = TRUE) +
  scale_color_viridis_c(name = "Total Family\nSize") +
  labs(
    x     = "Number of Affected Relatives",
    y     = "C9 Penetrance Score",
    title = "Family Burden vs. C9 Penetrance"
  ) +
  theme_minimal(base_size = 14)

# 4. Linear model: penetrance ~ Number_Affected
lm1 <- lm(Penetrance_score ~ Number_Affected, data = df)
summary(lm1)

Call:
lm(formula = Penetrance_score ~ Number_Affected, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.27106 -0.08102  0.01894  0.05402  0.34894 

Coefficients:
                Estimate Std. Error t value
(Intercept)      0.34091    0.06460   5.277
Number_Affected  0.07008    0.02874   2.438
                Pr(>|t|)    
(Intercept)     6.16e-05 ***
Number_Affected    0.026 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1515 on 17 degrees of freedom
Multiple R-squared:  0.2591,    Adjusted R-squared:  0.2155 
F-statistic: 5.946 on 1 and 17 DF,  p-value: 0.02602

# 5. Adjusted model with Total_family
lm2 <- lm(Penetrance_score ~ Number_Affected + Total_family, data = df)
summary(lm2)

Call:
lm(formula = Penetrance_score ~ Number_Affected + Total_family, 
    data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.27186 -0.07715  0.02024  0.06274  0.32648 

Coefficients:
                 Estimate Std. Error t value
(Intercept)     0.3328102  0.0757433   4.394
Number_Affected 0.0654414  0.0361390   1.811
Total_family    0.0006987  0.0031309   0.223
                Pr(>|t|)    
(Intercept)     0.000453 ***
Number_Affected 0.088985 .  
Total_family    0.826225    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1559 on 16 degrees of freedom
Multiple R-squared:  0.2614,    Adjusted R-squared:  0.1691 
F-statistic: 2.832 on 2 and 16 DF,  p-value: 0.08855

summary(lm1)

Call:
lm(formula = Penetrance_score ~ Number_Affected, data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.27106 -0.08102  0.01894  0.05402  0.34894 

Coefficients:
                Estimate Std. Error t value
(Intercept)      0.34091    0.06460   5.277
Number_Affected  0.07008    0.02874   2.438
                Pr(>|t|)    
(Intercept)     6.16e-05 ***
Number_Affected    0.026 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1515 on 17 degrees of freedom
Multiple R-squared:  0.2591,    Adjusted R-squared:  0.2155 
F-statistic: 5.946 on 1 and 17 DF,  p-value: 0.02602

t_res <- t.test(Penetrance_score ~ family_history, data = df)
print(t_res)

    Welch Two Sample t-test

data:  Penetrance_score by family_history
t = -2.371, df = 16.29, p-value = 0.03039
alternative hypothesis: true difference in means between group Low and group High is not equal to 0
95 percent confidence interval:
 -0.29811478 -0.01688522
sample estimates:
 mean in group Low mean in group High 
            0.3825             0.5400 

df %>%
  group_by(family_history) %>%
  summarise(
    mean = mean(Penetrance_score),
    sd   = sd(Penetrance_score)
  ) %>%
  ggplot(aes(x = family_history, y = mean, fill = family_history)) +
  geom_col(width = 0.5) +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2) +
  labs(
    x = "Family History",
    y = "Mean Penetrance ± SD",
    title = "Group Means with Standard Deviations"
  ) +
  theme_minimal() +
  theme(legend.position = "none")

library(dplyr)

df <- df %>%
  mutate(
    K_category = if_else(Penetrance_score >= 0.5, "HighK", "LowK"),
    K_category = factor(K_category, levels = c("LowK","HighK"))
  )


library(tidyr)

heat_data <- df %>%
  group_by(ClinicalSyndrome, NeuropathSubtype) %>%
  summarise(
    total   = n(),
    highK   = sum(K_category == "HighK"),
    propHigh = highK / total,
    .groups = "drop"
  )

G1;Error in `group_by()`:
! Must group by variables found in `.data`.
Column `ClinicalSyndrome` is not found.
Column `NeuropathSubtype` is not found.
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.
g