nichenet实战silicosis

############################################https://www.jianshu.com/p/5db596333dfb
#https://cloud.tencent.com/developer/article/1819551
.libPaths()
.libPaths("G:/R_big_packages/")
.libPaths(c("G:/R_big_packages/","D:/Win10 System/Documents/R/win-library/4.1","C:/Program Files/R/R-4.1.0/library"))
.libPaths()

installed.packages()[, c("Package", "LibPath")]
Sys.chmod(.libPaths(),'777')
#library(devtools)
#install_github("saeyslab/nichenetr")

#https://www.jianshu.com/p/c0f3c3728a8d
#如果不能正常加载包的话,就首先unloadNamespace这些包！！！！！！！！！！！！！！！！！！！！！！！！！！！！
for (packagename in c("dplyr","recipes","tidyselect",
                      "dbplyr","tidyr","ggpubr","ggplot2","usethis","broom")) {
    
  unload(packagename)
}


library(nichenetr)
library(Seurat) # please update to Seurat V4
library(tidyverse)
library(ggplot2)

#BiocManager::install("glue")
#library("glue")
sessionInfo()





'''
NicheNet分析实践
数据来源
本文的分析数据和代码来自NicheNet官方分析单细胞数据的教程，
https://github.com/saeyslab/nichenetr/blob/master/vignettes/seurat_wrapper.md
演示数据集源自Medaglia et al. 2017 “Spatial Reconstruction of Immune Niches by Combining Photoactivatable Reporters and scRNA-Seq.” Science, December, eaao4277.
https://doi.org/10.1126/science.aao4277.
我们将使用Medaglia等人的小鼠NICHE-seq数据，探索淋巴细胞性脉络膜脑膜炎病毒（LCMV）感染之前
和之后72小时的腹股沟淋巴结T细胞区域的细胞间通讯。
在该数据集中，观察到稳态下的CD8 T细胞与LCMV感染后的CD8 T细胞之间存在差异表达。NicheNet可用于观察淋巴结中的几种免疫细胞群
（即单核细胞，树突状细胞，NK细胞，B细胞，CD4 T细胞）如何调节和诱导这些观察到的基因表达变化。

'''


getwd()
path="G:/silicosis/sicosis/NicheNet/6-30_silisosis"
dir.create(path)
setwd(path)
getwd()

1#读入单细胞数据
#seuratObj = readRDS(url("https://zenodo.org/record/3531889/files/seuratObj.rds"))
#save(seuratObj,file = "G:/silicosis/sicosis/NicheNet/seuratobj.rds")
load(file = "G:/silicosis/sicosis/NicheNet/seuratobj.rds")
[e-mail] %>% head()    ##aggregate是处理条件，SS相当于control，LCMV相当于case。
seuratObj@assays$RNA@data[1107:1108,501:505]

2##读入nichenet先验数据 
#ligand_target_matrix <- readRDS(url("https://zenodo.org/record/3260758/files/ligand_target_matrix.rds"))
#save(ligand_target_matrix,file ="G:/silicosis/sicosis/NicheNet/ligand_target_matrix.rds" )
load(file ="G:/silicosis/sicosis/NicheNet/ligand_target_matrix.rds")
head(ligand_target_matrix)[1:4,1:5] #说明是人源的数据

#lr_network = readRDS(url("https://zenodo.org/record/3260758/files/lr_network.rds"))
#save(lr_network,file = "G:/silicosis/sicosis/NicheNet/lr_network.rds")
load(file = "G:/silicosis/sicosis/NicheNet/lr_network.rds")
head(lr_network)

#weighted_networks = readRDS(url("https://zenodo.org/record/3260758/files/weighted_networks.rds"))
#save(weighted_networks,file ="G:/silicosis/sicosis/NicheNet/weighted_networks.rds" )
load(file = "G:/silicosis/sicosis/NicheNet/weighted_networks.rds")
weighted_networks_lr = weighted_networks$lr_sig %>% inner_join(lr_network %>% distinct(from,to), by = c("from","to"))
head(weighted_networks$lr_sig) # interactions and their weights in the ligand-receptor + signaling network
# weighted_networks列表包含两个数据框，lr_sig是配体-受体权重信号网络，gr是配体-靶基因权重调控网络


3 ### 
'''
NicheNet分析
我们希望使用NicheNet预测哪些配体可能影响CD8 T细胞在LCMV感染后的差异表达基因。
此例中‘CD8 T cell’是receiver细胞，‘CD4 T’, ‘Treg’, ‘Mono’, ‘NK’, ‘B’ and ‘DC’是sender细胞。
NicheNet提供了一个打包函数nichenet_seuratobj_aggregate，它可以一步完成seurat对象的配体调控网络分析。
'''
scRNA=seuratObj
table(scRNA$celltype,scRNA$aggregate)
Idents(scRNA)

Idents(scRNA) <- "celltype"
nichenet_output = nichenet_seuratobj_aggregate(seurat_obj = scRNA, 
                                               top_n_ligands = 20,
                                               receiver = "CD8 T", 
                                               sender = c("CD4 T","Treg", "Mono", "NK", "B", "DC"),
                                               condition_colname = "aggregate", 
                                               condition_oi = "LCMV", 
                                               condition_reference = "SS", 
                                               ligand_target_matrix = ligand_target_matrix, 
                                               lr_network = lr_network, 
                                               weighted_networks = weighted_networks, 
                                               organism = "mouse")
# top_n_ligands参数指定用于后续分析的高活性配体的数量             
#save(nichenet_output,file = "nichenet_output.rds")
load(file = "G:/silicosis/sicosis/NicheNet/nichenet_output.rds")


getwd()
4#NicheNet结果
## 查看配体活性分析结果
# 主要参考pearson指标，bona_fide_ligand=True代表有文献报道的配体-受体，

# bona_fide_ligand=False代表PPI预测未经实验证实的配体-受体。
head(nichenet_output)
names(nichenet_output)
x <- nichenet_output$ligand_activities
write.csv(x, "ligand_activities.csv", row.names = F)

5#
# 查看top20 ligands
nichenet_output$top_ligands
# 查看top20 ligands在各个细胞亚群中表达情况
p = DotPlot(scRNA, features = nichenet_output$top_ligands, cols = "RdYlBu") + RotatedAxis()
ggsave("top20_ligands.png", p, width = 12, height = 6)
# 按"aggregate"的分类对比配体的表达情况
p = DotPlot(scRNA, features = nichenet_output$top_ligands, split.by = "aggregate") + RotatedAxis()
ggsave("top20_ligands_compare.png", p, width = 12, height = 8)
# 用小提琴图对比配体的表达情况
p = VlnPlot(scRNA, features = nichenet_output$top_ligands, 
            split.by = "aggregate",split.plot = TRUE,  pt.size = 0)
ggsave("VlnPlot_ligands_compare.png", p, width = 12, height = 8)

6#
## 查看配体调控靶基因
p = nichenet_output$ligand_target_heatmap
ggsave("Heatmap_ligand-target.png", p, width = 12, height = 6)
# 更改热图的风格
p = nichenet_output$ligand_target_heatmap + 
  scale_fill_gradient2(low = "whitesmoke",  high = "royalblue", breaks = c(0,0.0045,0.009)) + 
  xlab("anti-LCMV response genes in CD8 T cells") + 
  ylab("Prioritized immmune cell ligands")
ggsave("Heatmap_ligand-target2.png", p, width = 12, height = 6)
# 查看top配体调控的靶基因及其评分
x <- nichenet_output$ligand_target_matrix
#x2 <- nichenet_output$ligand_target_df
write.csv(x, "ligand_target.csv", row.names = F)
# 查看被配体调控靶基因的表达情况
p = DotPlot(scRNA %>% subset(idents = "CD8 T"), 
            features = nichenet_output$top_targets, 
            split.by = "aggregate") + RotatedAxis()
ggsave("Targets_Expression_dotplot.png", p, width = 12, height = 6)
p = VlnPlot(scRNA %>% subset(idents = "CD8 T"), features = nichenet_output$top_targets, 
            split.by = "aggregate", pt.size = 0, combine = T, ncol = 8)
ggsave("Targets_Expression_vlnplot.png", p, width = 12, height = 8)


7#
## 查看受体情况
# 查看配体-受体互作
p = nichenet_output$ligand_receptor_heatmap
ggsave("Heatmap_ligand-receptor.png", p, width = 12, height = 6)
x <- nichenet_output$ligand_receptor_matrix
#x <- nichenet_output$ligand_receptor_df
write.csv(x, "ligand_receptor.csv", row.names = F)
# 查看受体表达情况
p = DotPlot(scRNA %>% subset(idents = "CD8 T"), 
            features = nichenet_output$top_receptors, 
            split.by = "aggregate") + RotatedAxis()
ggsave("Receptors_Expression_dotplot.png", p, width = 12, height = 6)
p = VlnPlot(scRNA %>% subset(idents = "CD8 T"), features = nichenet_output$top_receptors, 
            split.by = "aggregate", pt.size = 0, combine = T, ncol = 8)
ggsave("Receptors_Expression_vlnplot.png", p, width = 12, height = 8)
# 有文献报道的配体-受体
# Show ‘bona fide’ ligand-receptor links 
#that are described in the literature and not predicted based on PPI
p = nichenet_output$ligand_receptor_heatmap_bonafide
ggsave("Heatmap_ligand-receptor_bonafide.png", p, width = 8, height = 4)
x <- nichenet_output$ligand_receptor_matrix_bonafide
#x <- nichenet_output$ligand_receptor_df_bonafide
write.csv(x, "ligand_receptor_bonafide.csv", row.names = F)