单细胞测序（一）

本期讲述一下单细胞测序的基础操作

使用的数据是GEO数据库的GSE114417：GEO Accession viewer

整体设计：4个样本，两组：DAPI- sorted cells（DAPISCs）和 DAPI-CD45-CD31-EpCAM-PDPN+ sorted cells (CAFs)各两个样本

下载数据后解压

按名称整理到四个文件夹中，并进行重命名：

把每一个文件夹中的压缩包重命名成相同的格式，便于读取：

下面进行代码分析：

library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)
library(harmony)
library(cowplot)

###设置路径进行读取###

dir=c("CAFs.1/","CAFs.2/","DAPISCs.1/","DAPISCs.2/")
names(dir) = c('CAF1', 'CAF2', 'DapiNeg1',"DapiNeg2")

scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])#读取每一个样本的单细胞数据
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)#建立Seurat对象，设置的条件为每个基因至少表达于3个细胞，每个细胞至少表达300个基因，具体情况具体调整标准
}

for (i in 1:length(scRNAlist)) {
  scRNAlist[[i]] <- NormalizeData(scRNAlist[[i]])#进行归一化处理，移除测序深度的影响
  scRNAlist[[i]] <- FindVariableFeatures(scRNAlist[[i]], selection.method = "vst",nfeatures = 3000)#寻找高边基因，简单说就是找到一些在各个Cluster中变化幅度较大的基因，此处设置为3000个
}

#此处可以进行多线程运算用来加快速度，考虑到数据量不大，便没有进行

#关于SCTransform的评估与教程后续也会有所展示

scRNA.anchors <- FindIntegrationAnchors(object.list = scRNAlist,anchor.features = 3000)#此处使用锚定点算法来去除批次效应，也可以使用Harmony算法
scRNA1 <- IntegrateData(anchorset = scRNA.anchors)#整合得到新的Seurat对象

DefaultAssay(scRNA1) <- "integrated"#将默认处理的数据设定为整合后的对象
scRNA1 <-  ScaleData(scRNA1)#在进行PCA降维聚类之前要进行归一化处理
scRNA1 <- RunPCA(scRNA1, npcs = 20, verbose = T)#PCA聚类，此处设置的主成分为20，具体确定应根据“肘部法则”进行确定，保证结果更加准确高效

scRNA1 <- FindNeighbors(scRNA1, reduction = "pca", dims = 1:7)
scRNA1 <- FindClusters(scRNA1, resolution = 0.4)#进行谱聚类，基于共享最近邻图（SNN）和模块化优化来识别细胞Cluster
scRNA1 <- RunUMAP(scRNA1, reduction = "pca", dims = 1:7)#进行UMAP
scRNA1 <- RunTSNE(scRNA1, dims = 1:7)#进行tSNE

DefaultAssay(scRNA1) <- "RNA"

#此处进行简单讲解，tSNE是非线性降维，只是为了可视化而不是特征提取，UMAP点与点之间的距离可以说明Cluster的关系，而tSNE无法说明，同时UMAP相比于tSNE的优点：①尽可能多保留全局结构②耗时更短③可扩展到更大维度数据集

#下面进行对细胞亚群的注释，使用的是SingleR包

library(SingleR)

load("ref_Mouse_all.RData")#关于注释文件可以去网上搜索下载，有不少资源针对于各个物种的

refdata <- ref_Mouse

scRNA1 <- JoinLayers(scRNA1)#Seurat对象更新到V5后，新增了layers层，需要简单合并一下
testdata <- GetAssayData(scRNA1, slot="data")#提取RNA表达矩阵

clusters <- scRNA1@meta.data$seurat_clusters#提取Clusters

cellpred <- SingleR(test = testdata, ref = refdata, labels = refdata$label.main,
                    method = "cluster", clusters = clusters,
                    assay.type.test = "logcounts", assay.type.ref = "logcounts")#使用SingleR进行注释

celltype = data.frame(ClusterID=rownames(cellpred), celltype=cellpred$labels, stringsAsFactors = FALSE)#得到细胞Cluster的注释文件
scRNA1@meta.data$celltype = "NA"#在meta.data中新增一列
for(i in 1:nrow(celltype)){
  scRNA1@meta.data[which(scRNA1@meta.data$seurat_clusters == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}#将细胞类型写入Seurat对象中

DimPlot(scRNA1, reduction = "umap", group.by = "celltype",label = T)

DimPlot(scRNA1, reduction = "tsne", group.by = "celltype",label = T)