The TREM2 R47H variant is associated with liver-plasma-brain axis dyshomeostasis in the 5xFAD mouse model of Alzheimer's disease.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches

[1] § Methods › Statistical analysis ↔ TREM2-metabolomics-figures-final-v3.R, lines 4423–4472 · score 0.72 · Venn diagrams, ggVennDiagram, brain liver, brain plasma, Genes, heatmaps
[2] § Results › Effects of genotype along the liver-plasma-brain axis ↔ TREM2-metabolomics-figures-final-v3.R, lines 4423–4472 · score 0.56 · Venn diagrams, comparing brain liver, brain plasma, analytes, xFAD, Spearman

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 4,577 lines · 195 KB · MIT · 2 matches

rm(list=ls())
library(tidyverse)
library(readxl)
library(vegan)
library(pairwiseAdonis)
library(pheatmap)
library(ggpubr)
library(dplyr)
library(RColorBrewer)
library(FSA)
library(viridis)
#setwd("/media/gfaraci/model-AD/TREM2")
setwd("C:/Users/Gina Faraci/Documents/R/model-AD/TREM2 metabolome/From Julio")
#z-score function
z_score <- function(x){
(x - mean(x)) / sd(x)
}
#save_pheatmap function by mathzero: https://gist.github.com/mathzero/a2070a24a6b418740c44a5c023f5c01e
save_pheatmap <- function(x, filename, width=12, height=12){
stopifnot(!missing(x))
stopifnot(!missing(filename))
if(grepl(".png",filename)){
png(filename, width=width, height=height, units = "in", res=600)
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
else if(grepl(".pdf",filename)){
pdf(filename, width=width, height=height)
grid::grid.newpage()
grid::grid.draw(x$gtable)
dev.off()
}
else{
print("Filename did not contain '.png' or '.pdf'")
}
}
##############Preparing data##############
#Metadata
metadata = read.delim("metadata.txt")
metadata$Name = as.character(metadata$Name)
metadata$Age = as.character(metadata$Age)
#SampleIDs
SampleIDs = read_excel("u54-ad-sample-list.xlsx",
col_types = c("text", "text", "skip",
"skip", "skip", "skip"))
#Importing and wrangling data
plasma_lip_df = as.data.frame(t(read.csv("u54_ad_batch_corrected_data_TL_raw.csv", row.names = 1, check.names = F))) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "tm_id") %>%
inner_join(SampleIDs) %>%
select(!tm_id) %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_plasma"))
plasma_tar_df = as.data.frame(t(read.csv("u54_ad_batch_corrected_data_TM_raw.csv", row.names = 1, check.names = F))) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "tm_id") %>%
inner_join(SampleIDs) %>%
select(!tm_id) %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_plasma"))
liver_df = read_excel("20220519_mapstone_ad_model_liver_107_5500QTRAP_data_combined.xlsx", col_names = F)
names(liver_df) = liver_df[2,]
names(liver_df)[2:3] = c("Common name", "Mode")
liver_df = liver_df[-(1:2),]
liver_lip_df = liver_df %>%
filter(Mode == "Targeted lipidomics") %>%
select(!c(`Common name`, Mode)) %>%
column_to_rownames(var = "Sample IDs")
liver_lip_df = as.data.frame(t(liver_lip_df)) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "sample_id") %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_liver"))
liver_lip_df$value = as.numeric(liver_lip_df$value) #Abundance was character
liver_tar_df = liver_df %>%
filter(Mode == "Tarrgeted Metabolomics") %>%
select(!c(`Common name`, Mode)) %>%
column_to_rownames(var = "Sample IDs")
liver_tar_df = as.data.frame(t(liver_tar_df)) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "sample_id") %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_liver"))
liver_tar_df$value = as.numeric(liver_tar_df$value) #Abundance was character
brain_df = read_excel("20220519_mapstone_ad_model_cortex_159_5500QTRAP_combined_data.xlsx", col_names = F)
names(brain_df) = brain_df[2,]
brain_df = brain_df[-(1:2),]
brain_lip_df = brain_df %>%
select(!`9732`) %>% #This sample ID is duplicated for whatever reason
filter(Mode == "Targeted lipidomics") %>%
select(!c(`Sample IDs`, Mode)) %>%
column_to_rownames(var = "Metabolites")
brain_lip_df = as.data.frame(t(brain_lip_df)) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "sample_id") %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_brain"))
brain_tar_df = brain_df %>%
select(!`9732`) %>% #This sample ID is duplicated for whatever reason
filter(Mode == "Targeted Metabolomics") %>%
select(!c(`Sample IDs`, Mode)) %>%
column_to_rownames(var = "Metabolites")
brain_tar_df = as.data.frame(t(brain_tar_df)) %>%
replace(is.na(.), 0) %>%
rownames_to_column(var = "sample_id") %>%
pivot_longer(!sample_id) %>%
mutate(sample_id = paste0(sample_id, "_brain"))
#Combining dataframes
#Keeping targeted metabolomics and lipidomics seperate.
lip_df_combined = bind_rows(plasma_lip_df, brain_lip_df, liver_lip_df) %>%
pivot_wider(id_cols = "sample_id", names_from = "name", values_from = "value") %>%
separate_wider_delim(cols = sample_id, names = c("Name", "Mode"), delim = "_", cols_remove = F) %>%
inner_join(metadata) %>%
filter(!sample_id %in% c("9379_plasma", "7772_plasma", "20740_brain"), !Age == "18") %>%
relocate(CageID:Background, .after = sample_id)
lip_df_combined[,-(1:12)] = lip_df_combined[,-(1:12)] %>% #No normalization
replace(is.na(.), 0)
#Use the following lines instead for PERMANOVAs and PCAs
# lip_df_combined[,-(1:12)] = log10(lip_df_combined[,-(1:12)]) %>% #Normalization
# replace(is.na(.), 0)
#Export for future microbiome studies
write.table(lip_df_combined, "./lip_df_combined.txt", sep = "\t", row.names = TRUE, quote=FALSE)
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
s.test <- as.data.frame(sapply(lip_df_combined[,-(1:12)], as.numeric))
shapiro_results <- apply(s.test, 2, shapiro.test)
s.list <- numeric(length(967))
for (i in 1:967){
if (shapiro_results[[i]]$p.value > 0.05){
s.list[[i]] <- 1
} else {
s.list[[i]] <- 0
}
}
# data is NOT normally distributed
tar_df_combined = bind_rows(plasma_tar_df, brain_tar_df, liver_tar_df) %>%
pivot_wider(id_cols = "sample_id", names_from = "name", values_from = "value") %>%
separate_wider_delim(cols = sample_id, names = c("Name", "Mode"), delim = "_", cols_remove = F) %>%
inner_join(metadata) %>%
filter(!Age == "18") %>%
relocate(CageID:Background, .after = sample_id)
tar_df_combined[,-(1:12)] = tar_df_combined[,-(1:12)] %>% #No normalization
replace(is.na(.), 0)
#Use the following lines instead for PERMANOVAs and PCAs
# tar_df_combined[,-(1:12)] = log10(tar_df_combined[,-(1:12)]) %>% #Normalization
# replace(is.na(.), 0)
#Export for future microbiome studies
write.table(tar_df_combined, "./tar_df_combined.txt", sep = "\t", row.names = TRUE, quote=FALSE)
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
s.test <- as.data.frame(sapply(tar_df_combined[,-(1:12)], as.numeric))
shapiro_results <- apply(s.test, 2, shapiro.test)
s.list <- numeric(length(324))
for (i in 1:324){
if (shapiro_results[[i]]$p.value > 0.05){
s.list[[i]] <- 1
} else {
s.list[[i]] <- 0
}
}
# data is NOT normally distributed
##############
##############Metadata analysis (misc)##############
##Count sample numbers (Table 1) AND
##Are there sig. differences in mouse weight between genotypes?
metadata_wt = read_excel("metadata-weights.xlsx")
weight_df = metadata_wt[metadata_wt$Age != 18, ] #182 samples
wt_counts <- weight_df[weight_df$Group2 == "WT 5xFAD", ] #65 WT
wt_4mo_counts <- wt_counts[wt_counts$Age == "4", ] #27 4mo
wt_4mo_f_counts <- wt_4mo_counts[wt_4mo_counts$Sex == "F", ] #13 females
wt_4mo_m_counts <- wt_4mo_counts[wt_4mo_counts$Sex == "M", ] #14 males
wt_12mo_counts <- wt_counts[wt_counts$Age == "12", ] #38 12mo
wt_12mo_f_counts <- wt_12mo_counts[wt_12mo_counts$Sex == "F", ] #17 females
wt_12mo_m_counts <- wt_12mo_counts[wt_12mo_counts$Sex == "M", ] #21 males
fad_counts <- weight_df[weight_df$Group2 == "HEMI 5xFAD", ] #65 5xFAD
#note, 20737 and 20735 we did not collect data for (both 12mo females)
fad_4mo_counts <- fad_counts[fad_counts$Age == "4", ] #28 4mo
fad_4mo_f_counts <- fad_4mo_counts[fad_4mo_counts$Sex == "F", ] #14 females
fad_4mo_m_counts <- fad_4mo_counts[fad_4mo_counts$Sex == "M", ] #14 males
fad_12mo_counts <- fad_counts[fad_counts$Age == "12", ] #37->35 12mo
fad_12mo_f_counts <- fad_12mo_counts[fad_12mo_counts$Sex == "F", ] #19->17 females
fad_12mo_m_counts <- fad_12mo_counts[fad_12mo_counts$Sex == "M", ] #18 males
t2_counts <- weight_df[weight_df$Group2 == "WT 5xFAD_TREM2", ] #27 TREM2
t2_4mo_counts <- t2_counts[t2_counts$Age == "4", ] #10 4mo
t2_4mo_f_counts <- t2_4mo_counts[t2_4mo_counts$Sex == "F", ] #5 females
t2_4mo_m_counts <- t2_4mo_counts[t2_4mo_counts$Sex == "M", ] #5 males
t2_12mo_counts <- t2_counts[t2_counts$Age == "12", ] #17 12mo
t2_12mo_f_counts <- t2_12mo_counts[t2_12mo_counts$Sex == "F", ] #8 females
t2_12mo_m_counts <- t2_12mo_counts[t2_12mo_counts$Sex == "M", ] #9 males
fadt2_counts <- weight_df[weight_df$Group2 == "HEMI 5xFAD_TREM2", ] #25 5xFAD, TREM2
fadt2_4mo_counts <- fadt2_counts[fadt2_counts$Age == "4", ] #10 4mo
fadt2_4mo_f_counts <- fadt2_4mo_counts[fadt2_4mo_counts$Sex == "F", ] #5 females
fadt2_4mo_m_counts <- fadt2_4mo_counts[fadt2_4mo_counts$Sex == "M", ] #5 males
fadt2_12mo_counts <- fadt2_counts[fadt2_counts$Age == "12", ] #15 12mo
fadt2_12mo_f_counts <- fadt2_12mo_counts[fadt2_12mo_counts$Sex == "F", ] #8 females
fadt2_12mo_m_counts <- fadt2_12mo_counts[fadt2_12mo_counts$Sex == "M", ] #7 males
plasma_counts <- tar_df_combined[tar_df_combined$Mode == "plasma", ] #164 plasma samples
wt_plasma_counts <- plasma_counts[plasma_counts$Group2 == "WT 5xFAD", ] #58 WT samples
wt_4mo_plasma_counts <- wt_plasma_counts[wt_plasma_counts$Age == "4", ] #25 4mo
wt_4mo_f_plasma_counts <- wt_4mo_plasma_counts[wt_4mo_plasma_counts$Sex == "F", ] #12 females
wt_4mo_m_plasma_counts <- wt_4mo_plasma_counts[wt_4mo_plasma_counts$Sex == "M", ] #13 males
wt_12mo_plasma_counts <- wt_plasma_counts[wt_plasma_counts$Age == "12", ] #33 12mo
wt_12mo_f_plasma_counts <- wt_12mo_plasma_counts[wt_12mo_plasma_counts$Sex == "F", ] #15 females
wt_12mo_m_plasma_counts <- wt_12mo_plasma_counts[wt_12mo_plasma_counts$Sex == "M", ] #18 males
fad_plasma_counts <- plasma_counts[plasma_counts$Group2 == "HEMI 5xFAD", ] #56 5xFAD samples
fad_4mo_plasma_counts <- fad_plasma_counts[fad_plasma_counts$Age == "4", ] #26 4mo
fad_4mo_f_plasma_counts <- fad_4mo_plasma_counts[fad_4mo_plasma_counts$Sex == "F", ] #13 females
fad_4mo_m_plasma_counts <- fad_4mo_plasma_counts[fad_4mo_plasma_counts$Sex == "M", ] #13 males
fad_12mo_plasma_counts <- fad_plasma_counts[fad_plasma_counts$Age == "12", ] #30 12mo
fad_12mo_f_plasma_counts <- fad_12mo_plasma_counts[fad_12mo_plasma_counts$Sex == "F", ] #15 females
fad_12mo_m_plasma_counts <- fad_12mo_plasma_counts[fad_12mo_plasma_counts$Sex == "M", ] #15 males
t2_plasma_counts <- plasma_counts[plasma_counts$Group2 == "WT 5xFAD_TREM2", ] #27 TREM2 samples
t2_4mo_plasma_counts <- t2_plasma_counts[t2_plasma_counts$Age == "4", ] #10 4mo
t2_4mo_f_plasma_counts <- t2_4mo_plasma_counts[t2_4mo_plasma_counts$Sex == "F", ] #5 females
t2_4mo_m_plasma_counts <- t2_4mo_plasma_counts[t2_4mo_plasma_counts$Sex == "M", ] #5 males
t2_12mo_plasma_counts <- t2_plasma_counts[t2_plasma_counts$Age == "12", ] #17 12mo
t2_12mo_f_plasma_counts <- t2_12mo_plasma_counts[t2_12mo_plasma_counts$Sex == "F", ] #8 females
t2_12mo_m_plasma_counts <- t2_12mo_plasma_counts[t2_12mo_plasma_counts$Sex == "M", ] #9 males
fadt2_plasma_counts <- plasma_counts[plasma_counts$Group2 == "HEMI 5xFAD_TREM2", ] #23 5xFAD, TREM2 samples
fadt2_4mo_plasma_counts <- fadt2_plasma_counts[fadt2_plasma_counts$Age == "4", ] #10 4mo
fadt2_4mo_f_plasma_counts <- fadt2_4mo_plasma_counts[fadt2_4mo_plasma_counts$Sex == "F", ] #5 females
fadt2_4mo_m_plasma_counts <- fadt2_4mo_plasma_counts[fadt2_4mo_plasma_counts$Sex == "M", ] #5 males
fadt2_12mo_plasma_counts <- fadt2_plasma_counts[fadt2_plasma_counts$Age == "12", ] #13 12mo
fadt2_12mo_f_plasma_counts <- fadt2_12mo_plasma_counts[fadt2_12mo_plasma_counts$Sex == "F", ] #6 females
fadt2_12mo_m_plasma_counts <- fadt2_12mo_plasma_counts[fadt2_12mo_plasma_counts$Sex == "M", ] #7 males
liver_counts <- tar_df_combined[tar_df_combined$Mode == "liver", ] #105 liver samples
wt_liver_counts <- liver_counts[liver_counts$Group2 == "WT 5xFAD", ] #34 WT samples
wt_4mo_liver_counts <- wt_liver_counts[wt_liver_counts$Age == "4", ] #17 4mo
wt_4mo_f_liver_counts <- wt_4mo_liver_counts[wt_4mo_liver_counts$Sex == "F", ] #8 females
wt_4mo_m_liver_counts <- wt_4mo_liver_counts[wt_4mo_liver_counts$Sex == "M", ] #9 males
wt_12mo_liver_counts <- wt_liver_counts[wt_liver_counts$Age == "12", ] #17 12mo
wt_12mo_f_liver_counts <- wt_12mo_liver_counts[wt_12mo_liver_counts$Sex == "F", ] #7 females
wt_12mo_m_liver_counts <- wt_12mo_liver_counts[wt_12mo_liver_counts$Sex == "M", ] #10 males
fad_liver_counts <- liver_counts[liver_counts$Group2 == "HEMI 5xFAD", ] #31 5xFAD samples
fad_4mo_liver_counts <- fad_liver_counts[fad_liver_counts$Age == "4", ] #18 4mo
fad_4mo_f_liver_counts <- fad_4mo_liver_counts[fad_4mo_liver_counts$Sex == "F", ] #9 females
fad_4mo_m_liver_counts <- fad_4mo_liver_counts[fad_4mo_liver_counts$Sex == "M", ] #9 males
fad_12mo_liver_counts <- fad_liver_counts[fad_liver_counts$Age == "12", ] #13 12mo
fad_12mo_f_liver_counts <- fad_12mo_liver_counts[fad_12mo_liver_counts$Sex == "F", ] #6 females
fad_12mo_m_liver_counts <- fad_12mo_liver_counts[fad_12mo_liver_counts$Sex == "M", ] #7 males
t2_liver_counts <- liver_counts[liver_counts$Group2 == "WT 5xFAD_TREM2", ] #21 TREM2 samples
t2_4mo_liver_counts <- t2_liver_counts[t2_liver_counts$Age == "4", ] #10 4mo
t2_4mo_f_liver_counts <- t2_4mo_liver_counts[t2_4mo_liver_counts$Sex == "F", ] #5 females
t2_4mo_m_liver_counts <- t2_4mo_liver_counts[t2_4mo_liver_counts$Sex == "M", ] #5 males
t2_12mo_liver_counts <- t2_liver_counts[t2_liver_counts$Age == "12", ] #11 12mo
t2_12mo_f_liver_counts <- t2_12mo_liver_counts[t2_12mo_liver_counts$Sex == "F", ] #2 females
t2_12mo_m_liver_counts <- t2_12mo_liver_counts[t2_12mo_liver_counts$Sex == "M", ] #9 males
fadt2_liver_counts <- liver_counts[liver_counts$Group2 == "HEMI 5xFAD_TREM2", ] #19 5xFAD, TREM2 samples
fadt2_4mo_liver_counts <- fadt2_liver_counts[fadt2_liver_counts$Age == "4", ] #10 4mo
fadt2_4mo_f_liver_counts <- fadt2_4mo_liver_counts[fadt2_4mo_liver_counts$Sex == "F", ] #5 females
fadt2_4mo_m_liver_counts <- fadt2_4mo_liver_counts[fadt2_4mo_liver_counts$Sex == "M", ] #5 males
fadt2_12mo_liver_counts <- fadt2_liver_counts[fadt2_liver_counts$Age == "12", ] #9 12mo
fadt2_12mo_f_liver_counts <- fadt2_12mo_liver_counts[fadt2_12mo_liver_counts$Sex == "F", ] #2 females
fadt2_12mo_m_liver_counts <- fadt2_12mo_liver_counts[fadt2_12mo_liver_counts$Sex == "M", ] #7 males
brain_counts <- tar_df_combined[tar_df_combined$Mode == "brain", ] #133 brain samples
wt_brain_counts <- brain_counts[brain_counts$Group2 == "WT 5xFAD", ] #45 WT samples
wt_4mo_brain_counts <- wt_brain_counts[wt_brain_counts$Age == "4", ] #17 4mo
wt_4mo_f_brain_counts <- wt_4mo_brain_counts[wt_4mo_brain_counts$Sex == "F", ] #8 females
wt_4mo_m_brain_counts <- wt_4mo_brain_counts[wt_4mo_brain_counts$Sex == "M", ] #9 males
wt_12mo_brain_counts <- wt_brain_counts[wt_brain_counts$Age == "12", ] #28 12mo
wt_12mo_f_brain_counts <- wt_12mo_brain_counts[wt_12mo_brain_counts$Sex == "F", ] #12 females
wt_12mo_m_brain_counts <- wt_12mo_brain_counts[wt_12mo_brain_counts$Sex == "M", ] #16 males
fad_brain_counts <- brain_counts[brain_counts$Group2 == "HEMI 5xFAD", ] #42 5xFAD samples
fad_4mo_brain_counts <- fad_brain_counts[fad_brain_counts$Age == "4", ] #17 4mo
fad_4mo_f_brain_counts <- fad_4mo_brain_counts[fad_4mo_brain_counts$Sex == "F", ] #9 females
fad_4mo_m_brain_counts <- fad_4mo_brain_counts[fad_4mo_brain_counts$Sex == "M", ] #8 males
fad_12mo_brain_counts <- fad_brain_counts[fad_brain_counts$Age == "12", ] #25 12mo
fad_12mo_f_brain_counts <- fad_12mo_brain_counts[fad_12mo_brain_counts$Sex == "F", ] #12 females
fad_12mo_m_brain_counts <- fad_12mo_brain_counts[fad_12mo_brain_counts$Sex == "M", ] #13 males
t2_brain_counts <- brain_counts[brain_counts$Group2 == "WT 5xFAD_TREM2", ] #24 TREM2 samples
t2_4mo_brain_counts <- t2_brain_counts[t2_brain_counts$Age == "4", ] #10 4mo
t2_4mo_f_brain_counts <- t2_4mo_brain_counts[t2_4mo_brain_counts$Sex == "F", ] #5 females
t2_4mo_m_brain_counts <- t2_4mo_brain_counts[t2_4mo_brain_counts$Sex == "M", ] #5 males
t2_12mo_brain_counts <- t2_brain_counts[t2_brain_counts$Age == "12", ] #14 12mo
t2_12mo_f_brain_counts <- t2_12mo_brain_counts[t2_12mo_brain_counts$Sex == "F", ] #8 females
t2_12mo_m_brain_counts <- t2_12mo_brain_counts[t2_12mo_brain_counts$Sex == "M", ] #6 males
fadt2_brain_counts <- brain_counts[brain_counts$Group2 == "HEMI 5xFAD_TREM2", ] #22 5xFAD, TREM2 samples
fadt2_4mo_brain_counts <- fadt2_brain_counts[fadt2_brain_counts$Age == "4", ] #9 4mo
fadt2_4mo_f_brain_counts <- fadt2_4mo_brain_counts[fadt2_4mo_brain_counts$Sex == "F", ] #4 females
fadt2_4mo_m_brain_counts <- fadt2_4mo_brain_counts[fadt2_4mo_brain_counts$Sex == "M", ] #5 males
fadt2_12mo_brain_counts <- fadt2_brain_counts[fadt2_brain_counts$Age == "12", ] #13 12mo
fadt2_12mo_f_brain_counts <- fadt2_12mo_brain_counts[fadt2_12mo_brain_counts$Sex == "F", ] #8 females
fadt2_12mo_m_brain_counts <- fadt2_12mo_brain_counts[fadt2_12mo_brain_counts$Sex == "M", ] #5 males
##Table S1
library(ggVennDiagram)
set.seed(20231214)
#WT tissues
brain <- wt_brain_counts$Name
plasma <- wt_plasma_counts$Name
liver <- wt_liver_counts$Name
wt <- list(brain=brain,
plasma=plasma,
liver=liver)
ggVennDiagram(wt) +
scale_fill_gradient(low="grey90",high ="red")
ggVennDiagram(wt, show_intersect = TRUE)
#5xFAD tissues
brain <- fad_brain_counts$Name
plasma <- fad_plasma_counts$Name
liver <- fad_liver_counts$Name
fad <- list(brain=brain,
plasma=plasma,
liver=liver)
ggVennDiagram(fad) +
scale_fill_gradient(low="grey90",high ="red")
ggVennDiagram(fad, show_intersect = TRUE)
#TREM2 tissues
brain <- t2_brain_counts$Name
plasma <- t2_plasma_counts$Name
liver <- t2_liver_counts$Name
t2 <- list(brain=brain,
plasma=plasma,
liver=liver)
ggVennDiagram(t2) +
scale_fill_gradient(low="grey90",high ="red")
ggVennDiagram(t2, show_intersect = TRUE)
#5xFAD, TREM2 tissues
brain <- fadt2_brain_counts$Name
plasma <- fadt2_plasma_counts$Name
liver <- fadt2_liver_counts$Name
fadt2 <- list(brain=brain,
plasma=plasma,
liver=liver)
ggVennDiagram(fadt2) +
scale_fill_gradient(low="grey90",high ="red")
ggVennDiagram(fadt2, show_intersect = TRUE)
##Kruskal-Wallis test with post-hoc Dunn’s test (all ages)
#Extract relevant columns
kw_df <- weight_df[, c("Group2", "Age", "Sex", "Weight (g)")]
kw_df <- na.omit(kw_df)
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
shapiro.test(kw_df$`Weight (g)`)
# W = 0.98562, p-value = 0.104
# data IS normally distributed
anova <- aov(`Weight (g)` ~ Group2, data=kw_df)
summary(anova)
# Df Sum Sq Mean Sq F value Pr(>F)
# Group2 3 537 179.00 5.692 0.00101 **
# Residuals 153 4812 31.45
TukeyHSD(anova)
# diff lwr upr p adj
# HEMI 5xFAD_TREM2-HEMI 5xFAD -2.8470588 -6.4035394 0.7094217 0.1645015
# WT 5xFAD-HEMI 5xFAD 2.2177560 -0.6265834 5.0620954 0.1832385
# WT 5xFAD_TREM2-HEMI 5xFAD 2.3455338 -1.1214281 5.8124956 0.2980466
# WT 5xFAD-HEMI 5xFAD_TREM2 5.0648148 1.5409812 8.5886484 0.0015034 **
# WT 5xFAD_TREM2-HEMI 5xFAD_TREM2 5.1925926 1.1494550 9.2357302 0.0058232 **
# WT 5xFAD_TREM2-WT 5xFAD 0.1277778 -3.3056861 3.5612416 0.9996751
##Violin plot
TREM2_tar_df = kw_df
TREM2_tar_df <- TREM2_tar_df %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD_TREM2", "TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD", "WT")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD_TREM2", "5xFAD*TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD", "5xFAD"))
gwsC <- c("forestgreen", "gold", "steelblue1", "firebrick3")
p <- ggplot(TREM2_tar_df, aes(x=Group2, y=`Weight (g)`, fill=Group2)) +
geom_violin(trim=FALSE, scale="area") +
theme_bw() +
#ylim(-3.2,2.2) +
stat_summary(fun=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD"), values = gwsC) +
labs(y = "Weight (g)") +
scale_x_discrete(limits=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p
##Kruskal-Wallis test with post-hoc Dunn’s test (4-months)
#Extract relevant columns
kw_df <- weight_df[, c("Group2", "Age", "Sex", "Weight (g)")]
kw_df <- na.omit(kw_df)
kw_df = kw_df %>% filter(Age == 4)
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
shapiro.test(kw_df$`Weight (g)`)
# W = 0.97723, p-value = 0.1951
# data IS normally distributed
anova <- aov(`Weight (g)` ~ Group2, data=kw_df)
summary(anova)
# Df Sum Sq Mean Sq F value Pr(>F)
# Group2 3 78.5 26.16 0.842 0.475
# Residuals 71 2204.9 31.05
##Violin plot
TREM2_tar_df = kw_df
TREM2_tar_df <- TREM2_tar_df %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD_TREM2", "TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD", "WT")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD_TREM2", "5xFAD*TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD", "5xFAD"))
gwsC <- c("forestgreen", "gold", "steelblue1", "firebrick3")
p1 <- ggplot(TREM2_tar_df, aes(x=Group2, y=`Weight (g)`, fill=Group2)) +
geom_boxplot(trim=FALSE, scale="area") +
theme_bw() +
#ylim(-3.2,2.2) +
stat_summary(fun=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD"), values = gwsC) +
labs(y = "Weight (g)") +
scale_x_discrete(limits=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1
##Kruskal-Wallis test with post-hoc Dunn’s test (12-months)
#Extract relevant columns
kw_df <- weight_df[, c("Group2", "Age", "Sex", "Weight (g)")]
kw_df <- na.omit(kw_df)
kw_df = kw_df %>% filter(Age == 12)
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
shapiro.test(kw_df$`Weight (g)`)
# W = 0.98039, p-value = 0.2429
# data IS normally distributed
anova <- aov(`Weight (g)` ~ Group2, data=kw_df)
summary(anova)
# Df Sum Sq Mean Sq F value Pr(>F)
# Group2 3 648 216.00 7.488 0.000181 ***
# Residuals 78 2250 28.85
TukeyHSD(anova)
# diff lwr upr p adj
# HEMI 5xFAD_TREM2-HEMI 5xFAD -2.6681159 -7.3477490 2.011517 0.4443203
# WT 5xFAD-HEMI 5xFAD 4.4429952 0.4419933 8.443997 0.0235102 *
# WT 5xFAD_TREM2-HEMI 5xFAD 4.0416880 -0.4682562 8.551632 0.0950913
# WT 5xFAD-HEMI 5xFAD_TREM2 7.1111111 2.5703706 11.651852 0.0005502 ***
# WT 5xFAD_TREM2-HEMI 5xFAD_TREM2 6.7098039 1.7148159 11.704792 0.0038780 **
# WT 5xFAD_TREM2-WT 5xFAD -0.4013072 -4.7669635 3.964349 0.9950176
##Violin plot
TREM2_tar_df = kw_df
TREM2_tar_df <- TREM2_tar_df %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD_TREM2", "TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "WT 5xFAD", "WT")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD_TREM2", "5xFAD*TREM2")) %>%
mutate(., Group2 = stringr::str_replace(Group2, "HEMI 5xFAD", "5xFAD"))
gwsC <- c("forestgreen", "gold", "steelblue1", "firebrick3")
p2 <- ggplot(TREM2_tar_df, aes(x=Group2, y=`Weight (g)`, fill=Group2)) +
geom_boxplot(trim=FALSE, scale="area") +
theme_bw() +
#ylim(-3.2,2.2) +
stat_summary(fun=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD"), values = gwsC) +
labs(y = "Weight (g)") +
scale_x_discrete(limits=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2
ggarrange(p1, p2, ncol=2, nrow=1, common.legend = TRUE, legend="right")
ggsave("Figure_for_review.png", width = 11, height = 4.5)
##############
##############PERMANOVAs##############
#PERMANOVA (Supplemental Figures 1, 2, and 3)
#Repeat for all Mode-Age combinations (liver, plasma, brain - 4, 12)
#To calculate for males and females separately, add Sex == "M" or "F" to lines 155 and 160, and use ~ Group2 in lines 157 and 162
lip_perma_df = lip_df_combined %>%
filter(Mode == "brain", Age == 12) %>%
column_to_rownames(var = "sample_id")
pairwise.adonis2(select(lip_perma_df, !Name:Background) ~ Group2*Sex, data = lip_perma_df, nperm = 999, method = "euclidean")
tar_perma_df = tar_df_combined %>%
filter(Mode == "brain", Age == 12) %>%
column_to_rownames(var = "sample_id")
pairwise.adonis2(select(tar_perma_df, !Name:Background) ~ Group2*Sex, data = tar_perma_df, nperm = 999, method = "euclidean")
#MAIN EFFECTS PERMANOVA (Figure 1C)
#Repeat for all Modes (liver, plasma, brain)
lip_perma_df = na.omit(lip_df_combined) %>%
filter(Mode == "brain") %>% # Comment this line out when comparing all tissues, use ~ Group2+Sex+Age+Mode in line 168 below
column_to_rownames(var = "sample_id")
adonis2(select(lip_perma_df, !Name:Background) ~ Group2+Sex+Age, data = lip_perma_df, permutations = 9999, method = "euclidean")
tar_perma_df = na.omit(tar_df_combined) %>%
filter(Mode == "brain") %>% # Comment this line out when comparing all tissues, use ~ Group2+Sex+Age+Mode in line 173 below
column_to_rownames(var = "sample_id")
adonis2(select(tar_perma_df, !Name:Background) ~ Group2+Sex+Age, data = tar_perma_df, permutations = 9999, method = "euclidean")
##############
##############PCAs##############
#Metadata
tar_metadata <- tar_df_combined[,c(1:12)]
lip_metadata <- lip_df_combined[,c(1:12)]
tar_data_pcoa <- tar_df_combined %>% remove_rownames %>% column_to_rownames(var="sample_id")
lip_data_pcoa <- lip_df_combined %>% remove_rownames %>% column_to_rownames(var="sample_id")
##Targmet
dist_matrix = tar_data_pcoa %>%
select(!Name:Background) %>%
vegdist(., method = "euclidean")
#PCoA - euclidean, with eigenvalues, calculated and normalized variance
pcoa = cmdscale(dist_matrix, eig = T, k = nrow(tar_data_pcoa)-1, add = T)
pcoa_eig = eigenvals(pcoa)
pcoa_var = pcoa_eig/sum(pcoa_eig)
pcoa_var[1:3] #First 3 axes
pcoa_plot_df = as.data.frame(pcoa$points[,1:2]) %>%
rownames_to_column(var = "sample_id") %>%
inner_join(tar_metadata)
shape_values<-c(21,24,22,3)
gwsC <- c("firebrick3", "steelblue1", "gold", "forestgreen")
ggplot(data = pcoa_plot_df, aes(x = V1, y = V2, fill = Group2)) +
#stat_ellipse(aes(group = Mode, linetype = Mode, color = Group2), level = 0.95, show.legend = F, linewidth = 0.6) +
geom_point(size = 3, alpha = 0.75, aes(shape=Mode)) +
#geom_text(aes(label = sample_id)) +
theme_bw() +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
labs(fill = "Genotype") +
labs(shape = "Tissue") +
scale_shape_manual(values=shape_values) +
scale_color_manual(values=gwsC) + scale_fill_manual(values=gwsC)+
labs(x = bquote("PC1:"~.(round(pcoa_var[1]*100, digits = 1))~"%"),
y = bquote("PC2:"~.(round(pcoa_var[2]*100, digits = 1))~"%"),
title = "Targeted Metabolomics",
subtitle = "Log10 normalized")
ggsave("PCA-targmet.png", width = 6, height = 4.25)
#PERMDISP (repeat for brain, plasma, and liver)
permdisp = tar_data_pcoa %>%
filter(Mode == "liver")
dist_matrix = permdisp %>%
select(!Name:Background) %>%
vegdist(., method = "euclidean")
pcoa.permdisp <- betadisper(dist_matrix,
permdisp$Group2,
type = c("median", "centroid"),
bias.adjust = FALSE,
sqrt.dist = FALSE,
add = FALSE)
anova(pcoa.permdisp)
# Brain:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 0.963 0.32113 0.7277 0.5372
# Residuals 129 56.925 0.44128
#
# Plasma:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 7.96 2.6521 0.9803 0.4036
# Residuals 160 432.86 2.7054
#
# Liver:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 23.067 7.6889 3.1431 0.02852 *
# Residuals 101 247.075 2.4463
#
# TukeyHSD(pcoa.permdisp)
# diff lwr upr p adj
# HEMI 5xFAD_TREM2-HEMI 5xFAD -0.4917491 -1.682191 0.6986930 0.7030365
# WT 5xFAD-HEMI 5xFAD -0.7358761 -1.750528 0.2787763 0.2368776
# WT 5xFAD_TREM2-HEMI 5xFAD -1.3257969 -2.480558 -0.1710360 0.0176508 *
# WT 5xFAD-HEMI 5xFAD_TREM2 -0.2441270 -1.414441 0.9261875 0.9477112
# WT 5xFAD_TREM2-HEMI 5xFAD_TREM2 -0.8340478 -2.127721 0.4596250 0.3373487
# WT 5xFAD_TREM2-WT 5xFAD -0.5899208 -1.723921 0.5440794 0.5279217
#
# (Only the TREM2 vs 5xFAD genotype comparison was significant after post-hoc testing)
##Lipid
dist_matrix = lip_data_pcoa %>%
select(!Name:Background) %>%
vegdist(., method = "euclidean")
#PCoA - euclidean, with eigenvalues, calculated and normalized variance
pcoa = cmdscale(dist_matrix, eig = T, k = nrow(lip_data_pcoa)-1, add = T)
pcoa_eig = eigenvals(pcoa)
pcoa_var = pcoa_eig/sum(pcoa_eig)
pcoa_var[1:3] #First 3 axes
pcoa_plot_df = as.data.frame(pcoa$points[,1:2]) %>%
rownames_to_column(var = "sample_id") %>%
inner_join(tar_metadata)
shape_values<-c(21,24,22,3)
gwsC <- c("firebrick3", "steelblue1", "gold", "forestgreen")
ggplot(data = pcoa_plot_df, aes(x = V1, y = V2, fill = Group2)) +
#stat_ellipse(aes(group = Mode, linetype = Mode, color = Group2), level = 0.95, show.legend = F, linewidth = 0.6) +
geom_point(size = 3, alpha = 0.75, aes(shape=Mode)) +
#geom_text(aes(label = sample_id)) +
theme_bw() +
guides(fill = guide_legend(override.aes = list(shape = 21))) +
labs(fill = "Genotype") +
labs(shape = "Tissue") +
scale_shape_manual(values=shape_values) +
scale_color_manual(values=gwsC) + scale_fill_manual(values=gwsC)+
labs(x = bquote("PC1:"~.(round(pcoa_var[1]*100, digits = 1))~"%"),
y = bquote("PC2:"~.(round(pcoa_var[2]*100, digits = 1))~"%"),
title = "Lipidomics",
subtitle = "Log10 normalized")
ggsave("PCA-lip.png", width = 6, height = 4.25)
#PERMDISP (repeat for brain, plasma, and liver)
permdisp = lip_data_pcoa %>%
filter(Mode == "brain")
dist_matrix = permdisp %>%
select(!Name:Background) %>%
vegdist(., method = "euclidean")
pcoa.permdisp <- betadisper(dist_matrix,
permdisp$Group2,
type = c("median", "centroid"),
bias.adjust = FALSE,
sqrt.dist = FALSE,
add = FALSE)
anova(pcoa.permdisp)
# Brain:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 24.14 8.0474 1.1985 0.3131
# Residuals 128 859.50 6.7148
#
# Plasma:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 42.16 14.054 1.2023 0.3108
# Residuals 159 1858.60 11.689
#
# Liver:
# Response: Distances
# Df Sum Sq Mean Sq F value Pr(>F)
# Groups 3 25.69 8.5625 1.2408 0.299
# Residuals 101 696.99 6.9009
##############
##############Volcano plots and DA metabolite lists#############
##Figure2
#Liver Trem2 v WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "WT 5xFAD_TREM2"), Age == 4, Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`WT 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_volc
LTM4_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
LTM4_DF <-select(LTM4_DF, var, Log2FC, padj)
write.csv(LTM4_DF, "Figure2D.csv", row.names=FALSE)
ggsave("Figure2C.png", LTM4_volc, width = 7, height = 4)
#Plasma Trem2 vs WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "WT 5xFAD_TREM2"), Age == 4, Mode == "plasma") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`WT 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
PTM4_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
PTM4_volc
PTM4_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
PTM4_DF <-select(PTM4_DF, var, Log2FC, padj)
write.csv(PTM4_DF, "Figure2F.csv", row.names=FALSE)
ggsave("Figure2E.png", PTM4_volc, width = 7, height = 4)
#Brain Trem2 vs WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "WT 5xFAD_TREM2"), Age == 4, Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`WT 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_volc
BTM4_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
BTM4_DF <-select(BTM4_DF, var, Log2FC, padj)
write.csv(BTM4_DF, "Figure2H.csv", row.names=FALSE)
ggsave("Figure2G.png", BTM4_volc, width = 7, height = 4)
##SuppFigure5
#Liver targ.met 5xFAD*Trem2 vs Trem2
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2", "HEMI 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_12_T2_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_12_T2_5xT2_volc
# LTM4_12_T2_5xT2_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
# LTM4_12_T2_5xT2_DF <-select(LTM4_12_T2_5xT2_DF, var, Log2FC, padj)
# write.csv(LTM4_12_T2_5xT2_DF, "LTM4_12_T2_5xT2_DF_Fig3.csv", row.names=FALSE)
ggsave("SuppFig5A.png", LTM4_12_T2_5xT2_volc, width = 7, height = 4)
#Plasma targ.met Trem2 vs 5xFAD*Trem2 (no signifigant metabolites)
#Brain targ.met Trem2 vs 5xFAD*Trem2 (no signifigant metabolites)
#*Liver lipid 5xFAD*Trem2 vs Trem2
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2", "HEMI 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LL4_12_T2_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LL4_12_T2_5xT2_volc
# LL4_12_T2_5xT2_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
# LL4_12_T2_5xT2_DF <-select(LL4_12_T2_5xT2_DF, var, Log2FC, padj)
# write.csv(LL4_12_T2_5xT2_DF, "LL4_12_T2_5xT2_DF_Fig3.csv", row.names=FALSE)
ggsave("SuppFig5B.png", LL4_12_T2_5xT2_volc, width = 7, height = 4)
#Plasma lipid Trem2 vs 5xFAD*Trem2 (no signifigant metabolites)
#*Brain lipid Trem2 vs 5xFAD*Trem2 (no signifigant metabolites)
##Figure2
#Liver targ.met WT vs TREM2
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "WT 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`WT 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_12_WT_T2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = TRUE) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_12_WT_T2_volc
# LTM4_12_WT_T2_DF <- log2fc %>% filter(padj < 0.05)
# LTM4_12_WT_T2_DF <-select(LTM4_12_WT_T2_DF, var, Log2FC, padj)
# write.csv(LTM4_12_WT_T2_DF, "LTM4_12_WT_T2_DF_Fig3.csv", row.names=FALSE)
ggsave("Figure2A.png", LTM4_12_WT_T2_volc, width = 7, height = 4)
#Plasma targ.met WT vs TREM2 (no signifigant metabolites)
#Brain targ.met WT vs TREM2
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "WT 5xFAD_TREM2"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`WT 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_12_WT_T2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = TRUE) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_12_WT_T2_volc
# BTM4_12_WT_T2_DF <- log2fc %>% filter(padj < 0.05)
# BTM4_12_WT_T2_DF <-select(BTM4_12_WT_T2_DF, var, Log2FC, padj)
# write.csv(BTM4_12_WT_T2_DF, "BTM4_12_WT_T2_DF_Fig3.csv", row.names=FALSE)
ggsave("Figure2B.png", BTM4_12_WT_T2_volc, width = 7, height = 4)
#Liver lipid WT vs TREM2 (no signifigant metabolites)
#Plasma lipid WT vs TREM2 (no signifigant metabolites)
#Brain lipid WT vs TREM2 (no signifigant metabolites)
#Liver targ.met WT vs 5xFAD (no signifigant metabolites)
#Plasma targ.met WT vs 5xFAD (no signifigant metabolites)
#Brain targ.met WT vs 5xFAD (no signifigant metabolites)
#Liver lipid WT vs 5xFAD (no signifigant metabolites)
#Brain lipid WT vs 5xFAD (no signifigant metabolites)
#Plasma lipid WT vs 5xFAD (no signifigant metabolites)
##Figure4
#Liver targ.met 5xFAD*Trem2 vs 5xFAD
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_12_5x_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(title = "", subtitle = "", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_12_5x_5xT2_volc
ggsave("Figure4A.png", LTM4_12_5x_5xT2_volc, width = 4.5, height = 4)
#Plasma targ.met 5xFAD*Trem2 vs 5xFAD
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "plasma") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
PTM4_12_5x_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(title = "", subtitle = "", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
PTM4_12_5x_5xT2_volc
ggsave("Figure4B.png", PTM4_12_5x_5xT2_volc, width = 4.5, height = 4)
#Brain targ.met 5xFAD*Trem2 vs 5xFAD
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_12_5x_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(title = "", subtitle = "", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_12_5x_5xT2_volc
ggsave("Figure4C.png", BTM4_12_5x_5xT2_volc, width = 4.5, height = 4)
#Liver lipid 5xFAD vs 5xFAD_TREM2 (no signifigant metabolites)
#Plasma lipid 5xFAD vs 5xFAD_TREM2 (no signifigant metabolites)
#Brain lipid 5xFAD vs 5xFAD_TREM2 (no signifigant metabolites)
#Liver targ.met Trem2 vs 5xFAD (no signifigant metabolites)
##SuppFigure4
#Plasma targ.met 5xFAD vs Trem2
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2", "HEMI 5xFAD"), Mode == "plasma") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD` / log2fc$`WT 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
PTM4_12_T2_5x_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(title = "5xFAD vs. TREM2", subtitle = "Plasma polar metabolites", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
PTM4_12_T2_5x_volc
ggsave("SuppFig4A.png", PTM4_12_T2_5x_volc, width = 4.5, height = 4)
#Brain targ.met 5xFAD vs Trem2
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2", "HEMI 5xFAD"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD` / log2fc$`WT 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_12_T2_5x_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(title = "", subtitle = "Brain polar metabolites", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_12_T2_5x_volc
ggsave("SuppFig4B.png", BTM4_12_T2_5x_volc, width = 4.5, height = 4)
#Liver lipid 5xFAD vs Trem2
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2", "HEMI 5xFAD"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD` / log2fc$`WT 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LL4_12_T2_5x_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue", "#BF0D3E"), labels = c("N.S.", "Decreased", "Increased")) +
labs(title = "", subtitle = "Liver lipids", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LL4_12_T2_5x_volc
ggsave("SuppFig4C.png", LL4_12_T2_5x_volc, width = 4.5, height = 4)
#Plasma lipid Trem2 vs 5xFAD (no signifigant metabolites)
#Brain lipid Trem2 vs 5xFAD (no signifigant metabolites)
#Liver targ.met 5xFAD*TREM2 vs WT
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_12_WT_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(title = "5xFAD*TREM2 vs. WT", subtitle = "Liver polar metabolites", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_12_WT_5xT2_volc
ggsave("SuppFig4D.png", LTM4_12_WT_5xT2_volc, width = 4.5, height = 4)
#Plasma targ.met WT vs 5xFAD*TREM2 (no signifigant metabolites)
#Brain targ.met 5xFAD*TREM2 vs WT
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_12_WT_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "#BF0D3E"), labels = c("N.S.", "Increased")) +
labs(title = "", subtitle = "Brain polar metabolites", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_12_WT_5xT2_volc
ggsave("SuppFig4E.png", BTM4_12_WT_5xT2_volc, width = 4.5, height = 4)
#Liver lipid 5xFAD*Trem2 vs WT
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LL4_12_WT_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "steelblue", "#BF0D3E"), labels = c("N.S.", "Decreased", "Increased")) +
labs(title = "", subtitle = "Liver lipids", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LL4_12_WT_5xT2_volc
ggsave("SuppFig4F.png", LL4_12_WT_5xT2_volc, width = 4.5, height = 4)
#Plasma lipid WT vs 5xFAD*Trem2 (no signifigant metabolites)
#Brain lipid 5xFAD*Trem2 vs WT
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BL4_12_WT_5xT2_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2, show.legend = FALSE) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(title = "", subtitle = "Brain lipids", x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.2, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BL4_12_WT_5xT2_volc
ggsave("SuppFig4G.png", BL4_12_WT_5xT2_volc, width = 4.5, height = 4)
##Figure3
#Liver 5xFAD vs WT 12mo. lipid
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD", "HEMI 5xFAD"), Age == 12, Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD` / log2fc$`WT 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LL12_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LL12_volc
LL12_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
LL12_DF <-select(LL12_DF, var, Log2FC, padj)
write.csv(LL12_DF, "Figure3B.csv", row.names=FALSE)
ggsave("Figure3A.png", LL12_volc, width = 7, height = 4)
##Figure5
#Brain 5xFAD*TREM2 vs. 5xFAD 4mo. lipid
TREM2_lip_df = lip_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Age == 4, Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BL12_volc <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BL12_volc
BL12_DF <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
BL12_DF <-select(BL12_DF, var, Log2FC, padj)
write.csv(BL12_DF, "Figure5H.csv", row.names=FALSE)
ggsave("Figure5G.png", BL12_volc, width = 7, height = 4)
#Liver 5xFAD*Trem2 vs 5xFAD 4mo. targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Age == 4, Mode == "liver") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
LTM4_volc_2 <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "steelblue"), labels = c("N.S.", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
LTM4_volc_2
LTM4_DF_2 <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
LTM4_DF_2 <-select(LTM4_DF_2, var, Log2FC, padj)
write.csv(LTM4_DF_2, "Figure5B.csv", row.names=FALSE)
ggsave("Figure5A.png", LTM4_volc_2, width = 7, height = 4)
#Plasma 5xFAD*Trem2 vs 5xFAD 4mo. targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Age == 4, Mode == "plasma") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
PTM4_volc_2 <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
PTM4_volc_2
PTM4_DF_2 <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
PTM4_DF_2 <-select(PTM4_DF_2, var, Log2FC, padj)
write.csv(PTM4_DF_2, "Figure5D.csv", row.names=FALSE)
ggsave("Figure5C.png", PTM4_volc_2, width = 7, height = 4)
#Brain 5xFAD*Trem2 vs. 5xFAD 4mo. targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Age == 4, Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>%
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD_TREM2` / log2fc$`HEMI 5xFAD`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#Volcano plot
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
BTM4_volc_2 <- ggplot(data = log2fc) +
aes(x = Log2FC, y = -log10(padj), color = color) +
geom_hline(yintercept = -log10(0.05), lty = 2) +
geom_vline(xintercept = 0.58496250072, lty = 3) +
geom_vline(xintercept = -0.58496250072, lty = 3) +
geom_point(alpha = 0.5, size = 2) +
scale_color_manual(values = c("gray", "#BF0D3E", "steelblue"), labels = c("N.S.", "Increased", "Decreased")) +
labs(x = expression("log"[2]*"(FC)"), y = expression("-log"[10]*"(p-adj)"), color = "Abundance") +
theme_bw() +
#theme(plot.title = element_text(size = 18)) +
ggrepel::geom_text_repel(aes(x = Log2FC, y = -log10(padj), label = var), color = "black", size = 2.6, max.overlaps = 8) #+
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 2, vjust = -8.25, label = bquote("Total:"~.(sum(TREM2_tar_res_df$qval < 0.05)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.525, vjust = -6.75, label = bquote("Increased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef > 0)))) +
#annotate("text", x = Inf, y = -Inf, size = 3, hjust = 1.55, vjust = -5, label = bquote("Decreased:"~.(sum(TREM2_tar_res_df$qval < 0.05 & TREM2_tar_res_df$coef < 0))))
BTM4_volc_2
BTM4_DF_2 <- log2fc %>% filter((padj < 0.05 & Log2FC > 0.58496250072) | (padj < 0.05 & Log2FC < -0.58496250072))
BTM4_DF_2 <-select(BTM4_DF_2, var, Log2FC, padj)
write.csv(BTM4_DF_2, "Figure5F.csv", row.names=FALSE)
ggsave("Figure5E.png", BTM4_volc_2, width = 7, height = 4)
##############
##############5-MTHF abundance p-values and Violin plots##############
##Kruskal-Wallis test with post-hoc Dunn’s test
#Extract 4-month brain samples
TREM2_tar_df = tar_df_combined %>%
filter(Mode == "brain", Age == "4") %>%
column_to_rownames(var = "sample_id")
#Extract relevant columns
kw_df <- TREM2_tar_df[, c("Group2", "5-METHYLTETRAHYDROFOLIC ACID_pos_4")]
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
shapiro.test(kw_df$`5-METHYLTETRAHYDROFOLIC ACID_pos_4`)
# W = 0.49383, p-value = 2.953e-12
# data is NOT normally distributed
kruskal.test(`5-METHYLTETRAHYDROFOLIC ACID_pos_4` ~ Group2, data=kw_df)
# Kruskal-Wallis chi-squared = 32.924, df = 3,
# p-value = 3.342e-07
dunnTest(`5-METHYLTETRAHYDROFOLIC ACID_pos_4` ~ Group2, data=kw_df, method="bonferroni")
# Comparison Z P.unadj P.adj
# 1 HEMI 5xFAD - HEMI 5xFAD_TREM2 -3.9341319 8.349791e-05 5.009875e-04 ***
# 2 HEMI 5xFAD - WT 5xFAD 0.4664076 6.409238e-01 1.000000e+00 n.s.
# 3 HEMI 5xFAD_TREM2 - WT 5xFAD 4.3222065 1.544765e-05 9.268589e-05 ***
# 4 HEMI 5xFAD - WT 5xFAD_TREM2 -3.7570979 1.718953e-04 1.031372e-03 **
# 5 HEMI 5xFAD_TREM2 - WT 5xFAD_TREM2 0.2708968 7.864704e-01 1.000000e+00 n.s.
# 6 WT 5xFAD - WT 5xFAD_TREM2 -4.1585176 3.203195e-05 1.921917e-04 ***
#Extract 12-month brain samples
TREM2_tar_df = tar_df_combined %>%
filter(Mode == "brain", Age == "12") %>%
column_to_rownames(var = "sample_id")
#Extract relevant columns
kw_df <- TREM2_tar_df[, c("Group2", "5-METHYLTETRAHYDROFOLIC ACID_pos_4")]
#Shapiro-Wilk normality test
#if p < 0.05 the data is not normally distributed
#if p > 0.05 the data is normally distributed
shapiro.test(kw_df$`5-METHYLTETRAHYDROFOLIC ACID_pos_4`)
# W = 0.47607, p-value = 2.088e-15
# data is NOT normally distributed
kruskal.test(`5-METHYLTETRAHYDROFOLIC ACID_pos_4` ~ Group2, data=kw_df)
# Kruskal-Wallis chi-squared = 7.0677, df = 3,
# p-value = 0.06977
dunnTest(`5-METHYLTETRAHYDROFOLIC ACID_pos_4` ~ Group2, data=kw_df, method="bonferroni")
# Comparison Z P.unadj P.adj
# 1 HEMI 5xFAD - HEMI 5xFAD_TREM2 -1.9675186 0.04912346 0.29474076 n.s.
# 2 HEMI 5xFAD - WT 5xFAD 0.4890598 0.62479939 1.00000000 n.s.
# 3 HEMI 5xFAD_TREM2 - WT 5xFAD 2.4055679 0.01614735 0.09688412 n.s.
# 4 HEMI 5xFAD - WT 5xFAD_TREM2 -1.1300357 0.25846120 1.00000000 n.s.
# 5 HEMI 5xFAD_TREM2 - WT 5xFAD_TREM2 0.7673536 0.44287128 1.00000000 n.s.
# 6 WT 5xFAD - WT 5xFAD_TREM2 -1.5635354 0.11792672 0.70756034 n.s.
##Violin plots
#4-months
TREM2_tar_df = tar_df_combined %>%
filter(Mode == "brain", Age == "4") %>%
column_to_rownames(var = "sample_id")
#Adjust names and values
TREM2_tar_df[,-(1:12)] = log10(TREM2_tar_df[,-(1:12)]) %>% #Normalization
replace(is.na(.), 0)
TREM2_tar_df = subset(TREM2_tar_df, select = -c(Genotype))
names(TREM2_tar_df) [8] <- c("Genotype")
TREM2_tar_df <- TREM2_tar_df %>%
mutate(., Genotype = stringr::str_replace(Genotype, "WT 5xFAD_TREM2", "TREM2")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "WT 5xFAD", "WT")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "HEMI 5xFAD_TREM2", "5xFAD*TREM2")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "HEMI 5xFAD", "5xFAD"))
gwsC <- c("forestgreen", "gold", "steelblue1", "firebrick3")
p1 <- ggplot(TREM2_tar_df, aes(x=Genotype, y=`5-METHYLTETRAHYDROFOLIC ACID_pos_4`, fill=Genotype)) +
geom_violin(trim=FALSE, scale="area") +
theme_bw() +
ylim(-3.2,2.2) +
stat_summary(fun=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD"), values = gwsC) +
labs(y = expression("Log"[10]*" 5-MTHF Abundance")) +
scale_x_discrete(limits=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p1
ggsave("Figure8a.png", p, width = 8.5, height = 4.5)
#TREM2 and WT only (Figure2H)
TREM2_tar_df2 <- TREM2_tar_df %>%
filter(Genotype %in% c("TREM2", "WT"))
q <- ggplot(TREM2_tar_df2, aes(x=Genotype, y=`5-METHYLTETRAHYDROFOLIC ACID_pos_4`, fill=Genotype)) +
geom_violin(trim=FALSE, scale="area") +
theme_bw() +
stat_summary(fun.y=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT"), values = gwsC) +
labs(y = expression("Log"[10]*" 5-MTHF Abundance")) +
scale_x_discrete(limits=c("TREM2", "WT")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
q
ggsave("Figure2H.png", q, width = 5, height = 4)
#12-months
TREM2_tar_df = tar_df_combined %>%
filter(Mode == "brain", Age == "12") %>%
column_to_rownames(var = "sample_id")
#Adjust names and values
TREM2_tar_df[,-(1:12)] = log10(TREM2_tar_df[,-(1:12)]) %>% #Normalization
replace(is.na(.), 0)
TREM2_tar_df = subset(TREM2_tar_df, select = -c(Genotype))
names(TREM2_tar_df) [8] <- c("Genotype")
TREM2_tar_df <- TREM2_tar_df %>%
mutate(., Genotype = stringr::str_replace(Genotype, "WT 5xFAD_TREM2", "TREM2")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "WT 5xFAD", "WT")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "HEMI 5xFAD_TREM2", "5xFAD*TREM2")) %>%
mutate(., Genotype = stringr::str_replace(Genotype, "HEMI 5xFAD", "5xFAD"))
gwsC <- c("forestgreen", "gold", "steelblue1", "firebrick3")
p2 <- ggplot(TREM2_tar_df, aes(x=Genotype, y=`5-METHYLTETRAHYDROFOLIC ACID_pos_4`, fill=Genotype)) +
geom_violin(trim=FALSE, scale="area") +
theme_bw() +
ylim(-3.2,2.2) +
stat_summary(fun=mean, geom="point", shape=23, size=2, fill="yellow") +
geom_point(aes(shape=Sex), position=position_jitter(0.2)) +
scale_fill_manual(breaks=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD"), values = gwsC) +
labs(y = expression("Log"[10]*" 5-MTHF Abundance")) +
scale_x_discrete(limits=c("TREM2", "WT", "5xFAD*TREM2", "5xFAD")) +
theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
axis.title.x=element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank())
p2
ggsave("Figure8b.png", p, width = 8.5, height = 4.5)
ggarrange(p1, p2, ncol=2, nrow=1, common.legend = TRUE, legend="right")
ggsave("Figure8.png", width = 11, height = 4.5)
##############
##############Plasma lipid DA analysis by Sex##############
#Lipid plasma
TREM2_lip_df = lip_df_combined %>%
filter(Mode == "plasma") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_lip_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>%
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_lip_df$Sex)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_lip_df %>%
group_by(Sex) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>%
column_to_rownames(var = "Sex") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`M`/ log2fc$`F`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
log2fc = log2fc %>%
mutate(color = case_when(
padj < 0.05 & Log2FC > 0.58496250072 ~ "color2",
padj < 0.05 & Log2FC < -0.58496250072 ~ "color3",
padj > 0.05 ~ "color1",
Log2FC > -1 ~ "color1",
Log2FC < 1 ~ "color1"
))
# only significant (X significant DA metabolites)
SEX1 <- log2fc %>% filter(color %in% c("color2", "color3"))
write.csv(SEX1, "plasma-sex.csv", row.names=FALSE)
##############
##############Heatmaps##############
#Figure 6
###Top100 M vs. F
#Repeat for all Modes (liver-plasma-brain) for both lipids and polar metabolites
TREM2_tar_df = tar_df_combined %>% #use tar_df_combined for targ.met, and lip_df_combined for lipids
filter(Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`) %>% #use `1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1` for targ.met and `CE(12:0)`:`TAG56:1-FA18:1` for lipids
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Sex)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Sex) %>%
summarise(across(`1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1`, mean)) %>% #use 1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1` for targ.met and `CE(12:0)`:`TAG56:1-FA18:1` for lipids
column_to_rownames(var = "Sex") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`M`/ log2fc$`F`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#top 100 most differential based on Log2FC
top_100_tar = log2fc %>%
slice_max(order_by = abs(Log2FC), n = 100)
#Plot relative abundances (log10 transformed)
test = tar_df_combined %>% #use tar_df_combined for targ.met, and lip_df_combined for lipids
column_to_rownames(var = "sample_id") %>%
filter(Mode == "brain") %>%
select(top_100_tar$var)
tar_zscore_df = test
tar_zscore_df = log10(tar_zscore_df) %>% #Normalization
replace(is.na(.), 0)
rownames(tar_zscore_df) = rownames(test)
tar_zscore_df = as.data.frame(tar_zscore_df)
tar_zscore_df = replace(tar_zscore_df, is.na(tar_zscore_df), 0)
#Make heatmap
ann_colors = list(Sex = c(`F`="lightpink", `M`="lightblue"),
Genotype = c(`HEMI 5xFAD`="firebrick3", `HEMI 5xFAD_TREM2`="steelblue1", `WT 5xFAD`="gold", `WT 5xFAD_TREM2`="forestgreen"),
Age = c(`12`="azure3", `4`="steelblue4")
)
tar_merged_zscores = tar_zscore_df %>%
select(which(!colSums(.) == 0)) %>%
rownames_to_column(var = "sample_id") %>%
separate_wider_delim(cols = sample_id, names = c("Name", "Mode"), delim = "_", cols_remove = F) %>%
inner_join(metadata) %>%
relocate(CageID:Background, .after = sample_id) %>%
column_to_rownames(var = "sample_id") %>%
select(-Genotype) %>%
rename(Genotype=Group2) %>%
arrange(Mode, Sex, Genotype, Age)
hm2 <- pheatmap(t(select(tar_merged_zscores, !Name:Background)), border_color = "NA",
annotation_col = data.frame(rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Age) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Genotype) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Sex) %>% column_to_rownames(var = "sample_id")
),
annotation_colors = ann_colors,
show_colnames = F, fontsize = 12, fontsize_row = 5.5, cluster_rows = T, cluster_cols = F, cutree_rows = 5, main = "Brain - metabolites")
save_pheatmap(hm2, "HM-Figure6F.png", width=7.5, height=8.5)
#for plasma: use breaks to make tables more different
my.breaks <- c(seq(-1.5, 1.5, by=0.1))
hm3 <- pheatmap(t(select(tar_merged_zscores, !Name:Background)),
border_color = "NA",
color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(my.breaks)),
breaks = my.breaks,
annotation_col = data.frame(rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Age) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Genotype) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Sex) %>% column_to_rownames(var = "sample_id")
),
annotation_colors = ann_colors,
show_colnames = F, fontsize = 12, fontsize_row = 5.5, cluster_rows = T, cluster_cols = F, cutree_rows = 2, main = "Plasma - metabolites") #fontsize row 6 for lipids and 5.5 for targmet
save_pheatmap(hm3, "HM-Figure6D.png", width=7.5, height=8.5)
#Figure S6
###Top100 5xFAD vs. 5xFAD*TREM2
#Repeat for all Modes (liver-plasma-brain) for both lipids and polar metabolites
TREM2_tar_df = lip_df_combined %>% #use tar_df_combined for targ.met, and lip_df_combined for lipids
filter(Group2 %in% c("HEMI 5xFAD", "HEMI 5xFAD_TREM2"), Mode == "brain") %>%
column_to_rownames(var = "sample_id")
#Note: its named t test but its using a wilcox test
t_test_df = TREM2_tar_df %>%
select(`CE(12:0)`:`TAG56:1-FA18:1`) %>% #use `1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1` for targ.met and `CE(12:0)`:`TAG56:1-FA18:1` for lipids
map_df(~ broom::tidy(wilcox.test(. ~ TREM2_tar_df$Group2)), .id = 'var') %>%
select(var, p.value)
t_test_df$padj = p.adjust(t_test_df$p.value, method = "fdr") #changed to FDR, less stringent
#Need to calc fold change
log2fc = TREM2_tar_df %>%
group_by(Group2) %>%
summarise(across(`CE(12:0)`:`TAG56:1-FA18:1`, mean)) %>% #use 1-METHYLADENOSINE_pos_1`:`RIBOFLAVIN_neg_1` for targ.met and `CE(12:0)`:`TAG56:1-FA18:1` for lipids
column_to_rownames(var = "Group2") %>%
t(.)
log2fc = as.data.frame(log2fc)
log2fc$Ratio = log2fc$`HEMI 5xFAD`/ log2fc$`HEMI 5xFAD_TREM2`
log2fc$Log2FC = log(log2fc$Ratio, base = 2)
log2fc[is.na(log2fc)] <- 0 #change any NaN to 0
log2fc = log2fc %>%
rownames_to_column(var = "var") %>%
left_join(t_test_df)
#top 100 most differential based on Log2FC
top_100_tar = log2fc %>%
slice_max(order_by = abs(Log2FC), n = 100)
#Plot relative abundances (log10 transformed)
test = lip_df_combined %>% #use tar_df_combined for targ.met, and lip_df_combined for lipids
column_to_rownames(var = "sample_id") %>%
filter(Mode == "brain") %>%
select(top_100_tar$var)
tar_zscore_df = test
tar_zscore_df = log10(tar_zscore_df) %>% #Normalization
replace(is.na(.), 0)
rownames(tar_zscore_df) = rownames(test)
tar_zscore_df = as.data.frame(tar_zscore_df)
tar_zscore_df = replace(tar_zscore_df, is.na(tar_zscore_df), 0)
#Make heatmap
ann_colors = list(Sex = c(`F`="lightpink", `M`="lightblue"),
Genotype = c(`HEMI 5xFAD`="firebrick3", `HEMI 5xFAD_TREM2`="steelblue1", `WT 5xFAD`="gold", `WT 5xFAD_TREM2`="forestgreen"),
Age = c(`12`="azure3", `4`="steelblue4")
)
tar_merged_zscores = tar_zscore_df %>%
select(which(!colSums(.) == 0)) %>%
rownames_to_column(var = "sample_id") %>%
separate_wider_delim(cols = sample_id, names = c("Name", "Mode"), delim = "_", cols_remove = F) %>%
inner_join(metadata) %>%
relocate(CageID:Background, .after = sample_id) %>%
column_to_rownames(var = "sample_id") %>%
select(-Genotype) %>%
rename(Genotype=Group2) %>%
arrange(Mode, Genotype, Sex, Age)
hm1 <- pheatmap(t(select(tar_merged_zscores, !Name:Background)), border_color = "NA",
annotation_col = data.frame(rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Age) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Sex) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Genotype) %>% column_to_rownames(var = "sample_id")
),
annotation_colors = ann_colors,
show_colnames = F, fontsize = 12, fontsize_row = 5.5, cluster_rows = T, cluster_cols = F, cutree_rows = 4, main = "Brain - lipids") #fontsize row 6 for lipids and 5.5 for targmet
save_pheatmap(hm1, "HM-FigureS6C.png", width=7.5, height=8.5)
#for plasma: use breaks to make tables more different
my.breaks <- c(seq(-1.5, 1.5, by=0.1))
hm4 <- pheatmap(t(select(tar_merged_zscores, !Name:Background)),
border_color = "NA",
color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(my.breaks)),
breaks = my.breaks,
annotation_col = data.frame(rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Age) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Sex) %>% column_to_rownames(var = "sample_id"),
rownames_to_column(tar_merged_zscores, var = "sample_id") %>% select(sample_id, Genotype) %>% column_to_rownames(var = "sample_id")
),
annotation_colors = ann_colors,
show_colnames = F, fontsize = 12, fontsize_row = 5.5, cluster_rows = T, cluster_cols = F, cutree_rows = 2, main = "Plasma - lipids") #fontsize row 6 for lipids and 5.5 for targmet
save_pheatmap(hm4, "HM-FigureS6A.png", width=7.5, height=8.5)
##############
##############Correlation analyses - 4 months##############
#Import DA metabolites
liver_da <- read.csv("Figure7A.csv")
plasma_da <- read.csv("Figure7C.csv")
brain_da <- read.csv("Figure7E.csv")
############## . WT brain-liver metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.1477
# Significance: 0.1481
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.190 0.255 0.309 0.370
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#135 @ rho>0.4, 41 @ rho>0.6, 8 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#270 lipids @ rho>0.4, 94 @ rho>0.6, 14 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_brain-v-liver.png", width=10, height=5)
#if rho<0.4, then make it grey (find clusters of lipids and proteins)
bk1 <- c(seq(-0.78, -0.39, by=0.02)) #spearman rho ranges from -0.39 to 0.49
bk2 <- c(seq(0.39, 0.76, by=0.02))
bk <- c(bk1,bk2) #combine the break limits for purpose of graphing
my_palette <- c(colorRampPalette(colors = c("#440154", "#2E6DA4"))(n = length(bk1)-1),
"gray60",
c(colorRampPalette(colors = c("#9DDE4B", "#FDE725"))(n = length(bk2)-1)))
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = my_palette, show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE, breaks = bk)
save_pheatmap(cor, "WT_brain-v-liver_grey.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_brain.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_liver.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.09216
# Significance: 0.7906
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.144 0.186 0.223 0.270
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD brain-liver metabolites##############
#5xFAD 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.07397
# Significance: 0.2767
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.182 0.242 0.293 0.351
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#35 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#49 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5x_brain-v-liver.png", width=10, height=5)
write.csv(colnames(correlation_df), "5x_brain.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5x_liver.csv", row.names=FALSE)
# #if rho<0.4, then make it grey (find clusters of lipids and proteins)
# bk1 <- c(seq(-0.78, -0.39, by=0.02)) #spearman rho ranges from -0.39 to 0.49
# bk2 <- c(seq(0.39, 0.76, by=0.02))
# bk <- c(bk1,bk2) #combine the break limits for purpose of graphing
#
# my_palette <- c(colorRampPalette(colors = c("#440154", "#2E6DA4"))(n = length(bk1)-1),
# "gray60",
# c(colorRampPalette(colors = c("#9DDE4B", "#FDE725"))(n = length(bk2)-1)))
#
# cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
# gaps_col = FALSE, color = my_palette, show_colnames = TRUE, show_rownames = TRUE,
# fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
# cutree_rows = 2, fontsize_row = 6, silent = FALSE, breaks = bk)
#
# save_pheatmap(cor, "WT_brain-v-liver_grey.png", width=10, height=5)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.01922
# Significance: 0.3929
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.212 0.284 0.342 0.418
# Permutation: free
# Number of permutations: 9999
############## . TREM2 brain-liver metabolites##############
#TREM2 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1132
# Significance: 0.7244
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.229 0.303 0.364 0.442
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#120 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#251 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 8, cutree_cols = 6, angle_col = 90,
cutree_rows = 6, fontsize_row = 8, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "T2_brain-v-liver.png", width=14, height=28)
write.csv(colnames(correlation_df), "T2_brain.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "T2_liver.csv", row.names=FALSE)
# #if rho<0.4, then make it grey (find clusters of lipids and proteins)
# bk1 <- c(seq(-0.78, -0.39, by=0.02)) #spearman rho ranges from -0.39 to 0.49
# bk2 <- c(seq(0.39, 0.76, by=0.02))
# bk <- c(bk1,bk2) #combine the break limits for purpose of graphing
#
# my_palette <- c(colorRampPalette(colors = c("#440154", "#2E6DA4"))(n = length(bk1)-1),
# "gray60",
# c(colorRampPalette(colors = c("#9DDE4B", "#FDE725"))(n = length(bk2)-1)))
#
# cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
# gaps_col = FALSE, color = my_palette, show_colnames = TRUE, show_rownames = TRUE,
# fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
# cutree_rows = 2, fontsize_row = 6, silent = FALSE, breaks = bk)
#
# save_pheatmap(cor, "WT_brain-v-liver_grey.png", width=10, height=5)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2175
# Significance: 0.8302
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.353 0.448 0.505 0.566
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD*TREM2 brain-liver metabolites##############
#5xFAD*TREM2 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2762
# Significance: 0.9121
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.392 0.504 0.561 0.619
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#104 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#239 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 8, cutree_cols = 6, angle_col = 90,
cutree_rows = 6, fontsize_row = 8, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5xT2_brain-v-liver.png", width=14, height=28)
write.csv(colnames(correlation_df), "5xT2_brain.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5xT2_liver.csv", row.names=FALSE)
# #if rho<0.4, then make it grey (find clusters of lipids and proteins)
# bk1 <- c(seq(-0.78, -0.39, by=0.02)) #spearman rho ranges from -0.39 to 0.49
# bk2 <- c(seq(0.39, 0.76, by=0.02))
# bk <- c(bk1,bk2) #combine the break limits for purpose of graphing
#
# my_palette <- c(colorRampPalette(colors = c("#440154", "#2E6DA4"))(n = length(bk1)-1),
# "gray60",
# c(colorRampPalette(colors = c("#9DDE4B", "#FDE725"))(n = length(bk2)-1)))
#
# cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
# gaps_col = FALSE, color = my_palette, show_colnames = TRUE, show_rownames = TRUE,
# fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
# cutree_rows = 2, fontsize_row = 6, silent = FALSE, breaks = bk)
#
# save_pheatmap(cor, "WT_brain-v-liver_grey.png", width=10, height=5)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2175
# Significance: 0.8302
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.353 0.448 0.505 0.566
# Permutation: free
# Number of permutations: 9999
############## . . Comparing brain-liver metabolites##############
WT_brain <- read.csv("WT_brain.csv")
WT_liver <- read.csv("WT_liver.csv")
Fad_brain <- read.csv("5x_brain.csv")
Fad_liver <- read.csv("5x_liver.csv")
T2_brain <- read.csv("T2_brain.csv")
T2_liver <- read.csv("T2_liver.csv")
FadT2_brain <- read.csv("5xT2_brain.csv")
FadT2_liver <- read.csv("5xT2_liver.csv")
#venn diagrams
library(ggVennDiagram)
genes <- paste("gene",1:1000,sep="")
set.seed(20231214)
# b <- list(WT=WT_brain$x,
# `5xFAD`=Fad_brain$x,
# Trem2=T2_brain$x,
# `5xFAD, Trem2`=FadT2_brain$x)
# ggVennDiagram(b) +#, set_size = 4) +
# scale_fill_gradient(low="grey90",high ="red")
# ggsave("venn_diagram_brain.png", width = 6, height = 5)
#
# l <- list(WT=WT_liver$x,
# `5xFAD`=Fad_liver$x,
# Trem2=T2_liver$x,
# `5xFAD, Trem2`=FadT2_liver$x)
# ggVennDiagram(l) +#, set_size = 4) +
# scale_fill_gradient(low="grey90",high ="red")
# ggsave("venn_diagram_liver.png", width = 6, height = 5)
#Combine brain-liver lists to show in one venn diagram
WT <- unique(c(WT_brain$x, WT_liver$x)) #20
Fad <- unique(c(Fad_brain$x, Fad_liver$x)) #74
T2 <- unique(c(T2_brain$x, T2_liver$x)) #267
FadT2 <- unique(c(FadT2_brain$x, FadT2_liver$x)) #268
bl <- list(WT=WT,
`5xFAD`=Fad,
Trem2=T2,
`5xFAD, Trem2`=FadT2)
ggVennDiagram(bl) +#, set_size = 4) +
scale_fill_gradient(low="grey90",high ="red")
ggsave("venn_diagram_brain-liver.png", width = 6, height = 5)
############## . WT brain-plasma metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.01563
# Significance: 0.44
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.284 0.368 0.451 0.537
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#79 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#64 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_brain-v-plasma.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_brain-v-plasma.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_plasma-v-brain.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.04355
# Significance: 0.3567
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.255 0.341 0.403 0.458
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD brain-plasma metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.04648
# Significance: 0.396
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.304 0.393 0.467 0.548
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#71 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#69 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5x_brain-v-plasma.png", width=10, height=5)
write.csv(colnames(correlation_df), "5x_brain-v-plasma.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5x_plasma-v-brain.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.04355
# Significance: 0.3567
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.255 0.341 0.403 0.458
# Permutation: free
# Number of permutations: 9999
############## . TREM2 brain-plasma metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.0722
# Significance: 0.651
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.218 0.287 0.340 0.461
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#131 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#258 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "T2_brain-v-plasma.png", width=10, height=5)
write.csv(colnames(correlation_df), "T2_brain-v-plasma.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "T2_plasma-v-brain.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.04783
# Significance: 0.5925
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.212 0.278 0.335 0.438
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD*TREM2 brain-plasma metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.4512
# Significance: 0.0492
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.358 0.448 0.526 0.591
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#135 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#263 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5xT2_brain-v-plasma.png", width=10, height=5)
write.csv(colnames(correlation_df), "5xT2_brain-v-plasma.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5xT2_plasma-v-brain.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.4399
# Significance: 0.0532
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.349 0.449 0.524 0.590
# Permutation: free
# Number of permutations: 9999
############## . . Comparing brain-plasma metabolites##############
WT_brain <- read.csv("WT_brain-v-plasma.csv")
WT_plasma <- read.csv("WT_plasma-v-brain.csv")
Fad_brain <- read.csv("5x_brain-v-plasma.csv")
Fad_plasma <- read.csv("5x_plasma-v-brain.csv")
T2_brain <- read.csv("T2_brain-v-plasma.csv")
T2_plasma <- read.csv("T2_plasma-v-brain.csv")
FadT2_brain <- read.csv("5xT2_brain-v-plasma.csv")
FadT2_plasma <- read.csv("5xT2_plasma-v-brain.csv")
#venn diagrams
library(ggVennDiagram)
genes <- paste("gene",1:1000,sep="")
set.seed(20231214)
#Combine brain-liver lists to show in one venn diagram
WT <- unique(c(WT_brain$x, WT_plasma$x)) #124
Fad <- unique(c(Fad_brain$x, Fad_plasma$x)) #121
T2 <- unique(c(T2_brain$x, T2_plasma$x)) #280
FadT2 <- unique(c(FadT2_brain$x, FadT2_plasma$x)) #293
bp <- list(WT=WT,
`5xFAD`=Fad,
Trem2=T2,
`5xFAD, Trem2`=FadT2)
ggVennDiagram(bp) +#, set_size = 4) +
scale_fill_gradient(low="grey90",high ="red")
ggsave("venn_diagram_brain-plasma.png", width = 6, height = 5)
############## . . Extracting 5xFAD, TREM2 (only) brain-plasma metabolites##############
WT_brain <- read.csv("WT_brain-v-plasma.csv")
WT_plasma <- read.csv("WT_plasma-v-brain.csv")
Fad_brain <- read.csv("5x_brain-v-plasma.csv")
Fad_plasma <- read.csv("5x_plasma-v-brain.csv")
T2_brain <- read.csv("T2_brain-v-plasma.csv")
T2_plasma <- read.csv("T2_plasma-v-brain.csv")
FadT2_brain <- read.csv("5xT2_brain-v-plasma.csv")
FadT2_plasma <- read.csv("5xT2_plasma-v-brain.csv")
WT <- unique(c(WT_brain$x, WT_plasma$x)) #124
Fad <- unique(c(Fad_brain$x, Fad_plasma$x)) #121
T2 <- unique(c(T2_brain$x, T2_plasma$x)) #280
FadT2 <- unique(c(FadT2_brain$x, FadT2_plasma$x)) #293
all_exclusions <- c(WT, Fad, T2)
FadT2_only <- FadT2[!FadT2 %in% all_exclusions]
write.csv(FadT2_only, "5xT2_only_brain-v-plasma.csv", row.names=FALSE)
############## . WT plasma-liver metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep plasma data from Mouse IDs that we have liver data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have plasma data for
liver_df_cor <- liver_df[liver_df$Name %in% plasma_df$Name, ]
#Reassign rownames
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.106
# Significance: 0.293
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.297 0.390 0.471 0.548
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(plasma_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep plasma analytes (columns) with at least one rho value > threshold
#109 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#176 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_plasma-v-liver.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_plasma-v-liver.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_liver-v-plasma.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
plasma_subset <- plasma_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
plasma.dist <- vegdist(plasma_subset, method="euclidean")
liver.dist <- vegdist(liver_subset, method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1723
# Significance: 0.7679
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.292 0.382 0.463 0.550
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD plasma-liver metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep plasma data from Mouse IDs that we have liver data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have plasma data for
liver_df_cor <- liver_df[liver_df$Name %in% plasma_df$Name, ]
#Reassign rownames
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2044
# Significance: 0.8402
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.276 0.363 0.434 0.533
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(plasma_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep plasma analytes (columns) with at least one rho value > threshold
#126 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#219 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5x_plasma-v-liver.png", width=10, height=5)
write.csv(colnames(correlation_df), "5x_plasma-v-liver.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5x_liver-v-plasma.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
plasma_subset <- plasma_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
plasma.dist <- vegdist(plasma_subset, method="euclidean")
liver.dist <- vegdist(liver_subset, method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.1992
# Significance: 0.1571
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.259 0.354 0.436 0.524
# Permutation: free
# Number of permutations: 9999
############## . TREM2 plasma-liver metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep plasma data from Mouse IDs that we have liver data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have plasma data for
liver_df_cor <- liver_df[liver_df$Name %in% plasma_df$Name, ]
#Reassign rownames
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1271
# Significance: 0.7674
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.200 0.259 0.352 0.503
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(plasma_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep plasma analytes (columns) with at least one rho value > threshold
#254 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#269 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "T2_plasma-v-liver.png", width=10, height=5)
write.csv(colnames(correlation_df), "T2_plasma-v-liver.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "T2_liver-v-plasma.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
plasma_subset <- plasma_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
plasma.dist <- vegdist(plasma_subset, method="euclidean")
liver.dist <- vegdist(liver_subset, method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1227
# Significance: 0.767
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.205 0.269 0.353 0.508
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD*TREM2 plasma-liver metabolites##############
#WT 4 months targ.met
TREM2_tar_df = tar_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep plasma data from Mouse IDs that we have liver data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have plasma data for
liver_df_cor <- liver_df[liver_df$Name %in% plasma_df$Name, ]
#Reassign rownames
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.09592
# Significance: 0.3275
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.305 0.383 0.488 0.589
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(plasma_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep plasma analytes (columns) with at least one rho value > threshold
#250 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#266 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5xT2_plasma-v-liver.png", width=10, height=5)
write.csv(colnames(correlation_df), "5xT2_plasma-v-liver.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5xT2_liver-v-plasma.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
plasma_subset <- plasma_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
plasma.dist <- vegdist(plasma_subset, method="euclidean")
liver.dist <- vegdist(liver_subset, method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.1884
# Significance: 0.1885
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.275 0.346 0.400 0.455
# Permutation: free
# Number of permutations: 9999
############## . . Comparing plasma-liver metabolites##############
WT_liver <- read.csv("WT_liver-v-plasma.csv")
WT_plasma <- read.csv("WT_plasma-v-liver.csv")
Fad_liver <- read.csv("5x_liver-v-plasma.csv")
Fad_plasma <- read.csv("5x_plasma-v-liver.csv")
T2_liver <- read.csv("T2_liver-v-plasma.csv")
T2_plasma <- read.csv("T2_plasma-v-liver.csv")
FadT2_liver <- read.csv("5xT2_liver-v-plasma.csv")
FadT2_plasma <- read.csv("5xT2_plasma-v-liver.csv")
#venn diagrams
library(ggVennDiagram)
genes <- paste("gene",1:1000,sep="")
set.seed(20231214)
#Combine brain-liver lists to show in one venn diagram
WT <- unique(c(WT_liver$x, WT_plasma$x)) #229
Fad <- unique(c(Fad_liver$x, Fad_plasma$x)) #255
T2 <- unique(c(T2_liver$x, T2_plasma$x)) #314
FadT2 <- unique(c(FadT2_liver$x, FadT2_plasma$x)) #308
pl <- list(WT=WT,
`5xFAD`=Fad,
Trem2=T2,
`5xFAD, Trem2`=FadT2)
ggVennDiagram(pl) +#, set_size = 4) +
scale_fill_gradient(low="grey90",high ="red")
ggsave("venn_diagram_plasma-liver.png", width = 6, height = 5)
############## . WT brain-liver lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2867
# Significance: 0.9936
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.188 0.247 0.307 0.375
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#81 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#135 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_brain-v-liver-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_brain-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_liver-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1578
# Significance: 0.8656
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.187 0.255 0.307 0.360
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD brain-liver lipids##############
#5xFAD 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.07266
# Significance: 0.2829
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.191 0.252 0.305 0.362
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#243 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#227 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5x_brain-v-liver-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "5x_brain-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5x_liver-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.1887
# Significance: 0.1124
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.201 0.262 0.314 0.371
# Permutation: free
# Number of permutations: 9999
############## . TREM2 brain-liver lipids##############
#TREM2 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.05507
# Significance: 0.5745
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.278 0.350 0.422 0.508
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#495 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#586 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 8, cutree_cols = 6, angle_col = 90,
cutree_rows = 6, fontsize_row = 8, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "T2_brain-v-liver-lip.png", width=14, height=28)
write.csv(colnames(correlation_df), "T2_brain-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "T2_liver-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.07839
# Significance: 0.5926
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.273 0.361 0.448 0.528
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD*TREM2 brain-liver lipids##############
#5xFAD*TREM2 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep brain data from Mouse IDs that we have liver data for
brain_df_cor <- brain_df[brain_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have brain data for
liver_df_cor <- liver_df[liver_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="bray")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2273
# Significance: 0.816
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.339 0.434 0.520 0.586
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#494 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#703 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 8, cutree_cols = 6, angle_col = 90,
cutree_rows = 6, fontsize_row = 8, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5xT2_brain-v-liver-lip.png", width=14, height=28)
write.csv(colnames(correlation_df), "5xT2_brain-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5xT2_liver-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
liver_subset <- liver_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="bray")
liver.dist <- vegdist(liver_subset, method="bray")
mantel(brain.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2268
# Significance: 0.8229
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.338 0.438 0.522 0.591
# Permutation: free
# Number of permutations: 9999
############## . . Comparing brain-liver lipids##############
WT_brain <- read.csv("WT_brain-lip.csv")
WT_liver <- read.csv("WT_liver-lip.csv")
Fad_brain <- read.csv("5x_brain-lip.csv")
Fad_liver <- read.csv("5x_liver-lip.csv")
T2_brain <- read.csv("T2_brain-lip.csv")
T2_liver <- read.csv("T2_liver-lip.csv")
FadT2_brain <- read.csv("5xT2_brain-lip.csv")
FadT2_liver <- read.csv("5xT2_liver-lip.csv")
#venn diagrams
library(ggVennDiagram)
genes <- paste("gene",1:1000,sep="")
set.seed(20231214)
#Combine brain-liver lists to show in one venn diagram
WT <- unique(c(WT_brain$x, WT_liver$x)) #205
Fad <- unique(c(Fad_brain$x, Fad_liver$x)) #411
T2 <- unique(c(T2_brain$x, T2_liver$x)) #750
FadT2 <- unique(c(FadT2_brain$x, FadT2_liver$x)) #803
bl <- list(WT=WT,
`5xFAD`=Fad,
Trem2=T2,
`5xFAD, Trem2`=FadT2)
ggVennDiagram(bl) +#, set_size = 4) +
scale_fill_gradient(low="grey90",high ="red")
ggsave("venn_diagram_brain-liver-lip.png", width = 6, height = 5)
############## . WT brain-plasma lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.2525
# Significance: 0.0815
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.229 0.313 0.370 0.437
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#264 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#352 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_brain-v-plasma-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_brain-v-plasma-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_plasma-v-brain-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.4273
# Significance: 0.0157 *
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.232 0.308 0.382 0.464
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD brain-plasma lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.2319
# Significance: 0.9536
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.209 0.280 0.341 0.403
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#248 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#176 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5x_brain-v-plasma-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "5x_brain-v-plasma-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5x_plasma-v-brain-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1521
# Significance: 0.8285
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.219 0.289 0.346 0.408
# Permutation: free
# Number of permutations: 9999
############## . TREM2 brain-plasma lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.1145
# Significance: 0.688
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.282 0.358 0.426 0.488
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#474 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#652 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "T2_brain-v-plasma-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "T2_brain-v-plasma-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "T2_plasma-v-brain-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.003953
# Significance: 0.4918
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.270 0.345 0.404 0.474
# Permutation: free
# Number of permutations: 9999
############## . 5xFAD*TREM2 brain-plasma lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("HEMI 5xFAD_TREM2"), Age == 4) %>%
column_to_rownames(var = "sample_id")
brain_df = TREM2_tar_df %>%
filter(Mode == "brain")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
#Only keep brain data from Mouse IDs that we have plasma data for
brain_df_cor <- brain_df[brain_df$Name %in% plasma_df$Name, ]
#Only keep plasma data from Mouse IDs that we have brain data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% brain_df$Name, ]
#Reassign rownames
rownames(brain_df_cor) <- NULL
rownames(brain_df_cor) <- brain_df_cor$Name
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
#Arrange metadata
brain_df_cor <- brain_df_cor %>%
arrange(Name)
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
brain.dist <- vegdist(brain_df_cor[,-(1:11)], method="euclidean")
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.3506
# Significance: 0.894
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.439 0.547 0.657 0.751
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(brain_df_cor[,-(1:11)], plasma_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep brain analytes (columns) with at least one rho value > threshold
#505 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep plasma analytes (rows) with at least one rho value > threshold
#912 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "5xT2_brain-v-plasma-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "5xT2_brain-v-plasma-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "5xT2_plasma-v-brain-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
brain_subset <- brain_df_cor[,colnames(correlation_df)]
plasma_subset <- plasma_df_cor[,rownames(correlation_df)]
brain.dist <- vegdist(brain_subset, method="euclidean")
plasma.dist <- vegdist(plasma_subset, method="euclidean")
mantel(brain.dist, plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = brain.dist, ydis = plasma.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: -0.3514
# Significance: 0.8982
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.452 0.573 0.661 0.750
# Permutation: free
# Number of permutations: 9999
############## . . Comparing brain-plasma lipids##############
WT_brain <- read.csv("WT_brain-v-plasma-lip.csv")
WT_plasma <- read.csv("WT_plasma-v-brain-lip.csv")
Fad_brain <- read.csv("5x_brain-v-plasma-lip.csv")
Fad_plasma <- read.csv("5x_plasma-v-brain-lip.csv")
T2_brain <- read.csv("T2_brain-v-plasma-lip.csv")
T2_plasma <- read.csv("T2_plasma-v-brain-lip.csv")
FadT2_brain <- read.csv("5xT2_brain-v-plasma-lip.csv")
FadT2_plasma <- read.csv("5xT2_plasma-v-brain-lip.csv")
#venn diagrams
library(ggVennDiagram)
genes <- paste("gene",1:1000,sep="")
set.seed(20231214)
#Combine brain-liver lists to show in one venn diagram
WT <- unique(c(WT_brain$x, WT_plasma$x)) #528
Fad <- unique(c(Fad_brain$x, Fad_plasma$x)) #389
T2 <- unique(c(T2_brain$x, T2_plasma$x)) #815
FadT2 <- unique(c(FadT2_brain$x, FadT2_plasma$x)) #932
bp <- list(WT=WT,
`5xFAD`=Fad,
Trem2=T2,
`5xFAD, Trem2`=FadT2)
ggVennDiagram(bp) +#, set_size = 4) +
scale_fill_gradient(low="grey90",high ="red")
ggsave("venn_diagram_brain-plasma-lip.png", width = 6, height = 5)
############## . . Extracting 5xFAD, TREM2 (only) brain-plasma lipids##############
WT_brain <- read.csv("WT_brain-v-plasma-lip.csv")
WT_plasma <- read.csv("WT_plasma-v-brain-lip.csv")
Fad_brain <- read.csv("5x_brain-v-plasma-lip.csv")
Fad_plasma <- read.csv("5x_plasma-v-brain-lip.csv")
T2_brain <- read.csv("T2_brain-v-plasma-lip.csv")
T2_plasma <- read.csv("T2_plasma-v-brain-lip.csv")
FadT2_brain <- read.csv("5xT2_brain-v-plasma-lip.csv")
FadT2_plasma <- read.csv("5xT2_plasma-v-brain-lip.csv")
WT <- unique(c(WT_brain$x, WT_plasma$x)) #528
Fad <- unique(c(Fad_brain$x, Fad_plasma$x)) #389
T2 <- unique(c(T2_brain$x, T2_plasma$x)) #815
FadT2 <- unique(c(FadT2_brain$x, FadT2_plasma$x)) #932
all_exclusions <- c(WT, Fad, T2)
FadT2_only <- FadT2[!FadT2 %in% all_exclusions]
write.csv(FadT2_only, "5xT2_only_brain-v-plasma-lip.csv", row.names=FALSE)
############## . WT plasma-liver lipids##############
#WT 4 months targ.met
TREM2_tar_df = lip_df_combined %>%
filter(Group2 %in% c("WT 5xFAD"), Age == 4) %>%
column_to_rownames(var = "sample_id")
plasma_df = TREM2_tar_df %>%
filter(Mode == "plasma")
liver_df = TREM2_tar_df %>%
filter(Mode == "liver")
#Only keep plasma data from Mouse IDs that we have liver data for
plasma_df_cor <- plasma_df[plasma_df$Name %in% liver_df$Name, ]
#Only keep liver data from Mouse IDs that we have plasma data for
liver_df_cor <- liver_df[liver_df$Name %in% plasma_df$Name, ]
#Reassign rownames
rownames(plasma_df_cor) <- NULL
rownames(plasma_df_cor) <- plasma_df_cor$Name
rownames(liver_df_cor) <- NULL
rownames(liver_df_cor) <- liver_df_cor$Name
#Arrange metadata
plasma_df_cor <- plasma_df_cor %>%
arrange(Name)
liver_df_cor <- liver_df_cor %>%
arrange(Name)
#Mantel significance test
set.seed <- 0123456789
plasma.dist <- vegdist(plasma_df_cor[,-(1:11)], method="euclidean")
liver.dist <- vegdist(liver_df_cor[,-(1:11)], method="euclidean")
mantel(plasma.dist, liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
# Call:
# mantel(xdis = plasma.dist, ydis = liver.dist, method = "spearman", permutations = 9999, na.rm = TRUE)
#
# Mantel statistic r: 0.01878
# Significance: 0.4128
#
# Upper quantiles of permutations (null model):
# 90% 95% 97.5% 99%
# 0.247 0.325 0.400 0.488
# Permutation: free
# Number of permutations: 9999
#correlation analysis
correlation <- t(cor(plasma_df_cor[,-(1:11)], liver_df_cor[,-(1:11)], method = c('spearman')))
correlation[is.na(correlation)] <- 0
#only keep plasma analytes (columns) with at least one rho value > threshold
#714 @ rho>0.7
correlation_df <- as.data.frame(correlation)
columns_to_keep <- sapply(correlation_df, function(x) {
max_val <- max(x, na.rm = TRUE)
min_val <- min(x, na.rm = TRUE)
return(max_val > 0.7 | min_val < -0.7)
})
correlation_df <- correlation_df[, columns_to_keep]
#only keep liver analytes (rows) with at least one rho value > threshold
#449 @ rho>0.7
correlation_df <- correlation_df[apply(correlation_df, 1, max) > 0.7 |
apply(correlation_df, 1, min) < -0.7, ]
#heatmap w/ thresholds
correlation <- as.matrix(correlation_df)
cor <- pheatmap(correlation, cluster_cols = TRUE, clustering_method = 'average', cluster_rows = TRUE,
gaps_col = FALSE, color = viridis(10), show_colnames = TRUE, show_rownames = TRUE,
fontsize = 8, fontsize_col = 6, cutree_cols = 2, angle_col = 90,
cutree_rows = 2, fontsize_row = 6, silent = FALSE)#, breaks = my.breaks)
save_pheatmap(cor, "WT_plasma-v-liver-lip.png", width=10, height=5)
write.csv(colnames(correlation_df), "WT_plasma-v-liver-lip.csv", row.names=FALSE)
write.csv(rownames(correlation_df), "WT_liver-v-plasma-lip.csv", row.names=FALSE)
#Mantel significance test (analytes in heatmap only)
set.seed <- 0123456789
plasma_subse

TREM2-metabolomics-figures-final-v3.R at commit c38d0be, under MIT · at the source

Overview

Authors: Gina Faraci^1,2,3, Benjamin Goodfriend¹, Joseph Bishop¹, Michael Vu¹, Julio Avelar-Barragan², Sage J.B. Dunham², Jason A. Rothman⁴, Katrine L. Whiteson², Amrita K. Cheema⁵, Giedre Milinkeviciute³, Andrea J. Tenner^2,3,6,7, Frank M. LaFerla^3,6, Grant R. MacGregor^8,9, Kim N. Green^3,6, Mark Mapstone^1,3

ORCID iDs: Gina Faraci, Jason A. Rothman, Katrine L. Whiteson

Department of Neurology, University of California, Irvine, CA 92617, USA
Department of Molecular Biology and Biochemistry, University of California, Irvine, CA 92697, USA
Institute for Memory Impairments and Neurological Disorders, University of California, Irvine, CA 92697, USA
Department of Microbiology and Plant Pathology, University of California, Riverside, CA 92521, USA
Department of Oncology, Georgetown University, Washington, DC 20007, USA
Department of Neurobiology and Behavior, University of California, Irvine, CA 92697, USA
Department of Pathology and Laboratory Medicine, University of California, Irvine, CA 92617, USA
Transgenic Mouse Facility, ULAR, Office of Research, University of California, Irvine, CA 92697, USA
Department of Developmental and Cell Biology, University of California, Irvine, CA 92697, USA

Institutions: University of California, Irvine (United States); University of California, Riverside (United States); Georgetown University (United States)

Journal: Neurobiology of aging, volume 165, pages 24-37

Dates: published online 29 April 2026; in print September 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.1016/j.neurobiolaging.2026.04.007 · PMID 42102578 · PMCID PMC13249471 · OpenAlex W7160013979

Open access: hybrid, a free copy (OpenAlex)

Status: code verified

Categories: genetics / omics (modality), human (organism), mouse (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)

Methods: Statistics, Connectivity, Machine learning

Keywords: Alzheimer’s disease, Metabolomics, Lipid, Mouse, Brain, Plasma, Liver

MeSH: Alzheimer Disease*, Brain*, Genetic Association Studies*, Genetic Variation*, Homeostasis*, Liver*, Membrane Glycoproteins*, Receptors, Immunologic*, Animals, Disease Models, Animal, Female, Humans, Lipid Metabolism, Male, Mice, Transgenic, Triglycerides (* major topic)

Topic: Neuroinflammation and Neurodegeneration Mechanisms (Neurology, Neuroscience), according to OpenAlex

Funding: NCI NIH HHS (P30 CA062203); University of California, Irvine; National Institute on Aging (U54 AG054349); National Institutes of Health; NIA NIH HHS (U54 AG054349); National Institutes of Health National Cancer Institute (P30CA062203)

Citations: not cited yet (Europe PMC); 70 references in the paper

Abstract

The human Triggering Receptor Expressed on Myeloid cells 2 (TREM2) gene is expressed predominantly by microglia in the brain and the R47H coding variant of TREM2 is associated with increased risk for late-onset Alzheimer’s disease (LOAD). We performed lipidomic and metabolomic analysis of liver, plasma, and brain in 4- and 12-month-old Trem2R47H homozygous (n = 27), 5xFAD hemizygous (n = 63), 5xFAD hemizygous, Trem2R47H homozygous (n = 25), and wild type (n = 65) mice. Lipid and metabolite abundances differed significantly across tissue types with the most differences seen in the liver and plasma of Trem2R47H mice at the 4-month timepoint. Cross-tissue correlation analyses revealed increased metabolic crosstalk along the liver-plasma-brain axis in Trem2R47H mice. Plasma triacylglyceride levels were significantly lower in females compared to males regardless of genotype, and 5-methyltetrahydrofolic acid levels were elevated in the brains of animals homozygous for the Trem2R47H variant. Together, these findings demonstrate early dyshomeostasis of the liver-plasma-brain axis of Trem2R47H mice which impacts several key metabolic pathways involving lipids, cellular energy metabolism, and brain folate metabolism.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

GinaFaraci/TREM2-Metabolomics

License: MIT

State: the link answers, verified on 30 September 2026

Evidence: files inventoried

Commit: c38d0be42ec0194add0c2f8af1f784ea86dfd452, 16 March 2026

Languages: R (1)

Size: 10 files, 1 script

Software Heritage: not archived

Found in: “Data availability”

Holds: README, license file

Not found: CITATION.cff, environment file, tests, continuous integration, documentation

Tools: broom (1 file), ggpubr (1 file), pheatmap (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 30 September 2026: the link answers

30 September 2026: the link answers

3 files

TREM2-metabolomics-figures-final-v3.R, R, 4,577 lines, 2 matches
LICENSE, License, 21 lines
README.md, Text, 5 lines

At the source: github.com/GinaFaraci/TREM2-Metabolomics

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

doi:10.7910/dvn/xzzcde, at the source; found in “Data availability”

Data availability

The datasets generated and/or analyzed during the current study are publicly available via Harvard Dataverse (https://doi.org/10.7910/DVN/XZZCDE). An identical version of these materials, including analysis code, is also available at GitHub (https://github.com/GinaFaraci/TREM2-Metabolomics).

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 15 authors, 7 keywords, 16 MeSH terms, 6 funders, 68 references.

Cite

This paper

Faraci, G., Goodfriend, B., Bishop, J., Vu, M., Avelar-Barragan, J., Dunham, S. J., Rothman, J. A., Whiteson, K. L., Cheema, A. K., Milinkeviciute, G., Tenner, A. J., LaFerla, F. M., MacGregor, G. R., Green, K. N., & Mapstone, M. (2026). The TREM2 R47H variant is associated with liver-plasma-brain axis dyshomeostasis in the 5xFAD mouse model of Alzheimer's disease. Neurobiology of aging, 165, 24-37. https://doi.org/10.1016/j.neurobiolaging.2026.04.007

BibTeX

@article{faraci2026trem2, author = {Faraci, Gina and Goodfriend, Benjamin and Bishop, Joseph and Vu, Michael and Avelar-Barragan, Julio and Dunham, Sage J.B. and Rothman, Jason A. and Whiteson, Katrine L. and Cheema, Amrita K. and Milinkeviciute, Giedre and Tenner, Andrea J. and LaFerla, Frank M. and MacGregor, Grant R. and Green, Kim N. and Mapstone, Mark}, title = {{The TREM2 R47H variant is associated with liver-plasma-brain axis dyshomeostasis in the 5xFAD mouse model of Alzheimer's disease}}, journal = {Neurobiology of aging}, year = {2026}, month = apr, volume = {165}, pages = {24--37}, publisher = {Elsevier BV}, issn = {0197-4580}, doi = {10.1016/j.neurobiolaging.2026.04.007}, url = {https://doi.org/10.1016/j.neurobiolaging.2026.04.007}, pmid = {42102578}, pmcid = {PMC13249471} }

RIS

TY - JOUR AU - Faraci, Gina AU - Goodfriend, Benjamin AU - Bishop, Joseph AU - Vu, Michael AU - Avelar-Barragan, Julio AU - Dunham, Sage J.B. AU - Rothman, Jason A. AU - Whiteson, Katrine L. AU - Cheema, Amrita K. AU - Milinkeviciute, Giedre AU - Tenner, Andrea J. AU - LaFerla, Frank M. AU - MacGregor, Grant R. AU - Green, Kim N. AU - Mapstone, Mark TI - The TREM2 R47H variant is associated with liver-plasma-brain axis dyshomeostasis in the 5xFAD mouse model of Alzheimer's disease T2 - Neurobiology of aging J2 - Neurobiol Aging PY - 2026 DA - 2026/04/29 VL - 165 SP - 24 EP - 37 SN - 0197-4580 PB - Elsevier BV DO - 10.1016/j.neurobiolaging.2026.04.007 UR - https://doi.org/10.1016/j.neurobiolaging.2026.04.007 LA - en ER -

CSL-JSON

{ "id": "10.1016/j.neurobiolaging.2026.04.007", "type": "article-journal", "title": "The TREM2 R47H variant is associated with liver-plasma-brain axis dyshomeostasis in the 5xFAD mouse model of Alzheimer's disease", "container-title": "Neurobiology of aging", "author": [ { "family": "Faraci", "given": "Gina" }, { "family": "Goodfriend", "given": "Benjamin" }, { "family": "Bishop", "given": "Joseph" }, { "family": "Vu", "given": "Michael" }, { "family": "Avelar-Barragan", "given": "Julio" }, { "family": "Dunham", "given": "Sage J.B." }, { "family": "Rothman", "given": "Jason A." }, { "family": "Whiteson", "given": "Katrine L." }, { "family": "Cheema", "given": "Amrita K." }, { "family": "Milinkeviciute", "given": "Giedre" }, { "family": "Tenner", "given": "Andrea J." }, { "family": "LaFerla", "given": "Frank M." }, { "family": "MacGregor", "given": "Grant R." }, { "family": "Green", "given": "Kim N." }, { "family": "Mapstone", "given": "Mark" } ], "container-title-short": "Neurobiol Aging", "volume": "165", "page": "24-37", "DOI": "10.1016/j.neurobiolaging.2026.04.007", "PMID": "42102578", "PMCID": "PMC13249471", "ISSN": "0197-4580", "publisher": "Elsevier BV", "URL": "https://doi.org/10.1016/j.neurobiolaging.2026.04.007", "language": "en", "issued": { "date-parts": [ [ 2026, 4, 29 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1186/s40478-026-02287-x [code]: Cellular signatures of melanocortin pathway genes across the locus coeruleus.
Journal: Acta neuropathologica communications
In common: tidyverse, Alzheimer's / dementia, genetics / omics, mouse, 1 other category, 5 references
[2] doi:10.1186/s12974-026-03848-6 [code]: Iba1 deficiency impairs microglial synaptic remodeling and neuronal survival after axonal injury.
Journal: Journal of neuroinflammation
In common: pheatmap, tidyverse, genetics / omics, mouse, cellular / molecular, 4 references
[3] doi:10.1007/s00401-026-03036-z: Mechanisms of increased Alzheimer's disease pathology with R47H and R62H TREM2 variants.
Journal: Acta neuropathologica
In common: Alzheimer's / dementia, genetics / omics, cellular / molecular, 5 references
[4] doi:10.1038/s44400-026-00100-z [code]: Peripheral blood microarray-based transcriptomic and epigenetic analyses identify immune, inflammation, and metabolic dysregulation in Alzheimer's disease.
Journal: NPJ dementia
In common: broom, pheatmap, ggpubr, 1 other tool, Alzheimer's / dementia, genetics / omics, 1 reference
[5] doi:10.1093/bioinformatics/btag592 [code]: Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: broom, pheatmap, ggpubr, 1 other tool, genetics / omics, 1 reference
[6] doi:10.1016/j.isci.2026.115638 [code]: Genome-scale metabolic modeling uncovers cell-type specific signatures associated with APOE variants.
Journal: iScience
In common: tidyverse, Alzheimer's / dementia, genetics / omics, cellular / molecular, 3 references
[7] doi:10.1016/j.cell.2026.05.047 [code]: An emergent disease-associated motor neuron state precedes cell death in ALS.
Journal: Cell
In common: broom, pheatmap, ggpubr, 1 other tool, genetics / omics, mouse, cellular / molecular
[8] doi:10.1016/j.celrep.2026.117073 [code]: Single-cell epigenomics uncovers heterochromatin instability and transcription factor dysfunction during mouse brain aging.
Journal: Cell reports
In common: broom, pheatmap, ggpubr, 1 other tool, genetics / omics, mouse, cellular / molecular
[9] doi:10.1186/s12974-026-03897-x [code]: Ischemic injury triggers a protective microglial phenotype in models of Aβ pathology.
Journal: Journal of neuroinflammation
In common: pheatmap, ggpubr, tidyverse, Alzheimer's / dementia, mouse, cellular / molecular, 1 reference
[10] doi:10.1007/s12035-026-05859-z [code]: Uncovering Spatiotemporal and Functional Dynamics of Long Non-coding RNAs During Alzheimer's Progression in the Human Brain at Single-Cell Resolution.
Journal: Molecular neurobiology
In common: pheatmap, ggpubr, tidyverse, Alzheimer's / dementia, genetics / omics, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.