Rodent fMRI‐ Functional Connectivity Between Networks or Regions - CoBrALab/documentation GitHub Wiki
Disclaimer: This wiki was created my Marina Broomfield. Please cite this page and give credit when using or copying the following code.
Functional connectivity allows us to examine the time series of the desires network/region and see whether the activation patter of our zone of interest is correlated, anti correlated, or completely independent of other brain region's brain activity. An example of this could be though of as a question. What happens to the somatosensory network activity when the thalamic network is activated? When the thalamic network is activated, this might cause the somatosensory network to be active at the same time or it might cause it be suppressed, finally it could also be the case that the thalamic network activity does not impact the somatosensory network at all. Functional connectivity analysis allows us to visualize these patterns through each region's activity over time.
With functional connectivity analysis we can also see whether there are differences across each of our groups and examine how the activity of the brain might be affected by each condition.
First it is important to obtain the cleaned dual regression time course of each subject for all networks/regions of interest. Create a single .csv file for each subject with each columns representing networks/regions and rows the activity across time. A csv file per subject containing the activity of each voxel over time for all Independent components is generated by RABIES and can be found in:
/diagnosis/diagnosis_QC_females/analysis_datasink/dual_regression_timecourse_csv/_split_name_sub-C01S015_ses-1_task-rest_acq-EPI_bold/
Make sure you only extract the data from the desired ICs. A sample code on how to do this in bash can be seen below:
#!/bin/bash
# Define the source pattern and target directories
#files="/data/scratch2/marina/FINALexp/scanning/functional/diagnosis/diagnosis_all_females/analysis_datasink/dual_regression_timecourse_csv/*/*.csv"
tgt_dir="/data/scratch2/marina/FINALexp/functional_correlations/DR_QC/QC_females"
# Iterate over each CSV file matching the pattern
for file in /data/scratch2/marina/FINALexp/scanning/functional/diagnosis/diagnosis_QC_females/analysis_datasink/dual_regression_timecourse_csv/*/*.csv; do
# Extract the relative path of the file
rel_path="${file#$src_dir/}"
# Extract the directory path and create it in the target directory
dir_path=$(dirname "$rel_path")
mkdir -p "$tgt_dir/$dir_path"
# Extract the filename and add the 'merged_' prefix
filename=$(basename "$file")
new_filename="merged_$filename.csv"
# Extract the desired columns and save to the new location with the new filename
awk -F, 'BEGIN {OFS=","} {print $1, $2, $3, $9}' "$file" > "$tgt_dir/$new_filename"
done
Following this step, the fisher z correlation across each network/region for each group is calculated. Once the correlation values are obtained a t-test is performed to compare differences in correlations across groups. A Benjamini Hochberg correction is then applied to correct for multiple comparisons.
Disclaimer: The following code was created by Dr. Yohan Yee & Marina Broomfield. Please cite if used.
library(tidyverse)
library(psych)
library(corrplot)
library(pheatmap)
library(ggnewscale)
library(viridis)
library('colorspace')
# Load the data
make_cor_data <- function(f, s) {
out <- read_csv(f, col_names = F) %>%
cor() %>%
fisherz() %>%
as_tibble() %>%
mutate(seed_structure=colnames(.)) %>%
gather(target_structure, z, starts_with("X")) %>%
mutate(file=basename(f), sex=s)
return(out)
}
df_females <- list.files("/data/scratch2/marina/FINALexp/functional_correlations/SBC_QC/seed_females_QC", full.names = T) %>%
map_dfr(make_cor_data, s="F")
df_males <- list.files("/data/scratch2/marina/FINALexp/functional_correlations/SBC_QC/seed_males_QC", full.names = T) %>%
map_dfr(make_cor_data, s="M")
df <- bind_rows(df_females, df_males)
#####################################################
# Compute the average
df_mean <- df %>%
group_by(seed_structure, target_structure, sex) %>%
summarize(mean_z=mean(z)) %>%
ungroup() %>%
mutate(mean_r=fisherz2r(mean_z))
zmat_females <- df_mean %>%
filter(sex=="F") %>%
select(seed_structure, target_structure, value=mean_z) %>%
spread(target_structure, value) %>%
select(-seed_structure) %>%
as.matrix() %>%
`rownames<-`(colnames(.))
zmat_males <- df_mean %>%
filter(sex=="M") %>%
select(seed_structure, target_structure, value=mean_z) %>%
spread(target_structure, value) %>%
select(-seed_structure) %>%
as.matrix() %>%
`rownames<-`(colnames(.))
zmat <- matrix(0, nrow=12, ncol = 12, dimnames = list(rownames(zmat_females), colnames(zmat_females)))
zmat[upper.tri(zmat)] <- zmat_females[upper.tri(zmat_females)]
zmat[lower.tri(zmat)] <- zmat_males[lower.tri(zmat_males)]
#diag(zmat) <- 1
corrplot(zmat, method = "color", is.corr = F)
pheatmap(zmat, cluster_rows = F, cluster_cols = F)
df_mean %>%
ggplot(aes(x=seed_structure, y=target_structure)) +
geom_tile(aes(fill=mean_z)) +
coord_fixed(ratio=1) +
facet_wrap(~ sex) +
theme_bw()
strucs <- paste("X", 1:12, sep="")
strucs
# Create a named vector for the new labels
#DR
#new_labels <- c("X1" = "SM", "X2" = "SMS", "X3" = "THAL", "X4" = "STR")
new_labels=c("X1"="dSTR","X2"="FAC","X3"="ICC","X4"="INSU","X5"="MC1",
"X6"="MC2","X7"="RPC","X8"="SMS","X9"="SMS2","X10"="THAL",
"X11"="vHIP","X12"="NAC")
df_mean %>%
mutate(seed_order = as.integer(str_replace(seed_structure, "X", "")),
target_order = as.integer(str_replace(target_structure, "X", ""))) %>%
arrange(seed_order, target_order) %>%
select(-mean_r) %>%
spread(sex, mean_z) %>%
mutate(plot_value = case_when(target_order > seed_order ~ `M`,
seed_order > target_order ~ `F`,
TRUE ~ 0),
plot_sex = case_when(target_order > seed_order ~ "M",
seed_order > target_order ~ "F",
TRUE ~ NA)) %>%
ggplot(aes(x = factor(seed_structure, levels = strucs), y = factor(target_structure, levels = rev(strucs)))) +
geom_tile(aes(fill = plot_value, color = plot_sex), linewidth = 8) +
geom_tile(aes(fill = plot_value), linewidth = 8) +
scale_fill_continuous_diverging(palette = "Cork", name = "Fisher-\ntransformed\ncorrelation",
limits = c(-0.2, 0.2)) +
scale_color_manual(values = c("M" = "blue", "F" = "red"), name = "Sex", na.translate = FALSE) +
scale_x_discrete(labels = new_labels) + # Apply new labels to x-axis
scale_y_discrete(labels = new_labels) + # Apply new labels to y-axis
xlab("Seed region") +
ylab("Target region") +
ggtitle("Seed Correlation by Sex") +
coord_fixed(ratio = 1) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "right",
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
axis.title = element_text(size = 12),
axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels by 45 degrees
axis.text = element_text(size = 10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
# 'palette' should be one of: Blue-Red, Blue-Red 2, Blue-Red 3, Red-Green, Purple-Green, Purple-Brown, Green-Brown, Blue-Yellow 2, Blue-Yellow 3,
# Green-Orange, Cyan-Magenta, Tropic, Broc, Cork, Vik, Berlin, Lisbon, Tofino
######################################
df_ttest <- df %>%
mutate(z_fixed=case_when(seed_structure==target_structure ~ 0,
TRUE ~ z)) %>%
group_split(seed_structure, target_structure) %>%
map_dfr(function(x) {
x %>%
t.test(z_fixed ~ sex, data = .) %>%
broom::tidy() %>%
rename(diff=estimate, mean_F=estimate1, mean_M=estimate2) %>%
mutate(seed_structure=x$seed_structure[[1]], target_structure=x$target_structure[[1]])
}) %>%
mutate(seed_order=as.integer(str_replace(seed_structure, "X", "")),
target_order=as.integer(str_replace(target_structure, "X", "")))
df_ttest_labels <- df_ttest %>%
mutate(seed_order=as.integer(str_replace(seed_structure, "X", "")),
target_order=as.integer(str_replace(target_structure, "X", ""))) %>%
filter(target_order > seed_order) %>%
mutate(significant=(p.value < 0.05)) %>%
mutate(sigtext=case_when(p.value < 0.001 ~ "***",
p.value < 0.01 ~ "**",
p.value < 0.05 ~ "*",
p.value < 0.1 ~ "-",
TRUE ~ "")) %>%
mutate(seed_structure=factor(seed_structure, levels=strucs),
target_structure=factor(target_structure, levels=rev(strucs)))
# Check if the labels appear correctly
print(df_ttest_labels)
# Apply Benjamini-Hochberg correction
df_ttest <- df_ttest %>%
mutate(p.adjusted = p.adjust(p.value, method = "BH"))
df_ttest_labels <- df_ttest %>%
mutate(seed_order = as.integer(str_replace(seed_structure, "X", "")),
target_order = as.integer(str_replace(target_structure, "X", ""))) %>%
filter(target_order > seed_order) %>%
mutate(significant = (p.adjusted < 0.05)) %>%
mutate(sigtext = case_when(p.adjusted < 0.001 ~ "***",
p.adjusted < 0.01 ~ "**",
p.adjusted < 0.05 ~ "*",
p.adjusted < 0.1 ~ "-",
TRUE ~ "")) %>%
mutate(seed_structure = factor(seed_structure, levels = strucs),
target_structure = factor(target_structure, levels = rev(strucs)))
# Create df_ttest_labels with Benjamini-Hochberg correction and new labels
df_ttest_labels <- df_ttest %>%
mutate(seed_order = as.integer(str_replace(seed_structure, "X", "")),
target_order = as.integer(str_replace(target_structure, "X", ""))) %>%
filter(target_order > seed_order) %>%
mutate(significant = (p.adjusted < 0.05)) %>%
mutate(sigtext = case_when(p.adjusted < 0.001 ~ "***",
p.adjusted < 0.01 ~ "**",
p.adjusted < 0.05 ~ "*",
p.adjusted < 0.1 ~ "-",
TRUE ~ "")) %>%
mutate(seed_structure = factor(seed_structure, levels = strucs),
target_structure = factor(target_structure, levels = rev(strucs))) %>%
mutate(seed_structure = recode(seed_structure, !!!new_labels),
target_structure = recode(target_structure, !!!new_labels))
library(RColorBrewer)
# Calculate the maximum absolute t-value from the 'statistic' column
max_abs_t <- max(abs(df_ttest$statistic), na.rm = TRUE)
#DR
new_labels <- c("X1" = "SM", "X2" = "SMS", "X3" = "THAL", "X4" = "STR")
#new_labels=c("X1"="dSTR","X2"="FAC","X3"="ICC","X4"="INSU","X5"="MC1",
# "X6"="MC2","X7"="RPC","X8"="SMS","X9"="SMS2","X10"="THAL",
# "X11"="vHIP","X12"="NAC")
df_ttest %>%
ggplot(aes(x = factor(seed_structure, levels = strucs), y = factor(target_structure, levels = rev(strucs)))) +
geom_tile(aes(fill = statistic), color = 'white', linewidth = 0.2) +
scale_fill_gradient2(
low = '#0000FF', mid = '#FFFFFF', high = '#FF0000',
midpoint = 0, name = "t-statistic",
limits = c(-max_abs_t, max_abs_t),
breaks = c(-max_abs_t, 0, max_abs_t),
labels = c(sprintf("%.1f", -max_abs_t), "0", sprintf("%.1f", max_abs_t))
) +
geom_text(aes(label = sigtext), data = df_ttest_labels, size = 6, color = "green") +
xlab("Seed region") +
ylab("Target region") +
coord_fixed(ratio = 1) +
ggtitle("Sex Differences in Network Correlation") +
scale_x_discrete(labels = new_labels) + # Apply new labels to x-axis
scale_y_discrete(labels = new_labels) + # Apply new labels to y-axis
theme_bw() +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "right",
legend.title = element_text(size = 10),
legend.text = element_text(size = 8),
axis.title = element_text(size = 12),
axis.text.x = element_text(angle = 45, hjust = 1), # Rotate x-axis labels by 45 degrees
axis.text = element_text(size = 10),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
library(gridExtra)
library(png)
library(grid)
# Define the output file path
output_file <- "/data/scratch2/marina/FINALexp/functional_correlations/SBC_QC/stats_table.png"
# Create the table as a graphical object
table_plot <- tableGrob(df_ttest_labels)
# Open a PNG device
png(filename = output_file, width = 12000, height = 20000, res = 600) # Adjust width, height, and resolution as needed
# Draw the table on the PNG device
grid.draw(table_plot)
# Close the device to save the file
dev.off()
# Confirmation message
cat("Normality test statistics table saved as PNG at:", output_file, "\n")