参考1、后续分析 https://mp.weixin.qq.com/s/0A9hiXqV1RfWz8UAnratAA
参考2、一文打通单细胞转录组细胞类型丰度变化分析https://mp.weixin.qq.com/s/DlUGRNv4gi8Tp9OztSN-0A
参考3、详细解释1https://mp.weixin.qq.com/s/D0nGq3XGZRIlso7dKo0qkQ
参考4、详细解释2https://mp.weixin.qq.com/s/zSEsjWb0Uz3p5wAd7Gu61g
跑这个函数sce_milo <- calcNhoodDistance(sce_milo, d=50)的时候会出现这个报错：Error: (converted from warning) as(<dgTMatrix>, "dgCMatrix")不再有用。 官方给的解决办法：不用跑这步。
代码如下：
BiocManager::install("miloR")
# 差异表达分析我们都很熟悉，就是看不同组之间、或者实验处理后与对照组相比基因表达的变化。
# 而差异丰度分析是检测不同条件下celltype组成的变化，通过不同组、实验条件、发育阶段、
# 时间阶段比较，发现新的细胞亚群或者细胞比例显著差异的亚群。
#setwd('D:/KS项目/公众号文章/miloR单细胞差异丰度')
#devtools::install_github("MarioniLab/miloR", ref="devel")
#devtools::install_github("MarioniLab/miloR", ref="devel")
library(miloR)#UMAP降维数据不会出错，TSNE不行,要有重复
library(Seurat)
library(ggplot2)
library(SingleCellExperiment)
#remotes::install_github('satijalab/seurat-wrappers')
library(SeuratWrappers)
library(ggbeeswarm)
library(scater)
library(scales)
library(forcats)
library(data.table)
library(stringr)
library(dplyr)
#================================================================================
# 只有两组（例如normal组和疾病组）
#================================================================================
### rm(list=ls())
#rm(immune.5)
immune.5=immune.2
### dev.new()
#加载单细胞数据
#load("D:/KS项目/公众号文章/miloR单细胞差异丰度/immune.5.RData")
DimPlot(immune.5,label = T)
unique(immune.5$group)
# [1] "disease" "healthy"
immune.5$group<- factor(immune.5$group, levels=c("HC","Act"))#样本分组
#miloR输入对象是SingleCellExperiment，所以我们是常用的seurat对象的话，转化一下
immune.5 <- as.SingleCellExperiment(immune.5)
immune.5_milo <- miloR::Milo(immune.5) #milo object构建
immune.5_milo <- miloR::buildGraph(immune.5_milo, k = 30, d = 50) #构建k-nearest neighbour graph
#k是构建graph是需要考虑的邻近的数量，这个参数可自行调整
#定义具有代表性的邻域
immune.5_milo <- makeNhoods(immune.5_milo,
prop = 0.1, #随机抽取洗吧比例，一般情况0.1-0.2足够
k = 30, #建议使用与buildGraph一样的k值
d=50, #KNN降维数，建议使用与buildGraph一样的d值
refined = TRUE,
refinement_scheme="graph")
immune.5_milo <- countCells(immune.5_milo, meta.data = data.frame(colData(immune.5_milo)),
sample="orig.ident")#sample是样本重复
#分组，类似于做差异分析的时候构建的design
traj_design <- data.frame(colData(immune.5_milo))[,c("orig.ident", "group")]#分别是重复样本ID和分组
traj_design$orig.ident <- as.factor(traj_design$orig.ident)
traj_design <- distinct(traj_design)
rownames(traj_design) <- traj_design$orig.ident
#计算cell距离
#immune.5_milo <- calcNhoodDistance(immune.5_milo, d=50)
#差异分析结果、实验分组两者之前差异
da_results <- testNhoods(immune.5_milo,
design = ~ group,
design.df = traj_design,
fdr.weighting="graph-overlap") #分组
immune.5_milo <- buildNhoodGraph(immune.5_milo) #构建graph for 可视化
#save(immune.5_milo, file = "immune.5_milo-Covid-HC-1128.RData")
#rm(immune.5_milo)
#可视化，这些可视化函数都可以与ggplot2互通，所以图片修饰上可以按照ggplot2的方式
#可视化1----------------------------------------------------------------------------
#UMAP图
#plotUMAP(immune.5_milo, colour_by = "celltype")
p <- plotNhoodGraphDA(immune.5_milo, da_results,layout="UMAP", alpha=0.1) +
scale_fill_gradient2(low="#070091",#修改颜色
mid="lightgrey",
high="#910000",
name="log2FC",
limits=c(-5,5),
oob=squish)
ggsave("milo_PMN_3.pdf", p, width = 10, height = 8, units = "in")
#可视化2----------------------------------------------------------------------------
#蜂群图展示celltype logFC变化
da_results <- annotateNhoods(immune.5_milo, da_results, coldata_col = "celltype")
p1 <- plotDAbeeswarm(da_results, group.by = "celltype") +
scale_color_gradient2(low="#070091",
mid="lightgrey",
high="#910000",
limits=c(-5,5),
oob=squish) +
labs(x="", y="Log2 Fold Change") +
theme_bw(base_size=10)+
theme(axis.text = element_text(colour = 'black'))
ggsave("milo_蜂群图_3.pdf", p1, width = 10, height = 8, units = "in")
