单细胞基因的相关性计算--以pbmc3k为例，计算S100A8正相关和负相关的基因

参考：https://www.jianshu.com/p/9850c4d88a67

# 1.以pbmc3k为例，加载Seurat数据+设置工作目录----#####
setwd("/public/home/chidm/test/Gene相关性")
load("/public/home/chidm/test/pbmc3k/pbmc3k_final.RData")
pbmc
 
exprSet <- pbmc@assays$RNA@data
exprSet<-as.data.frame(t(exprSet))#转置

# 2. 计算某个基因和其它基因的相关性（以S100A8为例）-----#####
y <- as.numeric(exprSet[,"S100A8"])
colnames<-colnames(exprSet)
cor_data_df<- data.frame(colnames)
for(i in 1:length(colnames)){
  test<-cor.test(as.numeric(exprSet[,i]),y,type="spearman") #可更换pearson
  cor_data_df[i,2]<- test$estimate
  cor_data_df[i,3]<- test$p.value
  
}
names(cor_data_df)<-c("symbol","correlation","pvalue")
cor_data_df %>% head()

# 3. 筛选有意义的正相关和负相关的基因-----####
library(dplyr)
library(tidyr)
cor_data_sig_pos <- cor_data_df %>%
  dplyr::filter(pvalue <0.01)%>%dplyr::filter(correlation >0)%>%
  dplyr::arrange(desc(correlation))
 
cor_data_sig_neg <- cor_data_df %>%
  dplyr::filter(pvalue <0.01)%>%dplyr::filter(correlation <0)%>%
  dplyr::arrange(desc(abs(correlation)))

# 4. 随机选取正相关和负相关基因，分别作图验证----######
#1）S100A9正相关####
# install.packages("ggstatsplot") ##linux安装过程中报错，需先安装mpfr工具
# $conda install anaconda::mpfr
library(ggstatsplot)
png(paste0("相关性densigram-正相关.png"), width = 6, height = 6, res = 400, units = "in")
ggscatterstats(data = exprSet ,
               y = S100A8,
               x = S100A9,
               centrality.para = "mean",               
               margins = "both",
               xfill = "#CC79A7",
               yfill = "#009E73",
               marginal.type = "densigram", # #类型可以换成density,boxplot,violin,densigram                
               title = "Relationship between S100A8 and S100A9")
dev.off()
#2）S100A9负相关####
library(ggstatsplot)
png(paste0("相关性densigram-负相关.png"), width = 6, height = 6, res = 400, units = "in")
ggscatterstats(data = exprSet,
               y = MALAT1,
               x = S100A9,
               centrality.para = "mean",               
               margins = "both",
               xfill = "#CC79A7",
               yfill = "#009E73",
               marginal.type = "densigram", # #类型可以换成density,boxplot,violin,densigram                
               title = "Relationship between MALAT1 and S100A9")
dev.off()
 
# 这两个图的更简单画法（Seurat的FeatureScatter函数)
library(Seurat)
library(SeuratObject)
png(paste0("相关性-FeatureScatter.png"), width = 8, height = 4, res = 400, units = "in")
p1=FeatureScatter(pbmc, feature1 = "S100A8", feature2 = "S100A9")
p2=FeatureScatter(pbmc, feature1 = "S100A8", feature2 = "MALAT1")
p1|p2
dev.off()