NCounter Plots

#---------------------------------------- Use filter data to plot scatter Plot from long data with R and p value

#install.packages("reshape2")
#install.packages("ggsignif")
#install.packages("ggpubr")

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.4.2

library(readr)

Warning: package 'readr' was built under R version 4.4.2

library(dplyr)

Warning: package 'dplyr' was built under R version 4.4.2

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(Hmisc)  # For Spearman's rank correlation

Warning: package 'Hmisc' was built under R version 4.4.2

Attaching package: 'Hmisc'

The following objects are masked from 'package:dplyr':

    src, summarize

The following objects are masked from 'package:base':

    format.pval, units

# Read the CSV file
library(readr)
nCounter_Vera <- read_csv("nCounter_Vera.csv")

Rows: 3075 Columns: 45
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (9): HHKID, Visit, Subject Id, Gender, Age Group, FC ID, Date V0, Date ...
dbl (25): Age, Date Diff, C2i_No, Code, ACE, ACKR2, ACKR3, ACKR4, ACOX1, ACS...
lgl (11): Inclusion, Family Cohort, ECG V0, Lab V0, Hair V1, DXA V1, Fit Bit...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

wide <- nCounter_Vera |> dplyr::select(-c(1,2,3,7:21,23,24)) |> 
  relocate(C2i_No, .before = Gender)|> 
      dplyr::rename(Sex = Gender, 
                    SampleID = C2i_No)
wide 

# A tibble: 3,075 × 25
   SampleID Sex      Age `Age Group`   ACE ACKR2 ACKR3 ACKR4 ACOX1  ACSL1 ACSL3
      <dbl> <chr>  <dbl> <chr>       <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
 1  1200011 Female  31.6 30yr_39yr    5.87  6.99 21.5  11.5  1070.  7083.  87.5
 2  1200021 Female  29.6 20yr_29yr    3.26 10.2   8.48  1.74 1030.  4622.  64.2
 3  1200031 Female  75.2 70yr_80yr    2.51  2.51  2.51  2.51 1190.  9425.  44.5
 4  1200041 Female  75.9 70yr_80yr   14.0   2.38  2.38  4.47 1316.  7187.  49.8
 5  1200051 Male    57.6 50yr_59yr   26.1   1.28 45.3   3.04 1614.  6185.  91.4
 6  1200061 Female  60.1 60yr_69yr   19.9  13.7  30.7   5.98  938.  5056.  58.4
 7  1200071 Male    71.9 70yr_80yr   18.7  20.5  31.3  13.3  1055.  4894.  56.6
 8  1200081 Female  53.9 50yr_59yr   11.8  14.2  31.5   8.38  676.  5297.  63.9
 9  1200091 Female  52.8 50yr_59yr    3.41  6.13 18.4   6.13 1677. 11746.  68.8
10  1200101 Female  48.7 40yr_49yr    7.35 16.6  19.7  12.0  1433. 14858.  63.0
# ℹ 3,065 more rows
# ℹ 14 more variables: ACSL4 <dbl>, ACVR1 <dbl>, ADAR <dbl>, ADGRE5 <dbl>,
#   ADGRG3 <dbl>, ADORA2A <dbl>, AGT <dbl>, AHR <dbl>, AIF1 <dbl>, AIM2 <dbl>,
#   AKT1 <dbl>, AKT2 <dbl>, AKT3 <dbl>, ALOX12 <dbl>

library(tidyr)
long <- wide |> pivot_longer(
        cols = c(5:25),
        names_to = "Genes",
        values_to = "Measurements",
        values_drop_na = TRUE
      )
long

# A tibble: 64,575 × 6
   SampleID Sex      Age `Age Group` Genes Measurements
      <dbl> <chr>  <dbl> <chr>       <chr>        <dbl>
 1  1200011 Female  31.6 30yr_39yr   ACE           5.87
 2  1200011 Female  31.6 30yr_39yr   ACKR2         6.99
 3  1200011 Female  31.6 30yr_39yr   ACKR3        21.5 
 4  1200011 Female  31.6 30yr_39yr   ACKR4        11.5 
 5  1200011 Female  31.6 30yr_39yr   ACOX1      1070.  
 6  1200011 Female  31.6 30yr_39yr   ACSL1      7083.  
 7  1200011 Female  31.6 30yr_39yr   ACSL3        87.5 
 8  1200011 Female  31.6 30yr_39yr   ACSL4       707.  
 9  1200011 Female  31.6 30yr_39yr   ACVR1       264.  
10  1200011 Female  31.6 30yr_39yr   ADAR       1302.  
# ℹ 64,565 more rows

setwd("C:/Users/admin/Documents/lumx")
LongAnnotation <- read_csv("LongAnnotation.csv")

Rows: 2073 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): Genes, Pathway

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

combined_data <- long %>%
    left_join(LongAnnotation, by = "Genes")

Warning in left_join(., LongAnnotation, by = "Genes"): Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 10 of `x` matches multiple rows in `y`.
ℹ Row 8 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship =
  "many-to-many"` to silence this warning.

combined_data

# A tibble: 104,550 × 7
   SampleID Sex      Age `Age Group` Genes Measurements Pathway            
      <dbl> <chr>  <dbl> <chr>       <chr>        <dbl> <chr>              
 1  1200011 Female  31.6 30yr_39yr   ACE           5.87 Angiotensin_System 
 2  1200011 Female  31.6 30yr_39yr   ACKR2         6.99 Chemokine_Signaling
 3  1200011 Female  31.6 30yr_39yr   ACKR3        21.5  Chemokine_Signaling
 4  1200011 Female  31.6 30yr_39yr   ACKR4        11.5  Chemokine_Signaling
 5  1200011 Female  31.6 30yr_39yr   ACOX1      1070.   PPAR_Signaling     
 6  1200011 Female  31.6 30yr_39yr   ACSL1      7083.   PPAR_Signaling     
 7  1200011 Female  31.6 30yr_39yr   ACSL3        87.5  PPAR_Signaling     
 8  1200011 Female  31.6 30yr_39yr   ACSL4       707.   PPAR_Signaling     
 9  1200011 Female  31.6 30yr_39yr   ACVR1       264.   TGF_beta_Signaling 
10  1200011 Female  31.6 30yr_39yr   ADAR       1302.   DNA_Sensing        
# ℹ 104,540 more rows

Boxplot By Pathway

#Install the ggpubr package (if not already installed):
#install.packages("ggpubr")
#Load the ggpubr package:
library(ggpubr)

Warning: package 'ggpubr' was built under R version 4.4.2

library(dplyr)



# Create the original boxplot with adjusted y-axis
p1 <- ggplot(combined_data, aes(x = Pathway, y = Measurements)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6)) +
    labs(title = "Boxplot of Measurements by Pathways") +
    ylim(min(combined_data$Measurements) - 1, max(combined_data$Measurements) + 1)

p1

#pdf(file = "p1.pdf")
#p1
#dev.off()

#tried on add t-test/significance 

# Log-transform the measurements
combined_data <- combined_data %>%
    mutate(Log_Measurements = log(Measurements + 1))  # Adding 1 to avoid log(0)

# Create the boxplot for log-transformed measurements with adjusted y-axis
p2 <- ggplot(combined_data, aes(x = Pathway, y = Log_Measurements)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6))  +
    labs(title = "Boxplot of Log-Transformed Measurements by Pathways") +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1)

p2

pdf(file = "p2.pdf")
p2
dev.off()

png 
  2 

# Perform t-test by pathways comparing genders
ttest_results <- combined_data %>%
    group_by(Pathway) %>%
    summarise(
        t_stat = t.test(Measurements ~ Sex)$statistic,
        p_val = t.test(Measurements ~ Sex)$p.value
    )

ttest_results

# A tibble: 19 × 3
   Pathway                                    t_stat    p_val
   <chr>                                       <dbl>    <dbl>
 1 Angiotensin_System                         -1.09  2.74e- 1
 2 BCR_Signaling                              -1.41  1.59e- 1
 3 Chemokine_Signaling                        -4.04  5.44e- 5
 4 DNA_Sensing                                 4.89  1.03e- 6
 5 Inflammasomes                              11.5   6.29e-30
 6 Leukotriene_and_Prostaglandin_Inflammation -1.15  2.50e- 1
 7 Mononuclear_Cell_Migration                 -8.87  1.25e-18
 8 Myeloid_Activation                         -5.23  1.76e- 7
 9 Myeloid_Inflammation                        4.69  2.84e- 6
10 NLR_Signaling                              11.5   6.29e-30
11 NO_Signaling                               -0.958 3.38e- 1
12 Oxidative_Stress_Response                  -8.87  1.25e-18
13 PPAR_Signaling                             -0.278 7.81e- 1
14 TCR_Signaling                              -1.41  1.59e- 1
15 TGF_beta_Signaling                         10.2   3.90e-24
16 TH17_Differentiation                       -6.29  3.76e-10
17 TNF_Signaling                              -1.41  1.59e- 1
18 T_cell_Costimulation                        4.69  2.84e- 6
19 Type_I_Interferon_Signaling                 6.25  4.71e-10

# Create the boxplot with adjusted y-axis comparing genders
p3 <- ggplot(combined_data, aes(x = Pathway, y = Measurements, fill = Sex)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6))  +
    labs(title = "Boxplot of Measurements by Pathways Comparing Sex") +
    ylim(min(combined_data$Measurements) - 1, max(combined_data$Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))

p3

# Log-transform the measurements
combined_data <- combined_data %>%
    mutate(Log_Measurements = log(Measurements + 1))  # Adding 1 to avoid log(0)

# Create the boxplot with adjusted y-axis comparing genders
p4 <- ggplot(combined_data, aes(x = Pathway, y = Log_Measurements, fill = Sex)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6)) +
    labs(title = "Boxplot of Log-Transformed Measurements by Pathways Comparing Sex") +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))

p4

library(dplyr)

# Get unique pathways
unique_pathway <- unique(combined_data$Pathway)

combined_data$`Age Group` <- as.factor(combined_data$`Age Group`)

# Loop through each pathway and create a boxplot
for (pathway in unique_pathway) {
  pathway_data <- combined_data %>% filter(Pathway == pathway)
  
  p6 <- ggplot(pathway_data, aes(x = `Age Group`, y = Log_Measurements, fill = Sex)) +
   geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = paste("Log10 Measurements for Pathway:", pathway)) +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))
  
  print(p6)
}

Boxplot by Genes

# Create the original boxplot with adjusted y-axis grouped by genes
library(ggplot2)

ggplot2::ggplot(combined_data, aes(x = Genes, y = Measurements)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = "Boxplot of Measurements by Genes") +
    ylim(min(combined_data$Measurements) - 1, max(combined_data$Measurements) + 1)

library(ggplot2)
library(dplyr)
# Log-transform the measurements
combined_data <- combined_data %>%
    mutate(Log_Measurements = log(Measurements + 1))  # Adding 1 to avoid log(0)

# Create the boxplot for log-transformed measurements with adjusted y-axis grouped by genes
ggplot(combined_data, aes(x = Genes, y = Log_Measurements)) +
   geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = "Boxplot of Log-Transformed Measurements by Genes") +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1) 

# Create the boxplot with adjusted y-axis comparing genders
library(ggplot2)
ggplot(combined_data, aes(x = Genes, y = Measurements, fill = Sex)) +
   geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = "Boxplot of Measurements by Genes Comparing Sex") +
    ylim(min(combined_data$Measurements) - 1, max(combined_data$Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))

library(dplyr)
# Log-transform the measurements
combined_data <- combined_data %>%
    mutate(Log_Measurements = log(Measurements + 1))  # Adding 1 to avoid log(0)

# Perform t-test by pathways comparing genders
ttest_results <- combined_data %>%
    group_by(Genes) %>%
    summarise(
        t_stat = t.test(Log_Measurements ~ Sex)$statistic,
        p_val = t.test(Log_Measurements ~ Sex)$p.value
    )
ttest_results

# A tibble: 21 × 3
   Genes  t_stat    p_val
   <chr>   <dbl>    <dbl>
 1 ACE   -0.630  5.29e- 1
 2 ACKR2  2.39   1.71e- 2
 3 ACKR3 -7.07   1.84e-12
 4 ACKR4  1.55   1.21e- 1
 5 ACOX1 -4.52   6.44e- 6
 6 ACSL1 -0.0530 9.58e- 1
 7 ACSL3  0.473  6.37e- 1
 8 ACSL4  2.53   1.15e- 2
 9 ACVR1 11.3    3.51e-29
10 ADAR   9.18   5.95e-20
# ℹ 11 more rows

# Create the boxplot with adjusted y-axis comparing genders
ggplot(combined_data, aes(x = Genes, y = Log_Measurements, fill = Sex)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = "Boxplot of Log-Transformed Measurements by Genes Comparing Sex") +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))

# Get unique genes
unique_genes <- unique(combined_data$Genes)

combined_data$`Age Group` <- as.factor(combined_data$`Age Group`)

# Loop through each gene and create a boxplot
for (gene in unique_genes) {
  gene_data <- combined_data %>% filter(Genes == gene)
  
  p <- ggplot(gene_data, aes(x = `Age Group`, y = Log_Measurements, fill = Sex)) +
    geom_boxplot(size = 0.05) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8)) +
    labs(title = paste("Boxplot of Log-Transformed Measurements for Gene:", gene)) +
    ylim(min(combined_data$Log_Measurements) - 1, max(combined_data$Log_Measurements) + 1) +
    scale_fill_manual(values = c("Male" = "blue", "Female" = "pink"))
  
  print(p)
}

Heatmaps

library(ggplot2)
library(reshape2)

Warning: package 'reshape2' was built under R version 4.4.2

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

heatmap_data <- combined_data %>%
    dcast(Pathway ~ Genes, value.var = "Measurements", fun.aggregate = mean, na.rm = TRUE)

heatmap_data_melt <- melt(heatmap_data, id.vars = "Pathway")

ggplot(heatmap_data_melt, aes(x = variable, y = Pathway, fill = value)) +
    geom_tile() +
    scale_fill_gradient(low = "pink", high = "black") +
    labs(title = "Heatmap of Measurements by Pathways and Genes",
         x = "Genes",
         y = "Pathway",
         fill = "Measurements") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6), 
          axis.title = element_text(size = 8),
          axis.text.y = element_text(size = 6), 
          legend.title = element_text(size = 7))

# Log-transform the measurements
combined_data <- combined_data %>%
    mutate(Log_Measurements = log(Measurements + 1))  # Adding 1 to avoid log(0)

# Create the heatmap with log-transformed measurements
heatmap_data <- combined_data %>%
    dcast(Pathway ~ Genes, value.var = "Log_Measurements", fun.aggregate = mean, na.rm = TRUE)

heatmap_data_melt <- melt(heatmap_data, id.vars = "Pathway")

ggplot(heatmap_data_melt, aes(x = variable, y = Pathway, fill = value)) +
    geom_tile() +
    scale_fill_gradient(low = "pink", high = "black") +
    labs(title = "Heatmap of Log-Transformed Measurements by Pathways and Genes",
         x = "Genes",
         y = "Pathway",
         fill = "Log Measurements") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6), 
          axis.title = element_text(size = 8),
          axis.text.y = element_text(size = 6), 
          legend.title = element_text(size = 7))

library(tidyr)
library(dplyr)
library(ggplot2)
library(reshape2)
combined_data

# A tibble: 104,550 × 8
   SampleID Sex      Age `Age Group` Genes Measurements Pathway Log_Measurements
      <dbl> <chr>  <dbl> <fct>       <chr>        <dbl> <chr>              <dbl>
 1  1200011 Female  31.6 30yr_39yr   ACE           5.87 Angiot…             1.93
 2  1200011 Female  31.6 30yr_39yr   ACKR2         6.99 Chemok…             2.08
 3  1200011 Female  31.6 30yr_39yr   ACKR3        21.5  Chemok…             3.11
 4  1200011 Female  31.6 30yr_39yr   ACKR4        11.5  Chemok…             2.52
 5  1200011 Female  31.6 30yr_39yr   ACOX1      1070.   PPAR_S…             6.98
 6  1200011 Female  31.6 30yr_39yr   ACSL1      7083.   PPAR_S…             8.87
 7  1200011 Female  31.6 30yr_39yr   ACSL3        87.5  PPAR_S…             4.48
 8  1200011 Female  31.6 30yr_39yr   ACSL4       707.   PPAR_S…             6.56
 9  1200011 Female  31.6 30yr_39yr   ACVR1       264.   TGF_be…             5.58
10  1200011 Female  31.6 30yr_39yr   ADAR       1302.   DNA_Se…             7.17
# ℹ 104,540 more rows

# Convert Age_Group to a factor
combined_data$`Age Group` <- as.factor(combined_data$`Age Group`)

# Create a pivot table for the heatmap
heatmap_data <- combined_data %>%
    dcast(Pathway ~ `Age Group`, value.var = "Measurements", fun.aggregate = mean, na.rm = TRUE)

# Melt the data for ggplot2
heatmap_data_melt <- melt(heatmap_data, id.vars = "Pathway")

# Plot the heatmap
ggplot(heatmap_data_melt, aes(x = variable, y = Pathway, fill = value)) +
    geom_tile() +
    scale_fill_gradient(low = "white", high = "black") +
    labs(title = "Heatmap of Measurements by Age Group and Pathway",
         x = "Age Group",
         y = "Pathway",
         fill = "Measurements") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Log-transform the Measurements
combined_data$Log_Measurements <- log(combined_data$Measurements)

# Create a pivot table for the heatmap
heatmap_data <- combined_data %>%
    dcast(Pathway ~ `Age Group`, value.var = "Log_Measurements", fun.aggregate = mean, na.rm = TRUE)

# Melt the data for ggplot2
heatmap_data_melt <- melt(heatmap_data, id.vars = "Pathway")

# Plot the heatmap
ggplot(heatmap_data_melt, aes(x = variable, y = Pathway, fill = value)) +
    geom_tile() +
    scale_fill_gradient(low = "white", high = "black") +
    labs(title = "Heatmap of Log-Transformed Measurements by Age Group and Pathway",
         x = "Age Group",
         y = "Pathway",
         fill = "Log Measurements") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 6), 
          axis.title = element_text(size = 8),
          axis.text.y = element_text(size = 6), 
          legend.title = element_text(size = 7))

heatmap_data_melt

                                       Pathway  variable    value
1                           Angiotensin_System 20yr_29yr 2.346088
2                                BCR_Signaling 20yr_29yr 5.922350
3                          Chemokine_Signaling 20yr_29yr 2.131087
4                                  DNA_Sensing 20yr_29yr 6.294965
5                                Inflammasomes 20yr_29yr 5.503240
6   Leukotriene_and_Prostaglandin_Inflammation 20yr_29yr 1.106239
7                   Mononuclear_Cell_Migration 20yr_29yr 7.598491
8                           Myeloid_Activation 20yr_29yr 6.809164
9                         Myeloid_Inflammation 20yr_29yr 5.902315
10                               NLR_Signaling 20yr_29yr 5.503240
11                                NO_Signaling 20yr_29yr 2.398238
12                   Oxidative_Stress_Response 20yr_29yr 7.598491
13                              PPAR_Signaling 20yr_29yr 6.517714
14                        T_cell_Costimulation 20yr_29yr 5.902315
15                               TCR_Signaling 20yr_29yr 5.922350
16                          TGF_beta_Signaling 20yr_29yr 5.452984
17                        TH17_Differentiation 20yr_29yr 6.527679
18                               TNF_Signaling 20yr_29yr 5.922350
19                 Type_I_Interferon_Signaling 20yr_29yr 7.086691
20                          Angiotensin_System 30yr_39yr 2.412937
21                               BCR_Signaling 30yr_39yr 5.915640
22                         Chemokine_Signaling 30yr_39yr 2.197285
23                                 DNA_Sensing 30yr_39yr 6.349493
24                               Inflammasomes 30yr_39yr 5.571834
25  Leukotriene_and_Prostaglandin_Inflammation 30yr_39yr 1.206064
26                  Mononuclear_Cell_Migration 30yr_39yr 7.583406
27                          Myeloid_Activation 30yr_39yr 6.823015
28                        Myeloid_Inflammation 30yr_39yr 5.964276
29                               NLR_Signaling 30yr_39yr 5.571834
30                                NO_Signaling 30yr_39yr 2.464694
31                   Oxidative_Stress_Response 30yr_39yr 7.583406
32                              PPAR_Signaling 30yr_39yr 6.542996
33                        T_cell_Costimulation 30yr_39yr 5.964276
34                               TCR_Signaling 30yr_39yr 5.915640
35                          TGF_beta_Signaling 30yr_39yr 5.455047
36                        TH17_Differentiation 30yr_39yr 6.520399
37                               TNF_Signaling 30yr_39yr 5.915640
38                 Type_I_Interferon_Signaling 30yr_39yr 7.127151
39                          Angiotensin_System 40yr_49yr 2.362083
40                               BCR_Signaling 40yr_49yr 5.918711
41                         Chemokine_Signaling 40yr_49yr 2.150471
42                                 DNA_Sensing 40yr_49yr 6.333122
43                               Inflammasomes 40yr_49yr 5.565584
44  Leukotriene_and_Prostaglandin_Inflammation 40yr_49yr 1.172400
45                  Mononuclear_Cell_Migration 40yr_49yr 7.564925
46                          Myeloid_Activation 40yr_49yr 6.804775
47                        Myeloid_Inflammation 40yr_49yr 5.989927
48                               NLR_Signaling 40yr_49yr 5.565584
49                                NO_Signaling 40yr_49yr 2.449727
50                   Oxidative_Stress_Response 40yr_49yr 7.564925
51                              PPAR_Signaling 40yr_49yr 6.522277
52                        T_cell_Costimulation 40yr_49yr 5.989927
53                               TCR_Signaling 40yr_49yr 5.918711
54                          TGF_beta_Signaling 40yr_49yr 5.452005
55                        TH17_Differentiation 40yr_49yr 6.493111
56                               TNF_Signaling 40yr_49yr 5.918711
57                 Type_I_Interferon_Signaling 40yr_49yr 7.100660
58                          Angiotensin_System 50yr_59yr 2.356228
59                               BCR_Signaling 50yr_59yr 5.924444
60                         Chemokine_Signaling 50yr_59yr 2.232268
61                                 DNA_Sensing 50yr_59yr 6.316104
62                               Inflammasomes 50yr_59yr 5.522424
63  Leukotriene_and_Prostaglandin_Inflammation 50yr_59yr 1.247175
64                  Mononuclear_Cell_Migration 50yr_59yr 7.568239
65                          Myeloid_Activation 50yr_59yr 6.799370
66                        Myeloid_Inflammation 50yr_59yr 5.984336
67                               NLR_Signaling 50yr_59yr 5.522424
68                                NO_Signaling 50yr_59yr 2.420011
69                   Oxidative_Stress_Response 50yr_59yr 7.568239
70                              PPAR_Signaling 50yr_59yr 6.514966
71                        T_cell_Costimulation 50yr_59yr 5.984336
72                               TCR_Signaling 50yr_59yr 5.924444
73                          TGF_beta_Signaling 50yr_59yr 5.444069
74                        TH17_Differentiation 50yr_59yr 6.511957
75                               TNF_Signaling 50yr_59yr 5.924444
76                 Type_I_Interferon_Signaling 50yr_59yr 7.109784
77                          Angiotensin_System 60yr_69yr 2.361770
78                               BCR_Signaling 60yr_69yr 5.918631
79                         Chemokine_Signaling 60yr_69yr 2.300007
80                                 DNA_Sensing 60yr_69yr 6.345105
81                               Inflammasomes 60yr_69yr 5.590698
82  Leukotriene_and_Prostaglandin_Inflammation 60yr_69yr 1.215685
83                  Mononuclear_Cell_Migration 60yr_69yr 7.538955
84                          Myeloid_Activation 60yr_69yr 6.792548
85                        Myeloid_Inflammation 60yr_69yr 5.974649
86                               NLR_Signaling 60yr_69yr 5.590698
87                                NO_Signaling 60yr_69yr 2.465648
88                   Oxidative_Stress_Response 60yr_69yr 7.538955
89                              PPAR_Signaling 60yr_69yr 6.531396
90                        T_cell_Costimulation 60yr_69yr 5.974649
91                               TCR_Signaling 60yr_69yr 5.918631
92                          TGF_beta_Signaling 60yr_69yr 5.467647
93                        TH17_Differentiation 60yr_69yr 6.520051
94                               TNF_Signaling 60yr_69yr 5.918631
95                 Type_I_Interferon_Signaling 60yr_69yr 7.099512
96                          Angiotensin_System 70yr_80yr 2.419829
97                               BCR_Signaling 70yr_80yr 5.927711
98                         Chemokine_Signaling 70yr_80yr 2.275514
99                                 DNA_Sensing 70yr_80yr 6.391167
100                              Inflammasomes 70yr_80yr 5.662168
101 Leukotriene_and_Prostaglandin_Inflammation 70yr_80yr 1.269287
102                 Mononuclear_Cell_Migration 70yr_80yr 7.547959
103                         Myeloid_Activation 70yr_80yr 6.809017
104                       Myeloid_Inflammation 70yr_80yr 5.969582
105                              NLR_Signaling 70yr_80yr 5.662168
106                               NO_Signaling 70yr_80yr 2.503028
107                  Oxidative_Stress_Response 70yr_80yr 7.547959
108                             PPAR_Signaling 70yr_80yr 6.559378
109                       T_cell_Costimulation 70yr_80yr 5.969582
110                              TCR_Signaling 70yr_80yr 5.927711
111                         TGF_beta_Signaling 70yr_80yr 5.468173
112                       TH17_Differentiation 70yr_80yr 6.507793
113                              TNF_Signaling 70yr_80yr 5.927711
114                Type_I_Interferon_Signaling 70yr_80yr 7.120165

library(ggplot2)
library(ggpubr)
library(dplyr)
# Calculate Spearman's Rank Correlation and p-value by pathway
spearman_results <- combined_data %>%
    group_by(Pathway) %>%
    summarise(
        spearman_corr = cor(Age, Measurements, method = "spearman"),
        p_val = tryCatch(cor.test(Age, Measurements, method = "spearman")$p.value, warning = function(w) NA)
    )

# Create the regression plot with smaller legend key size
ggplot(combined_data, aes(x = Age, y = Measurements, color = Pathway)) +
    geom_point(size = 1, alpha = 0.6) +  # Increase point size and add transparency
    geom_smooth(method = "lm", se = FALSE, linetype = "dashed", size = 1) +  # Customize smoothing line
    stat_cor(method = "spearman", label.x.npc = "left", label.y.npc = "top", size = 1.5) +   # Adjust correlation label size
    labs(title = "Regression Plot with Spearman's Rank Correlation by Pathway",
         x = "Age",
         y = "Measurements",
         color = "Pathway") +
    theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.2, size = 8),  # Adjust axis text size
          axis.title = element_text(size = 10),
          axis.text.y = element_text(size = 8), 
          legend.title = element_text(size = 9),
          legend.key.size = unit(0.055, 'cm'),  # Adjust legend key size
          plot.title = element_text(size = 12, face = "bold"))  # Adjust plot title size and style

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

`geom_smooth()` using formula = 'y ~ x'

#Load the necessary libraries:

library(ggplot2)
library(dplyr)
#Read in your datasets (assuming you have already joined them as combined_data):



combined_data$Age <- as.numeric(combined_data$Age)
#Loop through each pathway and create regression plots:

pathways <- unique(combined_data$Pathway)

for (pathway in pathways) {
    pathway_data <- combined_data %>% filter(Pathway == pathway)
    
    # Calculate Spearman's Rank Correlation and p-value, suppressing warnings
    spearman_test <- suppressWarnings(cor.test(pathway_data$Age, pathway_data$Measurements, method = "spearman"))
    spearman_corr <- spearman_test$estimate
    p_value <- spearman_test$p.value
    
    # Create regression plot with sex information
    plot <- ggplot(pathway_data, aes(x = Age, y = Measurements, color = Sex)) +
        geom_point(size = 0.5) +
        geom_smooth(method = "lm", se = FALSE, size = 1.3) +
        annotate("text", x = Inf, y = Inf, label = paste("Spearman's rho:", round(spearman_corr, 2), "\np-value:", round(p_value, 4)), hjust = 1.1, vjust = 1.1, size = 3, color = "red") +
        scale_color_manual(values = c("Male" = "dodgerblue", "Female" = "firebrick2")) +
        labs(title = paste("Regression Plot for Pathway:", pathway),
             x = "Age",
             y = "Measurements") +
        theme_minimal()
    
    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

combined_data$Age <- as.numeric(combined_data$Age)
combined_data$Log_Measurements <- log(combined_data$Measurements)

# Loop through each pathway and create regression plots
pathways <- unique(combined_data$Pathway)

for (pathway in pathways) {
    pathway_data <- combined_data %>% filter(Pathway == pathway)
    
    # Calculate Spearman's Rank Correlation and p-value, suppressing warnings
    spearman_test <- suppressWarnings(cor.test(pathway_data$Age, pathway_data$Log_Measurements, method = "spearman"))
    spearman_corr <- spearman_test$estimate
    p_value <- spearman_test$p.value
    
    # Create regression plot with sex information
    plot <- ggplot(pathway_data, aes(x = Age, y = Log_Measurements, color = Sex)) +
        geom_point(size = 0.5) +
        geom_smooth(method = "lm", se = FALSE, size = 1.3) +
        annotate("text", x = Inf, y = Inf, label = paste("Spearman's rho:", round(spearman_corr, 2), "\np-value:", round(p_value, 4)), hjust = 1.1, vjust = 1.1, size = 3, color = "red") +
        scale_color_manual(values = c("Male" = "dodgerblue", "Female" = "firebrick2")) +
        labs(title = paste("Regression Plot for Pathway:", pathway),
             x = "Age",
             y = "Log Measurements") +
        theme_minimal()
    
    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

genes <- unique(combined_data$Genes)

for (gene in genes) {
    gene_data <- combined_data %>% filter(Genes == gene)
    
    # Calculate Spearman's Rank Correlation and p-value, suppressing warnings
    spearman_test <- suppressWarnings(cor.test(gene_data$Age, gene_data$Measurements, method = "spearman"))
    spearman_corr <- spearman_test$estimate
    p_value <- spearman_test$p.value
    
    # Create regression plot with sex information
    plot <- ggplot(gene_data, aes(x = Age, y = Measurements, color = Sex)) +
        geom_point(size = 0.5) +
        geom_smooth(method = "lm", se = FALSE, size = 1.3) +
        annotate("text", x = Inf, y = Inf, label = paste("Spearman's rho:", round(spearman_corr, 2), "\np-value:", round(p_value, 4)), hjust = 1.1, vjust = 1.1, size = 3, color = "red") +
        scale_color_manual(values = c("Male" = "dodgerblue", "Female" = "firebrick2")) +
        labs(title = paste("Regression Plot for Gene:", gene),
             x = "Age",
             y = "Measurements") +
        theme_minimal()
    
    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

combined_data$Age <- as.numeric(combined_data$Age)
combined_data$Log_Measurements <- log(combined_data$Measurements)

# Loop through each gene and create regression plots
genes <- unique(combined_data$Genes)

for (gene in genes) {
    gene_data <- combined_data %>% filter(Genes == gene)
    
    # Calculate Spearman's Rank Correlation and p-value, suppressing warnings
    spearman_test <- suppressWarnings(cor.test(gene_data$Age, gene_data$Log_Measurements, method = "spearman"))
    spearman_corr <- spearman_test$estimate
    p_value <- spearman_test$p.value
    
    # Create regression plot with sex information
    plot <- ggplot(gene_data, aes(x = Age, y = Log_Measurements, color = Sex)) +
        geom_point(size = 0.5) +
        geom_smooth(method = "lm", se = FALSE, size = 1.3) +
        annotate("text", x = Inf, y = Inf, label = paste("Spearman's rho:", round(spearman_corr, 2), "\np-value:", round(p_value, 4)), hjust = 1.1, vjust = 1.1, size = 3, color = "red") +
        scale_color_manual(values = c("Male" = "dodgerblue", "Female" = "firebrick2")) +
        labs(title = paste("Regression Plot for Gene:", gene),
             x = "Age",
             y = "Log Measurements") +
        theme_minimal()
    
    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

# Count the number of Female and Male for each pathway
count_female <- combined_data %>% filter(Sex == "Female") %>% nrow()
count_male <- combined_data %>% filter(Sex == "Male") %>% nrow()

# Calculate Spearman's rank correlation coefficient for Female and Male
female_data <- combined_data %>% filter(Sex == "Female")
male_data <- combined_data %>% filter(Sex == "Male")

for (pathway in pathways) {
    pathway_data <- combined_data %>% filter(Pathway == pathway)

    if (nrow(female_data) > 0) {
        spearman_female <- rcorr(female_data$Age, female_data$Measurements, type = "spearman")
        r_female <- spearman_female$r[1, 2]
        p_female <- spearman_female$P[1, 2]
    } else {
        r_female <- NA
        p_female <- NA
    }

    if (nrow(male_data) > 0) {
        spearman_male <- rcorr(male_data$Age, male_data$Measurements, type = "spearman")
        r_male <- spearman_male$r[1, 2]
        p_male <- spearman_male$P[1, 2]
    } else {
        r_male <- NA
        p_male <- NA
    }

    # Create the regression plot
    plot <- ggplot(pathway_data, aes(x = Age, y = Measurements, color = Sex)) +
        geom_point(position = position_jitter(width = 0.1), alpha = 0.5) +
        geom_smooth(method = "loess", se = FALSE) +
        labs(title = paste("Regression Plot of", pathway, "by Gender"),
             subtitle = paste("n (Female) =", count_female, ", n (Male) =", count_male),
             x = "Age",
             y = "Measurements in log10 scale") + 
        scale_y_log10() +
        scale_color_manual(values = c("Female" = "#ff7070", "Male" = "#6495ED")) +  # Set colors
        theme_bw() +
        annotate("text", x = max(pathway_data$Age, na.rm = TRUE) * 0.8, y = max(pathway_data$Measurements, na.rm = TRUE) * 1.1, 
                 label = paste("Female: R =", round(r_female, 2), ", p =", format(p_female, scientific = TRUE)), 
                 hjust = 0.5, vjust = 1, size = 3, color = "black") +
        annotate("text", x = max(pathway_data$Age, na.rm = TRUE) * 0.8, y = max(pathway_data$Measurements, na.rm = TRUE) * 1.05, 
                 label = paste("Male: R =", round(r_male, 2), ", p =", format(p_male, scientific = TRUE)), 
                 hjust = 0.5, vjust = 2, size = 3, color = "black")

    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

# Count the number of Female and Male for each pathway
count_female <- combined_data %>% filter(Sex == "Female") %>% nrow()
count_male <- combined_data %>% filter(Sex == "Male") %>% nrow()

# Calculate Spearman's rank correlation coefficient for Female and Male
female_data <- combined_data %>% filter(Sex == "Female")
male_data <- combined_data %>% filter(Sex == "Male")

for (gene in genes) {
    gene_data <- combined_data %>% filter(Genes == gene)

    if (nrow(female_data) > 0) {
        spearman_female <- rcorr(female_data$Age, female_data$Measurements, type = "spearman")
        r_female <- spearman_female$r[1, 2]
        p_female <- spearman_female$P[1, 2]
    } else {
        r_female <- NA
        p_female <- NA
    }

    if (nrow(male_data) > 0) {
        spearman_male <- rcorr(male_data$Age, male_data$Measurements, type = "spearman")
        r_male <- spearman_male$r[1, 2]
        p_male <- spearman_male$P[1, 2]
    } else {
        r_male <- NA
        p_male <- NA
    }

    # Create the regression plot
    plot <- ggplot(gene_data, aes(x = Age, y = Measurements, color = Sex)) +
        geom_point(position = position_jitter(width = 0.1), alpha = 0.5) +
        geom_smooth(method = "loess", se = FALSE) +
        labs(title = paste("Regression Plot of", gene, "by Gender"),
             subtitle = paste("n (Female) =", count_female, ", n (Male) =", count_male),
             x = "Age",
             y = "Measurements in log10 scale") + 
        scale_y_log10() +
        scale_color_manual(values = c("Female" = "#ff7070", "Male" = "#6495ED")) +  # Set colors
        theme_bw() +
        annotate("text", x = max(gene_data$Age, na.rm = TRUE) * 0.8, y = max(gene_data$Measurements, na.rm = TRUE) * 1.1, 
                 label = paste("Female: R =", round(r_female, 2), ", p =", format(p_female, scientific = TRUE)), 
                 hjust = 0.5, vjust = 1, size = 3, color = "black") +
        annotate("text", x = max(gene_data$Age, na.rm = TRUE) * 0.8, y = max(gene_data$Measurements, na.rm = TRUE) * 1.05, 
                 label = paste("Male: R =", round(r_male, 2), ", p =", format(p_male, scientific = TRUE)), 
                 hjust = 0.5, vjust = 2, size = 3, color = "black")

    print(plot)
}

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'