练习数据集的GEO编号：GSE139324 (Immune Landscape of Viral- and Carcinogen-Driven Head and Neck Cancer)。原数据集有63个scRNA的数据，本次选取10个练习。

library(Seurat)
library(harmony)
library(tidyverse)
library(patchwork)
rm(list = ls())

一. 多个单细胞样本的合并

1. 读取并合并数据

• 1.1 读取数据

## 批量读取数据
### 设置数据路径与样本名称
assays <- dir("GSE139324/")
dir <- paste0("GSE139324/", assays)
# 按文件顺序给样本命名，名称不要以数字开头，中间不能有空格 
samples_name = c('HNC01PBMC', 'HNC01TIL', 'HNC10PBMC', 'HNC10TIL', 'HNC20PBMC',
                 'HNC20TIL',  'PBMC1', 'PBMC2', 'Tonsil1', 'Tonsil2')

• 1.2 批量创建seurat对象

scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  #不设置min.cells过滤基因会导致CellCycleScoring报错：
  #Insufficient data values to produce 24 bins.  
  scRNAlist[[i]] <- CreateSeuratObject(counts, project=samples_name[i],
                                       min.cells=3, min.features = 200)
  #给细胞barcode加个前缀，防止合并后barcode重名
  scRNAlist[[i]] <- RenameCells(scRNAlist[[i]], add.cell.id = samples_name[i])   
  #计算线粒体基因比例
  if(T){    
    scRNAlist[[i]][["percent.mt"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^MT-") 
  }
  #计算核糖体基因比例
  if(T){
    scRNAlist[[i]][["percent.rb"]] <- PercentageFeatureSet(scRNAlist[[i]], pattern = "^RP[SL]")
  }
  #计算红细胞基因比例
  if(T){
    HB.genes <- c("HBA1","HBA2","HBB","HBD","HBE1","HBG1","HBG2","HBM","HBQ1","HBZ")
    HB.genes <- CaseMatch(HB.genes, rownames(scRNAlist[[i]]))
    scRNAlist[[i]][["percent.HB"]]<-PercentageFeatureSet(scRNAlist[[i]], features=HB.genes) 
  }
}

### 给列表命名并保存数据
dir.create("Integrate")
setwd("./Integrate")

names(scRNAlist) <- samples_name
#system.time(save(scRNAlist, file = "Integrate/scRNAlist0.Rdata")) 
system.time(saveRDS(scRNAlist, file = "scRNAlist0.rds"))

这一步得到了一个包含了十个样本Seurat对象的list

• 1.3 使用merge函数将scRNAlist合成一个Seurat对象

scRNA <- merge(scRNAlist[[1]], scRNAlist[2:length(scRNAlist)])
scRNA
# An object of class Seurat 
# 18818 features across 19738 samples within 1 assay 
# Active assay: RNA (18818 features, 0 variable features)
table(scRNA$orig.ident)
# HNC01PBMC  HNC01TIL HNC10PBMC  HNC10TIL HNC20PBMC 
#      1721      1298      1750      1383      1525 
#  HNC20TIL     PBMC1     PBMC2   Tonsil1   Tonsil2 
#      1148      2444      2436      3324      2709 
save(scRNA,file = 'scRNA_orig.Rdata')
#scRNAlist <- SplitObject(scRNA, split.by = "orig.ident") #分割Seurat对象

2. 数据质控

### 绘制质控小提琴图
# 设置可能用到的主题
theme.set2 = theme(axis.title.x=element_blank())
# 设置绘图元素
plot.featrures = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb", "percent.HB")
group = "orig.ident"
# 质控前小提琴图
plots = list()
for(i in seq_along(plot.featrures)){
  plots[[i]] = VlnPlot(scRNA, group.by=group, pt.size = 0,
                       features = plot.featrures[i]) + theme.set2 + NoLegend()}
violin <- wrap_plots(plots = plots, nrow=2)    
dir.create("QC")
ggsave("QC/vlnplot_before_qc.pdf", plot = violin, width = 9, height = 8)

### 设置质控标准
minGene=500
maxGene=4000
maxUMI=15000
pctMT=10
pctHB=1

### 数据质控并绘制小提琴图
scRNA <- subset(scRNA, subset = nCount_RNA < maxUMI & nFeature_RNA > minGene & 
                  nFeature_RNA < maxGene & percent.mt < pctMT & percent.HB < pctHB)
plots = list()
for(i in seq_along(plot.featrures)){
  plots[[i]] = VlnPlot(scRNA, group.by=group, pt.size = 0,
                       features = plot.featrures[i]) + theme.set2 + NoLegend()}
violin <- wrap_plots(plots = plots, nrow=2)     
ggsave("QC/vlnplot_after_qc.pdf", plot = violin, width = 10, height = 8) 

3. 查看批次效应（对merge后的Seurat对象进行标准化和降维）

# 降维聚类
scRNA <- NormalizeData(scRNA) %>% FindVariableFeatures(nfeatures = 3000) %>% ScaleData()
scRNA <- RunPCA(scRNA, verbose = F)
ElbowPlot(scRNA, ndims = 50)

pc.num=1:30
scRNA <- scRNA %>% RunTSNE(dims=pc.num) %>% RunUMAP(dims=pc.num)
scRNA <- FindNeighbors(scRNA, dims=pc.num) %>% FindClusters()
DimPlot(scRNA, label = T)

p <- DimPlot(scRNA, group.by = "orig.ident")
ggsave("UMAP_Samples.pdf", p, width = 8, height = 6)

可以看出，不同样本间还是有批次效应的存在。结合这张图和上面那张图来看，比如上面那张图中的6和7cluster在整合之后应该就是一个cluster

p <- DimPlot(scRNA, group.by = "orig.ident", split.by = "orig.ident", ncol = 4)
ggsave("UMAP_Samples_Split.pdf", p, width = 18, height = 12)

saveRDS(scRNA, "scRNA.rds")

二. 数据整合

数据整合的目的：
1. 画出tsne/umap图，以辅助细胞类型鉴定（把不同组之间的一类细胞整合在一起）
2. 轨迹分析可能需要用到umap图

⚠️数据整合只影响降维聚类，不影响差异分析（原始count矩阵没有改变）。

1. seurat锚点整合

参考：https://satijalab.org/seurat/articles/integration_large_datasets.html
Seurat2的整合主要用的是CCA（canonical correlation analysis，典型关联分析）的方法，Seurat3和Seurat4用的是CCA+MNN（mutual nearest neighbor，互近邻）

锚点整合操作速度很慢，且常常会过度整合，因此在实际操作中，跨物种整合的时候或不同的数据类型如ATAC、蛋白组的数据和单细胞的数据整合的时候，可以使用锚点整合。单纯单细胞数据的整合，可以使用Harmony。

• 1.1 读入数据，拆分样本

rm(list=ls())
library(future) #Seurat并行计算的一个包，不加载这个包不能进行并行计算
options(future.globals.maxSize = 50 * 1024^3) #将全局变量上限调至50G（锚点整合很占内存）

##重新创建没有处理的经过降维等处理的数据 
scRNA <- readRDS("scRNA.rds")
cellinfo <- subset(scRNA@meta.data, select = c("orig.ident", "percent.mt", "percent.rb", "percent.HB"))
scRNA <- CreateSeuratObject(scRNA@assays$RNA@counts, meta.data = cellinfo)
#做锚点整合需要把样本处理成单个的Seurat对象来两两组合
scRNAlist <- SplitObject(scRNA, split.by = "orig.ident")
#也可以按别的指标（metadata中的）来进行拆分，比如可以按不同的分组来拆分样本，再进行整合。

• 1.2 标准化（由于锚点整合会把单个样本两两组合，所以需要单独做标准化）

#SCTransform标准化(⚠️使用log标准化还是SCT标准化差别不大)
#如果用log标准化，后面FindIntegrationAnchors和IntegrateData函数的normalization.method参数选'LogNormalize'
scRNAlist <- parallel::mclapply(scRNAlist, FUN=function(x) SCTransform(x), mc.cores = 10) #10个对象最好写10个核，没有10个核少写几个也可以。top命令可以查看服务器有几个核，mc.core设置为1就每次处理一个对象。
# mclapply是lapply的多核版本

• 1.3 选择用于整合的高变基因（三步）

### FindAnchors 
### 每个样本的高变基因不完全一样，SelectIntegrationFeatures可以整合这些高变基因，选出3000个
scRNA.features <- SelectIntegrationFeatures(scRNAlist, nfeatures = 3000) 
### 将每个样本的高变基因都调整成上一步选出的3000个
scRNAlist <- PrepSCTIntegration(scRNAlist, anchor.features = scRNA.features) 
##寻找锚点，运行速度非常慢，至少需要1-2小时
plan("multisession", workers = 10)
scRNA.anchors <- FindIntegrationAnchors(object.list = scRNAlist,
                                        normalization.method = "SCT",  #如果前面是log标准化，这里改成LogNormalize
                                        anchor.features = scRNA.features)

• 1.4 锚点整合

### Integrate 运行速度慢
scRNA.sct.int <- IntegrateData(scRNA.anchors, normalization.method="SCT") #速度慢
plan("sequential") #把并行计算改为单核计算

• 1.5 降维，可视化

### redunction
scRNA <- RunPCA(scRNA.sct.int, npcs = 50, verbose = FALSE)
ElbowPlot(scRNA, ndims=50)
pc.num=1:20
scRNA <- scRNA %>% RunTSNE(dims=pc.num) %>% RunUMAP(dims=pc.num)

### Visual
p <- DimPlot(scRNA, group.by = "orig.ident")
ggsave("UMAP_Samples_integr.pdf", p, width = 8, height = 6)
p <- DimPlot(scRNA, group.by = "orig.ident", split.by = "orig.ident", ncol = 4)
ggsave("UMAP_Samples_Split_integr.pdf", p, width = 18, height = 12)

##save seurat object
saveRDS(scRNA, "scRNA_SCT_int.rds")   

2. harmony整合

Harmony整合的官网教程及其原理此前已经介绍过：Harmony

• 2.1 准备数据

rm(list = ls())

scRNA <- readRDS("scRNA.rds")
cellinfo <- subset(scRNA@meta.data, select = c("orig.ident", "percent.mt", "percent.rb", "percent.HB"))
scRNA <- CreateSeuratObject(scRNA@assays$RNA@counts, meta.data = cellinfo)

• 2.2 数据标准化（和锚点整合不同，不需拆分样本，直接标准化）

### SCT标准化数据
scRNA <- SCTransform(scRNA)

• 2.3 使用harmony整合数据

### PCA
scRNA <- RunPCA(scRNA, npcs=50, verbose=FALSE)

### 整合方法1：单个样本间进行整合（推荐，效果更好）
scRNA <- RunHarmony(scRNA, group.by.vars="orig.ident", assay.use="SCT", max.iter.harmony = 20) 
# group.by.vars参数是设置按哪个分组来整合
# max.iter.harmony设置迭代次数，默认是10。运行RunHarmony结果会提示在迭代多少次后完成了收敛。
#⚠️RunHarmony函数中有个lambda参数，默认值是1，决定了Harmony整合的力度。lambda值调小，整合力度变大，反之。（只有这个参数影响整合力度，调整范围一般在0.5-2之间）

###整合方法2：按其他分类如不同分组来校正（实测效果不如方法1）
if(F){
  scRNA$batches <- scRNA$orig.ident
  scRNA$batches <- recode(scRNA$batches, 
                          "HNC01PBMC" = "batch1", 
                          "HNC01TIL" = "batch2",
                          "HNC10PBMC" = "batch1",
                          "HNC10TIL" = "batch2",
                          "HNC20PBMC" = "batch1",
                          "HNC20TIL" = "batch2",
                          "PBMC1" = "batch1",
                          "PBMC2" = "batch1",
                          "Tonsil1" = "batch3",
                          "Tonsil2" = "batch3")
  scRNA2 <- RunHarmony(scRNA, group.by.vars="batches", assay.use="SCT")
}

• 2.4 降维聚类

ElbowPlot(scRNA, ndims = 50)
pc.num=1:30
scRNA <- RunTSNE(scRNA, reduction="harmony", dims=pc.num) %>% RunUMAP(reduction="harmony", dims=pc.num)
#scRNA2 <- RunTSNE(scRNA2, reduction="harmony", dims=pc.num) %>% RunUMAP(reduction="harmony", dims=pc.num)

p <- DimPlot(scRNA, group.by = "orig.ident")
ggsave("UMAP_Samples_harmony.pdf", p, width = 8, height = 6)
p <- DimPlot(scRNA, group.by = "orig.ident", split.by = "orig.ident", ncol = 4)
ggsave("UMAP_Samples_Split_harmony.pdf", p, width = 18, height = 12)

##save seurat object
saveRDS(scRNA, "scRNA_SCT_harmony.rds") 

3. 整合结果评估

library(SingleR)
source("sc_function.R")

• 3.1 Seurat锚点整合结果评估

scRNA <- readRDS("scRNA_SCT_int.rds")
scRNA <- FindNeighbors(scRNA, dims = 1:20) %>% FindClusters()
load("ref_Hematopoietic.RData")
DefaultAssay(scRNA) <- "RNA"
scRNA <- cell_identify(scRNA, ref_Hematopoietic) #cell_identify是自己写的函数，ref_Hematopoietic是SingleR中的参考数据集
p <- DimPlot(scRNA, group.by = "SingleR", label = T)
ggsave("SingleR_Seurat.pdf", p, width = 8, height = 6)

DefaultAssay(scRNA) <- "integrated"
p <- FeaturePlot(scRNA, features = c('CD3E', 'CD3G', 'CD79B', 'MS4A1', 'GNLY', 'NKG7', 
                                     'CD14', 'FCGR3A', 'CD68', 'S100A12','CD1C', 'CST3', 
                                     'FCER1A', 'GZMB', 'IL3RA', 'PPBP'), ncol = 4)
ggsave("Features_Seurat_int.pdf", p, width = 18, height = 16)

这些基因是各种免疫细胞的marker基因，用的是整合之后的表达值(scRNA@assays$integrated@data)，可以这个矩阵对表达值的校正有点过，因此$integrated中的数据不建议拿来做后续分析。

DefaultAssay(scRNA) <- "SCT"
p <- FeaturePlot(scRNA, features = c('CD3E', 'CD3G', 'CD79B', 'MS4A1', 'GNLY', 'NKG7', 
                                     'CD14', 'FCGR3A', 'CD68', 'S100A12','CD1C', 'CST3', 
                                     'FCER1A', 'GZMB', 'IL3RA', 'PPBP'), ncol = 4)
ggsave("Features_Seurat_sct.pdf", p, width = 18, height = 16)

用data的数据绘图显示一些细胞群同时有多种特征性marker，说明可能存在着过度整合。

• 3.2 harmony整合结果评估

scRNA <- readRDS("scRNA_SCT_harmony.rds")
scRNA <- FindNeighbors(scRNA, dims = 1:30) %>% FindClusters()
load("ref_Hematopoietic.RData")
scRNA <- cell_identify(scRNA, ref_Hematopoietic)

p <- DimPlot(scRNA, group.by = "SingleR", label = T)
ggsave("SingleR_Harmony.pdf", p, width = 8, height = 6)

p <- FeaturePlot(scRNA, features = c('CD3E', 'CD3G', 'CD79B', 'MS4A1', 'GNLY', 'NKG7', 
                                     'CD14', 'FCGR3A', 'CD68', 'S100A12','CD1C', 'CST3', 
                                     'FCER1A', 'GZMB', 'IL3RA', 'PPBP'), ncol = 4)
ggsave("Features_Harmony.pdf", p, width = 18, height = 16)

和上一张图对比可以看出来Harmony整合的效果比较好，细胞群的marker基因分的也比较开。

saveRDS(scRNA, "scRNA.classified.rds")

4. 最后的吐槽

锚点整合是真的很容易过整合啊
同样的数据，左边是Harmony整合，右边是锚点整合。锚点整合完Ccr2+的细胞群都没了呢🤷♀️