单细胞测序两组差异分析

20210829修改
之前是根据官网+别人帖子写的总结，自己做了一段时间，把之前的再完善一下

尝试使用seurat包进行两组间差异分析
使用的是seurat包自带的数据

创建seurat对象

#首先载入需要的包
library(Seurat)
#安装seurat-data包
install.packages('devtools')
library("devtools")
devtools::install_github('satijalab/seurat-data')
library(SeuratData)
library(patchwork)

安装数据集

# install dataset
InstallData("ifnb")

出现以下报错，可能与网速有关
然后我尝试手动下载了这个包，网址为：http://seurat.nygenome.org/src/contrib/ifnb.SeuratData_3.1.0.tar.gz
之后在Rstudio中进行手动安装，可参考这篇文章https://blog.csdn.net/lym152898/article/details/77572163/
安装成功！撒花！
奇奇怪怪，LoadData之后进行分组，出现报错

之后我检查了ifnb.SeuratData包的安装情况，确实是已经安装了的
手动打勾，进行library这个包（官方文档中并没有这句代码），但却成了

# 载入数据集
LoadData("ifnb")# 将数据集分成两组
ifnb.list <- SplitObject(ifnb, split.by = "stim")# 标化数据并找到数据的特征分子
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {x <- NormalizeData(x)x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})# 选择不同数据集重复的特征分子进行整合
features <- SelectIntegrationFeatures(object.list = ifnb.list)

进行整合（消除批次效应）

进一步使用 FindIntegrationAnchors()确定anchors（锚点，用于整合两个数据），并使用anchors与IntegrateData()整合这两个数据集

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features)immune.combined <- IntegrateData(anchorset = immune.anchors)

整合分析

我们对修正后的数据进行下游分析，原始数据仍然存在于“RNA”assay中

DefaultAssay(immune.combined) <- "integrated"#运行标准化可视化与聚类
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)#可视化
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
p1 + p2

为了根据分组进行可视化，使用split.by 命令

DimPlot(immune.combined, reduction = "umap", split.by = "stim")

鉴定保守细胞的特征标记物

FindConservedMarkers() 用于鉴定不受处理条件影响的细胞marker，该函数为每个数据集/组执行差异基因表达测试，并使用MetaDE R包中的元分析方法结合p值。例如，可以计算出在簇6 (NK细胞)中，那些不受刺激条件影响的保守标记基因。

#为了进行整合后的差异表达，我们又回到了原始数据

DefaultAssay(immune.combined) <- "RNA"
#需要先安装包
install.packages('BiocManager')
BiocManager::install('multtest')
install.packages('metap')nk.markers <- FindConservedMarkers(immune.combined, ident.1 = 6, grouping.var = "stim", verbose = FALSE)
head(nk.markers)

我们可以对这些marker基因进行可视化，并使用它们注释不同细胞

#可视化标记基因
FeaturePlot(immune.combined, features = c("CD3D", "SELL", "CREM", "CD8A", "GNLY", "CD79A", "FCGR3A",  "CCL2", "PPBP"), min.cutoff = "q9")immune.combined <- RenameIdents(immune.combined, `0` = "CD14 Mono", `1` = "CD4 Naive T", `2` = "CD4 Memory T", `3` = "CD16 Mono", `4` = "B", `5` = "CD8 T", `6` = "NK", `7` = "T activated", `8` = "DC", `9` = "B Activated", `10` = "Mk", `11` = "pDC", `12` = "Eryth", `13` = "Mono/Mk Doublets", `14` = "HSPC")DimPlot(immune.combined, label = TRUE)

使用DotPlot函数对marker基因进行可视化，其中split.by参数可用于查看不同条件下的保守细胞类型marker基因，显示特定细胞类群中基因的表达水平和细胞的百分比。在这里，我们展示了13个细胞类群中每个群中的2-3个强标记基因。

Idents(immune.combined) <- factor(Idents(immune.combined), levels = c("HSPC", "Mono/Mk Doublets", "pDC", "Eryth", "Mk", "DC", "CD14 Mono", "CD16 Mono", "B Activated", "B", "CD8 T", "NK", "T activated", "CD4 Naive T", "CD4 Memory T"))
markers.to.plot <- c("CD3D", "CREM", "HSPH1", "SELL", "GIMAP5", "CACYBP", "GNLY", "NKG7", "CCL5", "CD8A", "MS4A1", "CD79A", "MIR155HG", "NME1", "FCGR3A", "VMO1", "CCL2", "S100A9", "HLA-DQA1", "GPR183", "PPBP", "GNG11", "HBA2", "HBB", "TSPAN13", "IL3RA", "IGJ", "PRSS57")DotPlot(immune.combined, features = markers.to.plot, cols = c("blue", "red"), dot.scale = 8, split.by = "stim") + RotatedAxis()

鉴定不同刺激处理下的差异表达基因

接下来，我们将对不同刺激处理下的两组数据集进行比较分析。我们使用的方法是计算刺激组和对照组中细胞的平均表达，并绘制散点图识别那些异常表达的基因。在这里，我们计算了刺激组和对照组中naive T细胞和CD14单核细胞群体的平均表达并生成散点图，突出显示那些对干扰素刺激表现出强烈反应的基因。

library(ggplot2)
library(cowplot)
theme_set(theme_cowplot())
t.cells <- subset(immune.combined, idents = "CD4 Naive T")
Idents(t.cells) <- "stim"
avg.t.cells <- as.data.frame(log1p(AverageExpression(t.cells, verbose = FALSE)$RNA))
avg.t.cells$gene <- rownames(avg.t.cells)cd14.mono <- subset(immune.combined, idents = "CD14 Mono")
#好像是分组的作用，如果不加这句，运行后面求平均值会出现警告，只有一组
Idents(cd14.mono) <- "stim"
avg.cd14.mono <- as.data.frame(log1p(AverageExpression(cd14.mono, verbose = FALSE)$RNA))
avg.cd14.mono$gene <- rownames(avg.cd14.mono)genes.to.label = c("ISG15", "LY6E", "IFI6", "ISG20", "MX1", "IFIT2", "IFIT1", "CXCL10", "CCL8")
p1 <- ggplot(avg.t.cells, aes(CTRL, STIM)) + geom_point() + ggtitle("CD4 Naive T Cells")
p1 <- LabelPoints(plot = p1, points = genes.to.label, repel = TRUE)
p2 <- ggplot(avg.cd14.mono, aes(CTRL, STIM)) + geom_point() + ggtitle("CD14 Monocytes")
p2 <- LabelPoints(plot = p2, points = genes.to.label, repel = TRUE)
p1 + p2

如我们所看到，在两种细胞中有相同基因上调，这些基因可能代表了保守的干扰素应答途径。

首先我们在原数据中创建了一个列，包括细胞类型和刺激信息，并将该列定义为数据的ident。之后我们使用FindMarkers()来发现处理组与对照组B细胞中的不同基因。注意，这里显示的很多顶级基因与我们之前绘制的核心干扰素反应基因相同。此外，我们所看到的CXCL10是针对单核细胞和B细胞干扰素应答的，在列表中显示出重要性。

#加了一列细胞类型及刺激状态
immune.combined$celltype.stim <- paste(Idents(immune.combined), immune.combined$stim, sep = "_")
#加了一列细胞类型
immune.combined$celltype <- Idents(immune.combined)
#将细胞类型及刺激状态作为分组
Idents(immune.combined) <- "celltype.stim" #使用FindMarkers函数寻找差异表达基因
b.interferon.response <- FindMarkers(immune.combined, ident.1 = "B_STIM", ident.2 = "B_CTRL", verbose = FALSE)
head(b.interferon.response, n = 15)

使用FeaturePlot和VlnPlot函数可视化标记基因，并设置split.by按分组变量（此处为刺激条件）进行划分。其中，CD3D和GNLY是T细胞和NK / CD8 T细胞典型的marker，不受干扰素刺激的影响，在对照组和受刺激组中显示出相似的基因表达模式。另一方面，IFI6和ISG15是核心的干扰素响应基因，因此在所有细胞类型中均被上调。而CD14和CXCL10是干扰素应答基因。

FeaturePlot(immune.combined, features = c("CD3D", "GNLY", "IFI6"), split.by = "stim", max.cutoff = 3, cols = c("grey", "red"))
plots <- VlnPlot(immune.combined, features = c("LYZ", "ISG15", "CXCL10"), split.by = "stim", group.by = "celltype", pt.size = 0, combine = FALSE)
wrap_plots(plots = plots, ncol = 1)

在seurat官方文档中还介绍了一种标化方法，这种方法待我将前面的流程完全搞明白后再更新（这种方法我暂时还没有看-_-）

对于assay的解释
参考：https://www.bilibili.com/read/cv6411199/
直接输入Seurat object的名称，我们可以得到类似如下内容：
An object of class Seurat 13425 features across 39233 samples within 1 assay

Active assay: RNA (13425 features, 3000 variable features) 3 dimensional reductions calculated: pca, umap, tsne

这个告诉我们当前对象主体是13425（基因数）*39233（细胞数）的矩阵，有一个叫RNA的assay，在这个assay中，我们选择了3000个基因作为variable features（高表达基因，认为这种变化的基因分析起来更具有代表性，一般用来计算PCA），计算了三种降维：PCA, UMAP, t-SNE。

默认情况下，我们的seurat对象中是一个叫RNA的Assay。在我们处理数据的过程中，做整合（integration），或者做变换（SCTransform），或者做去除污染（SoupX），或者是融合velocity的数据等，我们可能会生成新的相关的Assay，用于存放这些处理之后的矩阵。

在之后的处理中，我们可以根据情况使用指定Assay下的数据。不指定Assay使用数据的时候， Seurat给我们调用的是Default Assay下的内容。可以通过对象名@active.assay查看当前Default Assay，通过DefaultAssay函数更改当前Default Assay。Assay数据中，counts为raw，data为normalized，scale为scaled。

1.调用Assay中的数据的方式为，以调取一个名为PBMC的Seurat对象中Assay integrate中的nomalized数据为例：
PBMC@assays$RNA@data

2.meta.data
元数据，对每个细胞的描述。一般计算的nFeature_RNA等信息就以metafeature的形式存在Seurat对象的metadata中。计算的分类信息一般以RNA_snn_res.x（x指使用的resolution）存放在metadata中。

调取metadata中metafeature值的方式有多种，以调取一个名为PBMC的对象中stim这个metafeature为例：

方法1：PBMC[[“stim”]]
方法2：PBMC$stim

3.reductions

降维之后的每个细胞的坐标信息。

以调取一个名为PBMC的对象中PCA embedding (也就是坐标)信息为例：
PBMC@reductions$pca@cell.embeddings

ownames(object) 获取的是全部基因
colnames(object)获取的是全部细胞id
VariableFeatures(object)获取当前object的Variable feature
levels(object)获取当前object的分类信息