09_brain_functions_proteins.Rmd

---
title: "brain_funtion_protein"
author: "Yuting Xie"
date: "2023-04-26"
output: html_document
---

```{r}
rm(list=ls())
pacman::p_load(readxl,magrittr,
               RColorBrewer,
               Rtsne, umap,
               pheatmap,vioplot,
               ggpubr, ggplot2,
               corrplot, stringr,preprocessCore,reshape2,dplyr)
```

#expression median/mean/sum(intensity)
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
FP <- read.csv('data_function_protein_0.6_20230329.csv')
FP$pro_iso <- gsub("[.]", "-", FP$pro_iso)
functions <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
FP <- FP[-unlist(sapply(functions, function(x){(which(FP$func==x))})), ]

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
info<-data.table::fread('FIleInfo_HBA_20220816_withoutFALSE_and_pool.csv', stringsAsFactors = F, data.table = F, fill = T)

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
df_nor<-data.table::fread('HBA_prot_matrix_normolized_20220811.csv',data.table = F)
rownames(df_nor)<-df_nor[,1]
df_nor<-df_nor[,-1]
df_nor1 <- df_nor[-grep("pool", rownames(df_nor)),]

df_nor1$matter<-info$W_G[match(rownames(df_nor1),info$sample)]
unique(df_nor1$matter) #"G" "W" NA
df_nor2<-df_nor1[grep("G",df_nor1$matter),]
df_nor3<-df_nor2[!is.na(df_nor2$matter),]
df_nor3<-df_nor3[,-ncol(df_nor3)]
df_nor3$pathology <- info$pathology[match(rownames(df_nor3),info$sample)]
df_nor4 <- df_nor3[which(df_nor3$pathology=="health"), ]

df_nor5 <- df_nor4[,match(unique(FP$pro_iso),colnames(df_nor4))]
# dfMedian<-apply(df_nor5, 2, function(x){median(x)})
# dfMean <- apply(df_nor5, 2, function(x){mean(x)})
# dfMedianMean <- data.frame(pro_iso=colnames(df_nor5), median=dfMedian, mean=dfMean)
df_nor6 <- 2^df_nor5 
dfSUM <- apply(df_nor5, 2, function(x){sum(x)})
dfSUM <- data.frame(pro_iso=colnames(df_nor6), sum=dfSUM, log2SUM=log2(dfSUM))

FP1 <- FP[,c(3,2)]
datF <- merge(FP1, dfSUM, by="pro_iso")
# df_label <- datF[order(datF$median, decreasing = TRUE),]
# df_label <- df_label[1:25, ]

ggplot(datF,aes(x=log2SUM, y=func))+
  geom_jitter(size = 0.5, width = 0.05, show.legend = FALSE, color="#05635C")+ #, color="#10798B"
  theme_bw()
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
ggsave("brain_function_sum_intensity_expression.pdf", width=6, height = 6)
```

#number of process
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
FP <- read.csv("data_function_protein_0.6_20230607.csv")
colnames(FP)[2] <- "pro_iso" 

FP$pro_iso <- gsub("[.]", "-", FP$pro_iso)
FP1 <- unique(FP[,c(2,3,4)])
# delet functions
functions <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
FP2 <- FP1    #[-unlist(sapply(functions, function(x){(which(FP1$func==x))})), ]

PFN <- as.data.frame(table(FP2$pro))
colnames(PFN) <- c("pro", "number")

datF <- merge(FP2, PFN, by="pro")

# df_label <- datF[order(as.numeric(datF$number), decreasing = FALSE),]
# df_label <- df_label[1:30, ]

library(ggplot2)
library(tidyverse)
p1<-ggplot(data=datF,aes(x=number, y=func))+
    geom_boxplot()

datF %>% 
  group_by(func) %>% 
  summarise(median_value=median(number)) %>%
  rename("func_1"="func") %>% 
  cbind(ggplot_build(p1)$data[[1]]) -> datF1

ggplot(datF,aes(x=number, y=func))+
  stat_boxplot(geom="errorbar",width=0.3,linewidth=0.5,position=position_dodge(0.8),color="#ADBAC2")+
  geom_boxplot(width=0.75, size=0.3,
                      # aes(colour = pathology),
                      # fill="black",
                      # alpha=0.8,
                      color="#ADBAC2",
                      fill="#ADBAC2",
                      outlier.fill = "#ADBAC2",
                      outlier.colour="#ADBAC2",  #"black",
                      outlier.shape = 19,
                      outlier.size = 0.01,
                      position = position_dodge(0.8))+
  geom_segment(data=datF1,size=0.3,
               aes(y=ymin,yend=ymax,
                   x=median_value,
                   xend=median_value),
               color="red")+
  theme_bw()+
  xlim(37, 0)
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
ggsave("process_number_boxplot_20230619.pdf", width=5, height = 6)

```

#number of process isoform and no_isoform
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
FP <- read.csv("data_function_protein_0.6_20230607.csv")
colnames(FP)[2] <- "pro_iso" 

FP$pro_iso <- gsub("[.]", "-", FP$pro_iso)
FP1 <- unique(FP[,c(2,3,4)])
# delet functions
# functions <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
FP2 <- FP1    #[-unlist(sapply(functions, function(x){(which(FP1$func==x))})), ]

iso <- FP2[grep("-", FP2$pro_iso), ]
iso1 <- as.data.frame(table(iso$pro_iso))
colnames(iso1) <- c("pro_iso", "number")
datFiso <- merge(iso, iso1, by="pro_iso")

noneiso <- FP2[-grep("-", FP2$pro_iso), ]
noneiso1 <- as.data.frame(table(noneiso$pro_iso))
colnames(noneiso1) <- c("pro_iso", "number")
datFnoneiso <- merge(noneiso, noneiso1, by="pro_iso")

#iso
library(ggplot2)
library(tidyverse)
p1<-ggplot(data=datFiso,aes(x=number, y=func))+
    geom_boxplot()

datFiso %>% 
  group_by(func) %>% 
  summarise(median_value=median(number)) %>%
  rename("func_1"="func") %>% 
  cbind(ggplot_build(p1)$data[[1]]) -> datFiso1
ggplot(datFiso,aes(x=number, y=func))+
  stat_boxplot(geom="errorbar",width=0.3,linewidth=0.5,position=position_dodge(0.8),color="#ADBAC2")+
  geom_boxplot(width=0.75, size=0.3,
                      # aes(colour = pathology),
                      # fill="black",
                      # alpha=0.8,
                      color="#ADBAC2",
                      fill="#ADBAC2",
                      outlier.fill = "#ADBAC2",
                      outlier.colour="#ADBAC2",  #"black",
                      outlier.shape = 19,
                      outlier.size = 0.01,
                      position = position_dodge(0.8))+
  geom_segment(data=datFiso1,size=0.3,
               aes(y=ymin,yend=ymax,
                   x=median_value,
                   xend=median_value),
               color="red")+
  theme_bw()+
  xlim(37, 0)
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
ggsave("process_number_boxplot_isoform_20230818.pdf", width=5, height = 6)

#none iso
library(ggplot2)
library(tidyverse)
p1<-ggplot(data=datFnoneiso,aes(x=number, y=func))+
    geom_boxplot()

datFnoneiso %>% 
  group_by(func) %>% 
  summarise(median_value=median(number)) %>%
  rename("func_1"="func") %>% 
  cbind(ggplot_build(p1)$data[[1]]) -> datFnoneiso1
ggplot(datFnoneiso,aes(x=number, y=func))+
  stat_boxplot(geom="errorbar",width=0.3,linewidth=0.5,position=position_dodge(0.8),color="#ADBAC2")+
  geom_boxplot(width=0.75, size=0.3,
                      # aes(colour = pathology),
                      # fill="black",
                      # alpha=0.8,
                      color="#ADBAC2",
                      fill="#ADBAC2",
                      outlier.fill = "#ADBAC2",
                      outlier.colour="#ADBAC2",  #"black",
                      outlier.shape = 19,
                      outlier.size = 0.01,
                      position = position_dodge(0.8))+
  geom_segment(data=datFnoneiso1,size=0.3,
               aes(y=ymin,yend=ymax,
                   x=median_value,
                   xend=median_value),
               color="red")+
  theme_bw()+
  xlim(37, 0)
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
ggsave("process_number_boxplot_none_isoform_20230818.pdf", width=5, height = 6)

datFiso$group <- "isoform" 
datFnoneiso$group <- "noneisoform"
datall <- rbind(datFiso, datFnoneiso)

ggplot(datall, aes(x=number, y=func, fill=group))+
  scale_fill_manual(values = c("#EE7785","#67D5B5"))+
  scale_colour_manual(values = c("white","white"))+  #aesthetics = "colour"
  stat_boxplot(geom="errorbar",width=0.3,size=0.3,position=position_dodge(0.8),color="black")+
  geom_boxplot(width=0.3, size=0.3,
                      aes(colour = pathology),
                      # fill="black",
                      alpha=0.8,
                      color="black",
                      outlier.fill = "black",
                      outlier.colour="black",  #"black",
                      outlier.shape = 19,
                      outlier.size = 0.01,
                      position = position_dodge(0.8))+
  # geom_segment(data=datFnoneiso1,size=0.3,
  #              aes(y=ymin,yend=ymax,
  #                  x=median_value,
  #                  xend=median_value),
  #              color="red")+
  theme_bw()
  # xlim(37, 0)
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230426_brain_function")
ggsave("process_number_boxplot_isoform_with_none_isoforms_20230818.pdf", width=5, height = 8)




```

# process_proteins heatmap (RE ratio)
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
ts <- read.csv("HBA_ts_9639pro_1521sam_Grey_BBA_3_2_2.5_20220921.csv", row.names = 1)
FuncPro <- read.csv("data_function_protein_0.6_20230329.csv")
FuncBBA <- read.csv("BDf_FDR05.csv")

FuncPro1 <- FuncPro[,c(3,2)]
FuncPro1$pro_iso <- gsub("[.]", "-", FuncPro1$pro_iso)
# delet functions
Dfunc <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
FuncPro2 <- FuncPro1[-unlist(sapply(Dfunc, function(x){(which(FuncPro1$func==x))})), ]

FuncBBA$BD <- gsub("[.]", "", FuncBBA$BD)

tsBBAname <- data.frame(names=rownames(ts))
tsBBAname$BBA <- sapply(tsBBAname$names, function(x){str_split(x, "_G")[[1]][1]})
tsBBAname$BBA <- sapply(tsBBAname$BBA, function(x){str_split(x, ",")[[1]][1]})
tsBBAname$BBA[which(tsBBAname$BBA=="Hippocampus_Hipp")] <- "Hipp"
tsBBAname$BBA[which(tsBBAname$BBA=="Amygdala_Amyg")] <- "Amyg"

FuncBBA$BBA[c(which(FuncBBA$BBA=="rHipp"), which(FuncBBA$BBA=="cHipp"))] <- "Hipp"
FuncBBA$BBA[c(which(FuncBBA$BBA=="lAmyg"), which(FuncBBA$BBA=="mAmyg"))] <- "Amyg"

tsBBAname$match <- FuncBBA$BBA[match(tsBBAname$BBA, FuncBBA$BBA)]
matchBBA <- tsBBAname$match[!is.na(tsBBAname$match)]

FuncBBA2 <- FuncBBA[unlist(sapply(matchBBA, function(x){which(FuncBBA$BBA==x)})), ]
ts2 <- ts[!is.na(tsBBAname$match), ]
rownames(ts2) <- matchBBA
colnames(ts2) <- gsub("_", "-", colnames(ts2))

Func <- unique(FuncPro2$func)
datHeat <- data.frame(proteins=1)
for(i in 1:length(Func)){
  d1 <- FuncPro2[which(FuncPro2$func==Func[i]), ]
  d2 <- FuncBBA2[which(FuncBBA2$BD==Func[i]), ]
  pro <- unique(d1$pro_iso)
  BBA <- unique(d2$BBA)
  
  d3 <- ts2[BBA, pro]
  d4 <- apply(d3, 2, function(x){length(which(x==4))/nrow(d3)})
  datHeat0 <- data.frame(proteins=colnames(d3), percent=d4)
  colnames(datHeat0)[2] <- Func[i]
  datHeat <- merge(datHeat, datHeat0, by="proteins", all=T)
}
datHeat1 <- datHeat[-1, ]
rownames(datHeat1) <- datHeat1[,1]
datHeat1 <- datHeat1[,-1]
datHeat2 <- as.data.frame(t(datHeat1))
datHeat2 <- datHeat2[order(rownames(datHeat2)), ]
datHeat2[is.na(datHeat2)] <- 0

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
library(pheatmap)
pheatmap(datHeat2,
         fontsize=15,
         color = colorRampPalette(c("#EDF9F3","#495FA3"))(10), 
         #annotation_col = colnames(pathway2),
         #annotation_colors = list(W_G=c(rep("#4472C4",26),rep("#EE781C",19),Region=rainbow(45))),
         # na_color = "white",
         border_color="#BCBCBC",
         trace = "none", 
         scale = "none", 
         cluster_rows = T,
         cluster_cols = T,
         show_rownames = T, 
         show_colnames = T, 
         fontsize_col=1,
         fontsize_row=6,
         filename = "Function_proteins_RE_Heatmap_2.pdf",
         # cellwidth=5,
         # cellheight = 5,
         # main = "Number of differentially expressed RNA between pairs of brain regions in grey matter_tsData",
         height = 5,width =8
         )

```

# process proteins heatmap (proteins' intensity in BBA of functions)
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
ts <- read.csv("HBA_ts_9639pro_1521sam_Grey_BBA_3_2_2.5_20220921.csv", row.names = 1)
ts$bba <- rownames(ts)
FuncPro <- read.csv("data_function_protein_0.6_20230607.csv")
FuncBBA <- read.csv("BDf_FDR05.csv")

FuncPro1 <- FuncPro[,c(3,2)]
colnames(FuncPro1)[2] <- "pro_iso" 
FuncPro1$pro_iso <- gsub("[.]", "-", FuncPro1$pro_iso)
# # delet functions
# Dfunc <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
# FuncPro2 <- FuncPro1[-unlist(sapply(Dfunc, function(x){(which(FuncPro1$func==x))})), ]
FuncPro2 <- FuncPro1

FuncBBA$BD <- gsub("[.]", "", FuncBBA$BD)

tsBBAname <- data.frame(names=rownames(ts))
tsBBAname$BBA <- sapply(tsBBAname$names, function(x){str_split(x, "_G")[[1]][1]})
tsBBAname$BBA <- sapply(tsBBAname$BBA, function(x){str_split(x, ",")[[1]][1]})
tsBBAname$BBA[which(tsBBAname$BBA=="Hippocampus_Hipp")] <- "Hipp"
tsBBAname$BBA[which(tsBBAname$BBA=="Amygdala_Amyg")] <- "Amyg"
tsBBAname$BBA[which(tsBBAname$BBA=="TE1.0 and TE1.2")] <- "TE1.0/TE1.2"

FuncBBA$BBA[c(which(FuncBBA$BBA=="rHipp"), which(FuncBBA$BBA=="cHipp"))] <- "Hipp"
FuncBBA$BBA[c(which(FuncBBA$BBA=="lAmyg"), which(FuncBBA$BBA=="mAmyg"))] <- "Amyg"

tsBBAname$match <- FuncBBA$BBA[match(tsBBAname$BBA, FuncBBA$BBA)]
matchBBA <- tsBBAname$match[!is.na(tsBBAname$match)]
tsBBAname$BBAname <- sapply(tsBBAname$names, function(x){str_split(x, "_G")[[1]][1]})

FuncBBA2 <- FuncBBA[unlist(sapply(matchBBA, function(x){which(FuncBBA$BBA==x)})), ]
ts2 <- ts[!is.na(tsBBAname$match), ]
rownames(ts2) <- matchBBA
colnames(ts2) <- gsub("_", "-", colnames(ts2))

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
info<-data.table::fread('FIleInfo_HBA_20220816_withoutFALSE_and_pool.csv', stringsAsFactors = F, data.table = F, fill = T)
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
df_nor<-data.table::fread('HBA_prot_matrix_normolized_20220811.csv',data.table = F)
rownames(df_nor)<-df_nor[,1]
df_nor<-df_nor[,-1]
df_nor1 <- df_nor[-grep("pool", rownames(df_nor)),]
df_nor1$matter<-info$W_G[match(rownames(df_nor1),info$sample)]
unique(df_nor1$matter) #"G" "W" NA
df_nor2<-df_nor1[grep("G",df_nor1$matter),]
df_nor3<-df_nor2[!is.na(df_nor2$matter),]
df_nor3<-df_nor3[,-ncol(df_nor3)]
df_nor3$pathology <- info$pathology[match(rownames(df_nor3),info$sample)]
df_nor4 <- df_nor3[which(df_nor3$pathology=="health"), ]
df_nor4 <- df_nor4[, -ncol(df_nor4)]
df_nor4$BBA <- info$BBA_Anatomical.and.modified.Cyto.architectonic.descriptions[match(rownames(df_nor4), info$sample)]
df_nor4$BBAmatch <- tsBBAname$match[match(df_nor4$BBA, tsBBAname$BBAname)]

Func <- unique(FuncPro2$func)
datHeat <- data.frame(proteins=1)
for(i in 1:length(Func)){
  # d1 <- FuncPro2[which(FuncPro2$func==Func[i]), ]
  d2 <- FuncBBA2[which(FuncBBA2$BD==Func[i]), ]
  pro <- unique(FuncPro2$pro_iso)
  BBA <- unique(d2$BBA)
  
  d3 <- df_nor4[unlist(sapply(BBA, function(x){which(df_nor4$BBAmatch==x)})), pro]
  # d3 <- df_nor4[unlist(sapply(BBA, function(x){which(df_nor4$BBAmatch==x)})), c(pro, "BBAmatch")]
  
  #median
  d4 <- apply(d3, 2, function(x){mean(x)})  #median(x)
  # d4 <- melt(d3)
  # d4.1 <- dcast(d4, BBAmatch~variable, median)
  # rownames(d4.1) <- d4.1[,1]
  # d4 <- d4.1[,-1] 
  
  datHeat0 <- data.frame(proteins=colnames(d3), m=d4)
  colnames(datHeat0)[2] <- Func[i]

  # datHeat0 <- as.data.frame(t(d4))
  # colnames(datHeat0) <- paste0(Func[i], "_", colnames(datHeat0))
  # datHeat0$proteins <- rownames(datHeat0)

  datHeat <- merge(datHeat, datHeat0, by="proteins", all=T)
  
#   setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
# library(pheatmap)
# pheatmap(datHeat0,
#          fontsize=15,
#          color = colorRampPalette(c("#EDF9F3","#495FA3"))(10), 
#          #annotation_col = colnames(pathway2),
#          #annotation_colors = list(W_G=c(rep("#4472C4",26),rep("#EE781C",19),Region=rainbow(45))),
#          # na_color = "white",
#          border_color="#BCBCBC",
#          trace = "none", 
#          scale = "none", 
#          cluster_rows = T,
#          cluster_cols = T,
#          show_rownames = T, 
#          show_colnames = T, 
#          fontsize_col=4,
#          fontsize_row=4,
#          # filename = paste0(Func[i], "_proteins_intensity_in_BBA_heatmap.pdf"),
#          # cellwidth=5,
#          # cellheight = 5,
#          # main = "Number of differentially expressed RNA between pairs of brain regions in grey matter_tsData",
#          height = 10,width = 10
#          )
 }
datHeat1 <- datHeat[-1, ]
rownames(datHeat1) <- datHeat1[,1]
datHeat1 <- datHeat1[,-1]
datHeat2 <- as.data.frame(t(datHeat1))
datHeat2 <- datHeat2[order(rownames(datHeat2)), ]
datHeat2[is.na(datHeat2)] <- 0
#scale proteins
datHeat2 <- scale(datHeat2)
# datHeat2$FUN <- rownames(datHeat2)

# FuncPro2$value <- c(1)
# FuncPro3 <- dcast(FuncPro2, func~pro_iso)
# rownames(FuncPro3) <- FuncPro3[,1]
# FuncPro3 <- FuncPro3[,-1]
# annotation_col <- data.frame(anno=apply(FuncPro3, 2, function(x){sum(!is.na(x))}))
# annotation_col$anno <- as.character(annotation_col$anno)

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
protLocation <- read.csv("173_proteins_location.csv", row.names = 1)
annotation_col <- protLocation
annotation_row <- read.csv("brain_functions_ABB.csv", row.names = 1)
# annotation_colors <- read.csv("function_location_colors.csv", row.names = 1)
annotation_colors <- list(Subcellular.location=c(Others="#E96B65", Cytoplasm="#F4AFAD", Secreted="#FEF3F0", `Cell membrane`="#1C919D", Mitochondrion="#0E5E69", Nucleus="#F4BE4A"), Function=c(Action="#EAD6B2", Cognition="#D57560", Emotion="#3E3E5A", Interoception="#7FAC9B", Perception="#EAC585"))

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
# write.csv(pro, "functions_proteins_330.csv")
library(pheatmap)
pheatmap(datHeat2,
         fontsize=15,
         color = colorRampPalette(c("#FFFF33","#F4C70B","#AD4462","#67277C"))(10),  #"#67277C", "#AD4462", "#F4C70B", "#FFFF33"
         annotation_col = annotation_col,
         # annotation_colors = list(W_G=c(rep("#4472C4",26),rep("#EE781C",19),Region=rainbow(45))),
         annotation_row = annotation_row,
         annotation_colors = annotation_colors,
         # na_color = "white",
         border_color=NA,    #"#EFEFEF",
         trace = "none", 
         scale = "none", 
         cluster_rows = T,
         cluster_cols = T,
         show_rownames = T, 
         show_colnames = T, 
         fontsize_col=2,
         fontsize_row=4,
         filename = "Heatmap_log2_proteins_intensity_mean_function_BBA_20230618_scale.pdf",
         # cellwidth=4,
         # cellheight = 6,
         # main = "",
         height = 7,width = 15
         )

#UAMP
# datUmap <- as.data.frame(t(datHeat2))
datUmap <- as.data.frame(datHeat2)
set.seed(1)
datUmap1 <- uwot::umap(datUmap, n_neighbors = 5, #2~100
                 metric = "cosine", spread = 1,
                 min_dist = 0.01, verbose = TRUE)  #t(datExprT) as.dist(datExprT)  scale="z"
  #scale:"Z" or "scale"    metric="cosine" correlation euclidean
datUmap2<-as.data.frame(datUmap1)
colnames(datUmap2) <- c("UMAP1","UMAP2")
datUmap2$proteins <- rownames(datUmap2)
Function <- c("Action", "Cognition", "Emotion", "Interoception", "Perception")

datUmap3 <- NULL
for(i in 1:length(Function)){
  d1 <- datUmap2[grep(Function[i], datUmap2$proteins),]
  d1$group <- Function[i]
  datUmap3 <- rbind(datUmap3, d1)
}

ggplot(datUmap2,aes(UMAP1,UMAP2,color=datUmap3$group))+  #OrganABB  #,color=datUmap3$group
          geom_point(size=2)+ theme_classic() +
          theme(panel.border = element_rect(fill=NA,color="black", size=0.3, linetype="solid"))+
          #xlim(-3,4)+ylim(-4,4)+ 
          # geom_hline(yintercept = 0,lty=2,col="red",lwd=0.2) + 
          # geom_vline(xintercept = 0,lty=2,col="blue",lwd=0.2) +
          theme(plot.title = element_text(hjust = 0.5)) + 
          geom_text(aes(label=proteins,vjust=1,hjust=-0.5,angle=15, color="black"),size=0.5,show.legend=T)+
          labs(title = "Umap")+
          scale_color_manual(values=c(Action="#EAD6B2", Cognition="#D57560", Emotion="#3E3E5A", Interoception="#7FAC9B", Perception="#EAC585"))
          # scale_fill_manual(values=c("#ff9a9e","#a3bded"))
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
ggsave("function_proteins_umap_20230618_with_label_function_cluster.pdf",
       height = 6,
       width = 8)

#tsne
set.seed(1)
df.PC <- Rtsne::Rtsne(datUmap,
                 dims = 2,
                 perplexity = 10,  #0-50
                 verbose = T,
                 theta = 0.1,
                 # pca=TRUE,
                 partial_pca=FALSE,
                 pca_scale=FALSE,
                 )
df.PC.1 <- as.data.frame(df.PC$Y)
rownames(df.PC.1)<-rownames(datUmap)
colnames(df.PC.1) <- c("PC1","PC2")
df.PC.1$proteins <- rownames(df.PC.1)

ggplot(df.PC.1,aes(PC1,PC2))+  #OrganABB
          geom_point(size=2)+ theme_classic() +
          theme(panel.border = element_rect(fill=NA,color="black", size=0.5, linetype="solid"))+
          #xlim(-3,4)+ylim(-4,4)+ 
          # geom_hline(yintercept = 0,lty=2,col="red",lwd=0.2) + 
          # geom_vline(xintercept = 0,lty=2,col="blue",lwd=0.2) +
          theme(plot.title = element_text(hjust = 0.5)) + 
          geom_text(aes(label=proteins,vjust=1,hjust=-0.8,angle=0),size=0.5,show.legend=T)+
          labs(title = "TSNE")
          # scale_color_manual(values = c(colors$colors))
          # geom_text_repel(label=paste(df.PC.1$OrganABB),size=1)
          # scale_color_hue(h=c(0,1000),c=50)
          #scale_fill_brewer(palette = 'YlOrRd')
          #scale_color_brewer(palette="Paired")

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
ggsave("function_proteins_TSNE.pdf",
       height = 6,
       width = 6) 
```

#functions_regions heatmap 
```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
datfunction <- read.csv("BDf_FDR05.csv")
datfunction$BD <- paste0("AAA", datfunction$BD)
datfunction$BD <- gsub("AAA.", "", datfunction$BD)

#不分左右
datF <- unique(datfunction[,c("BBA", "BD", "score")])
datF1 <-  dcast(datF, BD ~ BBA, mean)
rownames(datF1) <- datF1$BD
datF1 <- datF1[,-1]

datOrder <- unique(datfunction[,c(1:3)])
datOrder$LR <- sapply(datOrder$BBA_name, function(x){str_split(x, "_")[[1]][2]})
datOrder$BBA_LR <- paste0(datOrder$BBA,"_",datOrder$LR)

datF1 <- datF1[, unique(datOrder$BBA)]
datF1 <- log2(datF1)
# datF1[!is.na(datF1)] <- 2
# datF1[is.na(datF1)] <- 0

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
ts <- read.csv("HBA_ts_9639pro_1521sam_Grey_BBA_3_2_2.5_20220921.csv", row.names = 1)
tsBBAname <- data.frame(names=rownames(ts))
tsBBAname$BBA <- sapply(tsBBAname$names, function(x){str_split(x, "_G")[[1]][1]})
tsBBAname$BBA <- sapply(tsBBAname$BBA, function(x){str_split(x, ",")[[1]][1]})
tsBBAname$BBA[which(tsBBAname$BBA=="Hippocampus_Hipp")] <- "Hipp"
tsBBAname$BBA[which(tsBBAname$BBA=="Amygdala_Amyg")] <- "Amyg"
tsBBAname$BBA[which(tsBBAname$BBA=="TE1.0 and TE1.2")] <- "TE1.0/TE1.2"

colnames(datF1)[c(which(colnames(datF1)=="rHipp"), which(colnames(datF1)=="cHipp"))] <- "Hipp"
colnames(datF1)[c(which(colnames(datF1)=="lAmyg"), which(colnames(datF1)=="mAmyg"))] <- "Amyg"

com_BBA <- intersect(colnames(datF1), tsBBAname$BBA)

datF2 <- datF1[,com_BBA]

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
FuncPro <- read.csv("data_function_protein_0.6_20230607.csv")
rownames(datF2) <- gsub("[.]", "", rownames(datF2))
com_func <- intersect(unique(FuncPro$func), rownames(datF2))
datF3 <- datF2[com_func,]
datF3 <- datF3[order(rownames(datF3)), ]

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
pdf("BBAregion_function_heatmap.pdf", width=20, height=10)
pheatmap::pheatmap(datF3,
             color = colorRampPalette(c("#07725D", "#07725D"))(1),  
             border_color = "#EDEDED",
             na_col="white",
             # gaps_col = num,
             #legend_breaks=seq(0,246,by=6),
             scale = "none",  #row
             cluster_rows = FALSE,
             cluster_cols = FALSE,
             show_colnames = T,
             show_rownames = T,
             cellwidth = 10,
             cellheight = 10,
             display_numbers = F,
             fontsize_number = 3,
             fontsize_col=5,
             fontsize_row=5,
             # legend_breaks = seq(0.2, 4.2, by=0.2),
             breaks=NA
             # legend_labels = NA
             # annotation_col = anno,
             # annotation_colors = anno_colors
            )
dev.off()


#分左右
datF <- unique(datfunction[,c("BBA_name", "BD", "score")])
datF1 <-  dcast(datF, BD ~ BBA_name)
datF1 <- data.frame(datF1, row.names="BD")

datOrder <- unique(datfunction[,c(1:3)])
datOrder$LR <- sapply(datOrder$BBA_name, function(x){str_split(x, "_")[[1]][2]})
datOrder$BBA_LR <- paste0(datOrder$BBA,"_",datOrder$LR)

datF1 <- datF1[, datOrder$BBA_name]
colnames(datF1) <- datOrder$BBA_LR
datF1 <- log2(datF1)

gaps_col <- as.data.frame(table(datOrder$BBA))
rownames(gaps_col) <- gaps_col$Var1
gaps_col <- gaps_col[unique(datOrder$BBA),]
num <- c()
for(i in 1:nrow(gaps_col)){
  num0 <- sum(gaps_col$Freq[1:i])
  num <-c(num,num0)
}

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/figures/add")
pdf("brain_function_heatmap_L_R.pdf", width=20, height=10)
pheatmap::pheatmap(datF1,
             color = colorRampPalette(c("white", "#9A96B5"))(20),  
             border_color = NA,
             na_col="#EDEDED", 
             gaps_col = num,
             #legend_breaks=seq(0,246,by=6),
             scale = "none",  #row
             cluster_rows = F,
             cluster_cols = F,
             show_colnames = T,
             show_rownames = T,
             cellwidth = 4, 
             cellheight = 5,
             display_numbers = F,
             fontsize_number = 3,
             fontsize_col=3,
             fontsize_row=5
             # annotation_col = anno,
             # annotation_colors = anno_colors
            )
dev.off()





datfunction1 <- datfunction[,c(2,4,5)]
library(reshape2)
datfunction2 <- dcast(datfunction1, BD ~ BBA, mean)

# setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
# info<-data.table::fread('FIleInfo_HBA_20220816_withoutFALSE_and_pool.csv', stringsAsFactors = F, data.table = F, fill = T)
# info1 <- unique(info[,c(16,19)])
# datOrder <- unique(datfunction[,1:3])
# datOrder$lobe <- info1$Lobe[unlist(sapply(datOrder$BBA_name, function(x){grep(x, info1$BBA_name)}))]





```

#20230713 proteins expression (mean) in each person and in BBA of brain functions
```{r}
#the proteins of function
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
FuncPro <- read.csv("data_function_protein_0.6_20230607.csv")
FuncBBA <- read.csv("BDf_FDR05.csv")

FuncPro1 <- FuncPro[,c(3,2)]
colnames(FuncPro1)[2] <- "pro_iso" 
FuncPro1$pro_iso <- gsub("[.]", "-", FuncPro1$pro_iso)
# delet functions
# Dfunc <- c("PerceptionVisionColor", "Interoception", "EmotionAnxiety", "Emotion", "Cognition")
# FuncPro2 <- FuncPro1[-unlist(sapply(Dfunc, function(x){(which(FuncPro1$func==x))})), ]
# FuncPro2 <- FuncPro1

#the BBAs of function
FuncBBA$BD <- gsub("[.]", "", FuncBBA$BD)
FuncBBA$BBA[c(which(FuncBBA$BBA=="rHipp"), which(FuncBBA$BBA=="cHipp"))] <- "Hippocampus_Hipp"
FuncBBA$BBA[c(which(FuncBBA$BBA=="lAmyg"), which(FuncBBA$BBA=="mAmyg"))] <- "Amygdala_Amyg"
FuncBBA$BBA[c(which(FuncBBA$BBA=="TE1.0/TE1.2"))] <- "TE1.0 and TE1.2"
func <- unique(FuncPro1$func) 
FuncBBA2 <- FuncBBA[unlist(sapply(func, function(x){(which(FuncBBA$BD==x))})), ]
FuncBBA3 <- FuncBBA2[,c(2,4)]
FuncBBA4 <- dcast(FuncBBA3, BBA~BD)

#the data of proteins expression
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
info<-read.csv('FIleInfo_HBA_20220816_withoutFALSE_and_pool.csv')
info$BBA_ABB <- sapply(info$BBA_Anatomical.and.modified.Cyto.architectonic.descriptions, function(x){str_split(x, ",")[[1]][1]})

df_nor<-data.table::fread('HBA_prot_matrix_normolized_20220811.csv',data.table = F)
rownames(df_nor)<-df_nor[,1]
df_nor<-df_nor[,-1]
df_nor1 <- df_nor[-grep("pool", rownames(df_nor)),]
df_nor1$matter<-info$W_G[match(rownames(df_nor1),info$sample)]
unique(df_nor1$matter) #"G" "W" NA
df_nor2<-df_nor1[grep("G",df_nor1$matter),]
df_nor3<-df_nor2[!is.na(df_nor2$matter),]
df_nor3<-df_nor3[,-ncol(df_nor3)]
df_nor3$pathology <- info$pathology[match(rownames(df_nor3),info$sample)]
df_nor4 <- df_nor3[which(df_nor3$pathology=="health"), ]
df_nor4 <- df_nor4[, -ncol(df_nor4)]

df_nor4$BBA <- info$BBA_ABB[match(rownames(df_nor4), info$sample)]

#proteins for antibody
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
prolist <- read.csv("proteins_for_antibody.csv")
prolist1 <- colnames(df_nor4)[unlist(sapply(prolist$Uniprot, function(x){grep(x, colnames(df_nor4))}))]
prolistFun <- FuncPro1[unlist(sapply(prolist$Uniprot, function(x){grep(x, FuncPro1$pro_iso)})),]

#data for heat-map
pro <- unique(prolistFun$pro_iso)
datHeatall <- c(1)
for(i in 1:length(pro)){
  d1 <- df_nor4[,c(which(colnames(df_nor4)==pro[i]), ncol(df_nor4))]
  datHeat <- NULL
  for(j in 1:length(func)){
    region <- unique(FuncBBA2$BBA[which(FuncBBA2$BD==func[j])])
    d2 <- d1[unlist(sapply(region, function(x){which(d1$BBA==x)})), ]
    d2$person <- info$Donor_ID[match(rownames(d2), info$sample)]
    d3 <- d2[, c(1,3)]
    
    #mean
    d4 <- melt(d3)
    d4.1 <- dcast(d4, person~variable, mean)
    d4 <- as.data.frame(t(d4.1))
    colnames(d4) <- d4[1,]
    d4 <- d4[-1,]
    rownames(d4) <- func[j]
    datHeat <- rbind(datHeat, d4)
  }
  colnames(datHeat) <- paste0(pro[i], "_", colnames(datHeat))
  datHeatall <- cbind(datHeatall, datHeat)
}
datHeatall <- datHeatall[,-1]
datHeatall2 <- apply(datHeatall, 2, function(x){as.numeric(x)})
rownames(datHeatall2) <- rownames(datHeatall)

#scale proteins
datHeatall3 <- as.data.frame(scale(datHeatall2))

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
process_order <- read.csv("process_order.csv", header=FALSE)
datHeatall3 <- datHeatall3[process_order$V1, ]

# #scale person and function
# datHeatall3 <- NULL
# for(i in 1:length(pro)){
#   d1 <- datHeatall2[,grep(pro[i],colnames(datHeatall2))]
#   d2 <- as.data.frame(t(d1))
#   d3 <- scale(d2)
#   datHeatall3 <- rbind(datHeatall3, d3)
# }
# datHeatall3 <- as.data.frame(t(datHeatall3))

#heatmap
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
# write.csv(pro, "functions_proteins_330.csv")
library(pheatmap)
pheatmap(datHeatall3,
         fontsize=15,
         color = colorRampPalette(c("#FFFF33","#F4C70B","#AD4462","#67277C"))(10),  #"#FFFF33","#F4C70B"
         # annotation_col = annotation_col,
         # annotation_colors = list(W_G=c(rep("#4472C4",26),rep("#EE781C",19),Region=rainbow(45))),
         # annotation_row = annotation_row,
         # annotation_colors = annotation_colors,
         # na_color = "white",
         gaps_col = seq(0,ncol(datHeatall2),by=4),
         border_color=NA,    #"#EFEFEF",
         trace = "none", 
         scale = "none", 
         cluster_rows = F,
         cluster_cols = F,
         show_rownames = T, 
         show_colnames = T, 
         fontsize_col=5,
         fontsize_row=5,
         filename = "Heatmap_log2_proteins_intensity_mean_function_BBA_scale_pro_20230813.pdf",  #_scale_person_function
         # cellwidth=4,
         # cellheight = 6,
         # main = "",
         height = 7,width = 15
         )

#mean expression in BBA region
df_nor5 <- df_nor4[,pro]   #df_nor4[,pro]  df_nor4[,unique(prolist1)]
df_nor5$BBA <- info$BBA_Anatomical.and.modified.Cyto.architectonic.descriptions[match(rownames(df_nor5), info$sample)]
BBA <- unique(df_nor5$BBA)
datall <- NULL
for(i in 1:length(BBA)){
  d1 <- df_nor5[which(df_nor5$BBA==BBA[i]),]
  d2 <- d1[,-ncol(d1)]
  d2$person <- info$Donor_ID[match(rownames(d2), info$sample)]
  d3 <- melt(d2)
  d4 <- dcast(d3, person~variable, mean)
  d5 <- melt(d4)
  d5$ProPer <- paste0(d5$variable, "_", d5$person)
  d5$BBA <- BBA[i]
  datall <- rbind(datall, d5)
}
datall2 <- datall[,c(4,5,3)]
datall3 <- dcast(datall2, ProPer~BBA)
rownames(datall3) <- datall3[,1]
datall3 <- datall3[,-1]
datall4 <- as.data.frame(t(datall3))
# setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
# write.csv(datall4, "datheatmap_BBA_mean.csv")

#scale by proteins
datall4[is.na(datall4)] <- min(datall4, na.rm=T)/2
# datall4 <- scale(datall4)

#heatmap
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
# write.csv(pro, "functions_proteins_330.csv")
library(pheatmap)
pheatmap(datall4,
         fontsize=15,
         color = colorRampPalette(c("#FFFF33","#F4C70B","#AD4462","#67277C"))(10),  #"#FFFF33","#F4C70B"
         # annotation_col = annotation_col,
         # annotation_colors = list(W_G=c(rep("#4472C4",26),rep("#EE781C",19),Region=rainbow(45))),
         # annotation_row = annotation_row,
         # annotation_colors = annotation_colors,
         na_color = "white",
         gaps_col = seq(0,ncol(datall4),by=4),
         border_color=NA,    #"#EFEFEF",
         trace = "none", 
         scale = "none", 
         cluster_rows = T,
         cluster_cols = F,
         show_rownames = T, 
         show_colnames = T, 
         fontsize_col=5,
         fontsize_row=3,
         filename = "Heatmap_log2_proteins_intensity_mean_BBA_no_scale_cluster_rows_T_20230814.pdf",  #_scale_27_pro
         # cellwidth=4,
         # cellheight = 6,
         # main = "",
         height = 7,width = 15
         )
```

select samples for antibody
```{r}
#the data of proteins expression
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/WGCNA_analysis")
info<-read.csv('FIleInfo_HBA_20220816_withoutFALSE_and_pool.csv')
info$BBA_ABB <- sapply(info$BBA_Anatomical.and.modified.Cyto.architectonic.descriptions, function(x){str_split(x, ",")[[1]][1]})

df_nor<-data.table::fread('HBA_prot_matrix_normolized_20220811.csv',data.table = F)
rownames(df_nor)<-df_nor[,1]
df_nor<-df_nor[,-1]
df_nor1 <- df_nor[-grep("pool", rownames(df_nor)),]
df_nor1$matter<-info$W_G[match(rownames(df_nor1),info$sample)]
unique(df_nor1$matter) #"G" "W" NA
df_nor2<-df_nor1[grep("G",df_nor1$matter),]
df_nor3<-df_nor2[!is.na(df_nor2$matter),]
df_nor3<-df_nor3[,-ncol(df_nor3)]
df_nor3$pathology <- info$pathology[match(rownames(df_nor3),info$sample)]
df_nor4 <- df_nor3[which(df_nor3$pathology=="health"), ]
df_nor4 <- df_nor4[, -ncol(df_nor4)]
df_nor4$BBA <- info$BBA_ABB[match(rownames(df_nor4), info$sample)]

pro <- c("Q9Y291", "Q13405", "Q9BSW7")
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
region <- read.csv("select_regions_for_antibody.csv", header=FALSE)
df <- df_nor4[, pro]
df$BBA <- df_nor4$BBA
df1 <- df[unlist(sapply(region$V1, function(x){which(df$BBA==x)})),]
df1$sample <- rownames(df1)

info1 <- info[,c("sample", "SampleID", "Donor_ID", "Hemisphere")]
df2 <- merge(df1, info1, by="sample")

donor <- c("P3", "P4")
df3 <- df2[unlist(sapply(donor, function(x){which(df2$Donor_ID==x)})), ]
df3 <- df3[order(df3$BBA), ]
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/pro_list_for_antibody")
write.csv(df3, "select_process_regions_info.csv")
```


```{r}
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
FuncBBA <- read.csv("BDf_FDR05.csv")

func <- c("Emotion", "Preparation", "Inhibition")

SF <- FuncBBA[unlist(sapply(func, function(x){grep(x, FuncBBA$BD)})), c(2,4)]
SF <- unique(SF)
SF <- SF[order(SF$BD),]

setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
write.csv(SF, "select_process_regions.csv")

emotion <- FuncBBA[grep("Emotion", FuncBBA$BD), c(2,4)]
emotion <- unique(emotion)
region <- unique(emotion$BBA)
num <- NULL
for(i in 1:length(region)){
  d1 <- emotion[which(emotion$BBA==region[i]), ]
  num0 <- data.frame(region=region[i], num_funtion = length(unique(d1$BD)))
  num <- rbind(num, num0)
}
num <- num[order(num$num_funtion, decreasing = T),]
setwd("D:/1_guomics_xyt/1Project/hba_research_workflow/9_DataAnalysis_20220901/20230423_new_figures/20230531_brain_function_heatmap")
write.csv(num, "BBA_emotion_num.csv")
```