单细胞转录组高级分析

单细胞分析内容本文主要学习生信会客厅公众号的系列文章为后续单细胞分析提供思路，学习笔记记录于此，仅做个人学习使用。

高级分析一：多样本合并与批次校正

scRNA 数据目前校正批次效应的算法有很多：MNN, CCA+MNN, Harmony, Scanorama, scMerge 等。

在 Seurat 中，将 CCA 与 MNN 算法结合起来，并参考 SNN 算法的理念设计了“锚点”评分体系，不仅可以校正实验的批次效应，还能跨平台整合数据，例如将 10x 单细胞数据、BD 单细胞数据和 SMART 单细胞数据整合在一起；也能整合单细胞多组学数据，例如将单细胞 ATAC、空间转录组与单细胞转录组数据整合在一起。

Seurat 整合流程与原理

使用 CCA 分析将两个数据集降维到同一个低维空间，因为 CCA 降维之后的空间距离不是相似性而是相关性，所以相同类型与状态的细胞可以克服技术偏倚重叠在一起。
*使用 PCA 降维，细胞之间的距离体现的是转录特征相似性，批次效应引入的系统误差会使样本分离。而使用 CCA 降维，细胞之间的距离体现的是转录特征相关性，因此同类型且同状态的细胞可以跨越技术差异重叠在一起。
CCA 降维之后细胞在低维空间有了可以度量的“距离”，MNN(mutual nearest neighbor)算法以此找到两个数据集之间互相“距离”最近的细胞，Seurat 将这些相互最近邻细胞称为“锚点细胞”。
锚点过滤，理想情况下相同类型和状态的细胞才能构成配对锚点细胞，但是异常的情况也会出现。在 query 数据中存在的细胞类型在 reference 数据集没有相同类型的细胞，但是它也找到了锚点配对细胞。所以需要过滤
- 在 CCA 低维空间找到的锚点，返回到基因表达数据构建的高维空间中验证，如果它们的转录特征相似性高则保留，否则过滤此锚点。
- 检查锚点细胞所在数据集最邻近的 30 个细胞，查看它们重叠的锚点配对细胞的数量，重叠越多分值越高，代表锚点可靠性更高。
经过层层过滤剩下的锚点细胞对，可以认为它们是相同类型和状态的细胞，它们之间的基因表达差异是技术偏倚引起的。Seurat 计算它们的差异向量，然后用此向量校正这个锚点锚定的细胞子集的基因表达值。校正后的基因表达值即消除了技术偏倚，实现了两个单细胞数据集的整合。

高级分析二：转录调控网络分析

组织内细胞异质性的基础是细胞转录状态的差异，转录状态的特异性又是由转录因子主导的基因调控网络（GRNs）决定并维持稳定的。然而单细胞转录组数据具有背景噪音高、基因检出率低和表达矩阵稀疏性的特点，给传统统计学和生物信息学方法推断高质量的 GRNs 带来了挑战。
SCENIC 是一种专为单细胞数据开发的 GRNs 算法,创新之处在于引入了转录因子 motif 序列验证统计学方法推断的基因共表达网络，从而识别高可靠性的由转录因子主导的 GRNs。

SCENIC 分析流程

主要分析有四步：

GENIE3/GRNBoost：基于共表达情况鉴定每个 TF 的潜在靶点；
RcisTarget：基于 DNA-motif 分析选择潜在的直接结合靶点；
AUCell：分析每个细胞的 regulons 活性；
细胞聚类：基于 regulons 的活性鉴定稳定的细胞状态并对结果进行探索

SCENIC-network

# SCENIC分析需要输入单细胞表达矩阵和细胞meta信息
# 还要设置运行线程数，参考数据库目录、版本，细胞meta信息存放目录等配置信息。
library(Seurat)
library(tidyverse)
library(patchwork)
library(SCENIC)
rm(list=ls())
##==分析准备==##
dir.create("SCENIC")
dir.create("SCENIC/int")
scRNA <- readRDS("scRNA.rds")
setwd("~/project/10xDemo2/SCENIC")
##准备细胞meta信息
cellInfo <- data.frame(scRNA@meta.data)
colnames(cellInfo)[which(colnames(cellInfo)=="orig.ident")] <- "sample"
colnames(cellInfo)[which(colnames(cellInfo)=="seurat_clusters")] <- "cluster"
colnames(cellInfo)[which(colnames(cellInfo)=="celltype_Monaco")] <- "celltype"
cellInfo <- cellInfo[,c("sample","cluster","celltype")]
saveRDS(cellInfo, file="int/cellInfo.Rds")
##准备表达矩阵
#为了节省计算资源，随机抽取1000个细胞的数据子集
subcell <- sample(colnames(scRNA),1000)
scRNAsub <- scRNA[,subcell]
saveRDS(scRNAsub, "scRNAsub.rds")
exprMat <- as.matrix(scRNAsub@assays$RNA@counts)
##设置分析环境
mydbDIR <- "./cisTarget"
mydbs <- c("hg38__refseq-r80__500bp_up_and_100bp_down_tss.mc9nr.feather",
           "hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.feather")
names(mydbs) <- c("500bp", "10kb")
scenicOptions <- initializeScenic(org="hgnc",
                                  nCores=8,    # 线程数
                                  dbDir=mydbDIR,   # 设置数据库存放目录
                                  dbs = mydbs,          # 数据库文件
                                  datasetTitle = "HNSCC")
saveRDS(scenicOptions, "int/scenicOptions.rds")

# 第一步：计算转录因子与每个基因的相关性
# 当处理大型数据时，可选择采用GENIE3推断共表达模块时随机抽取少量细胞计算，计算regulons的活性时，所有细胞都代入运算；或者使用python版本的SCENIC
##基因过滤
#过滤标准是基因表达量之和>细胞数*3%，且在1%的细胞中表达
genesKept <- geneFiltering(exprMat, scenicOptions,
              minCountsPerGene = 3 * 0.01 * ncol(exprMat),
              minSamples = ncol(exprMat) * 0.01)
exprMat_filtered <- exprMat[genesKept, ]
##计算相关性矩阵
runCorrelation(exprMat_filtered, scenicOptions)
##TF-Targets相关性回归分析
exprMat_filtered_log <- log2(exprMat_filtered+1)
runGenie3(exprMat_filtered_log, scenicOptions, nParts = 20)  # nParts作用是把表达矩阵分成n份分开计算，目的是防止数据量大时内存不够。
#这一步消耗的计算资源非常大，个人电脑需要几个小时的运行时间
# 生成多个中间文件
# 1.2_corrMat.Rds：基因之间的相关性矩阵
# 1.3_GENIE3_weightMatrix_part_1.Rds等：GENIE3的中间结果
# 1.4_GENIE3_linkList.Rds：GENIE3最终结果，是把“1.3_”开头的文件合并在一起。

# 第二步：推断共表达模块
# 过滤低相关性的组合形成共表达基因集（模块）
# 存在6种过滤标准
# w001：以每个TF为核心保留weight>0.001的基因形成共表达模块；
# w005：以每个TF为核心保留weight>0.005的基因形成共表达模块；
# top50：以每个TF为核心保留weight值top50的基因形成共表达模块；
# top5perTarget：每个基因保留weight值top5的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
# top10perTarget：每个基因保留weight值top10的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；
# top50perTarget：每个基因保留weight值top50的TF得到精简的TF-Target关联表，然后把基因分配给TF构建共表达模块；

runSCENIC_1_coexNetwork2modules(scenicOptions)

# 结果文件中corr是runCorrelation(exprMat_filtered, scenicOptions)命令得到的，1代表激活，-1代表抑制，0代表中性，SCENIC只会采用corr值为1的数据用于后续分析。

# 第三步：Motif验证共表达模块
# 修剪共表达模块形成有生物学意义的调控单元（regulons）
# 1. 对每个共表达模块进行motif富集分析，保留显著富集的motif；此项分析依赖gene-motif评分（排行）数据库，其行为基因、列为motif、值为排名，就是下载的cisTarget数据库。 https://resources.aertslab.org/cistarget/
# 2. 使用数据库对motif进行TF注释，注释结果分高、低可信度 。数据库直接注释和同源基因推断的TF是高可信结果，使用motif序列相似性注释的TF是低可信结果。
# 3. 用保留的motif对共表达模块内的基因进行打分（同样依据cisTarget数据库），识别显著高分的基因（理解为motif离这些基因的TSS很近）；
# 4. 删除共表达模块内与motif评分不高的基因，剩下的基因集作者称为调控单元（regulon）。

##推断转录调控网络（regulon）

runSCENIC_2_createRegulons(scenicOptions)

#以上代码可增加参数coexMethod=c("w001", "w005", "top50", "top5perTarget", "top10perTarget", "top50perTarget"))
#默认6种方法的共表达网络都计算，可以少选几种方法以减少计算量

# 第四步：Regulon活性评分与可视化
# 评分的基础是基因的表达值，分数越高代表基因集的激活程度越高。
##regulons计算AUC值并进行下游分析
exprMat_all <- as.matrix(scRNA@assays$RNA@counts)
exprMat_all <- log2(exprMat_all+1)
runSCENIC_3_scoreCells(scenicOptions, exprMat=exprMat_all)

# a. 计算regulon在每个细胞中AUC值，最后得到一个以regulon为行细胞为列的矩阵。
结果目录：int/3.4_regulonAUC.Rds
# b. 计算每个regulon的AUC阈值，细胞中regulonAUC值>阈值，代表此regulon在细胞中处于激活状态，否则代表沉默状态。结果目录：int/3.5_AUCellThresholds.Rds
# c. 使用regulonAUC矩阵对细胞进行降维聚类
# d. 用heatmap图展示regulonAUC矩阵，用t-SNE图分别展示每个regulon的活性分布情况。结果目录：output/Step3_开头的系列文件

# Regulon活性二进制转换与可视化
# 对于细胞类型清晰的数据集而言，构建regulons并对其活性打分足够后续分析。但是，在很多情况下将评分转换为二进制的“开/关”，则既方便解释，又能最大化体现细胞类型的差异。将特定的regulon转换为“0/1”有利于探索和解释关键转录因子。将所有的regulons转换为“0/1”后创建二进制的活性矩阵，则可以用于细胞聚类，对消除技术偏倚特别有用。AUCell会自动计算可能的阈值进行“0/1”转换，作者建议在转换之前手工检查并调整这些阈值

#使用shiny互动调整阈值
aucellApp <- plotTsne_AUCellApp(scenicOptions, exprMat_all)
savedSelections <- shiny::runApp(aucellApp)
#保存调整后的阈值
newThresholds <- savedSelections$thresholds
scenicOptions@fileNames$int["aucell_thresholds",1] <- "int/newThresholds.Rds"
saveRDS(newThresholds, file=getIntName(scenicOptions, "aucell_thresholds"))
saveRDS(scenicOptions, file="int/scenicOptions.Rds")

# 二进制转换及衍生分析
# 将regulonAUC矩阵转换为二进制矩阵后，会重新降维聚类
runSCENIC_4_aucell_binarize(scenicOptions, exprMat=exprMat_all)

# 结果可视化
# 调用SCENIC的分析结果，使用seurat和pheatmap进行可视化

# Seurat可视化SCENIC结果
# 把SCENIC结果中最重要的regulonAUC矩阵导入Seurat
##导入原始regulonAUC矩阵进行初步处理
AUCmatrix <- readRDS("int/3.4_regulonAUC.Rds")
AUCmatrix <- AUCmatrix@assays@data@listData$AUC
AUCmatrix <- data.frame(t(AUCmatrix), check.names=F)
RegulonName_AUC <- colnames(AUCmatrix)
RegulonName_AUC <- gsub(' \\(','_',RegulonName_AUC)
RegulonName_AUC <- gsub('\\)','',RegulonName_AUC)
colnames(AUCmatrix) <- RegulonName_AUC
scRNAauc <- AddMetaData(scRNA, AUCmatrix)
scRNAauc@assays$integrated <- NULL
saveRDS(scRNAauc,'scRNAauc.rds')

##导入二进制regulonAUC矩阵进行初步处理
BINmatrix <- readRDS("int/4.1_binaryRegulonActivity.Rds")
BINmatrix <- data.frame(t(BINmatrix), check.names=F)
RegulonName_BIN <- colnames(BINmatrix)
RegulonName_BIN <- gsub(' \\(','_',RegulonName_BIN)
RegulonName_BIN <- gsub('\\)','',RegulonName_BIN)
colnames(BINmatrix) <- RegulonName_BIN
scRNAbin <- AddMetaData(scRNA, BINmatrix)
scRNAbin@assays$integrated <- NULL
saveRDS(scRNAbin, 'scRNAbin.rds')

## 利用Seurat可视化AUC
## SCENIC对象与Seurat对象一致
dir.create('scenic_seurat')
#FeaturePlot
p1 = FeaturePlot(scRNAauc, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p2 = FeaturePlot(scRNAbin, features='CEBPB_extended_2290g', label=T, reduction = 'tsne')
p3 = DimPlot(scRNA, reduction = 'tsne', group.by = "celltype_Monaco", label=T)
plotc = p1|p2|p3
ggsave('scenic_seurat/CEBPB_extended_2290g.png', plotc, width=14 ,height=4)

#RidgePlot&VlnPlot  山脊图及小提琴图
p1 = RidgePlot(scRNAauc, features = "CEBPB_extended_2290g", group.by="celltype_Monaco") +
               theme(legend.position='none')
p2 = VlnPlot(scRNAauc, features = "CEBPB_extended_2290g", pt.size = 0, group.by="celltype_Monaco") +
             theme(legend.position='none')
plotc = p1 + p2
ggsave('scenic_seurat/Ridge-Vln_CEBPB_extended_2290g.png', plotc, width=10, height=8)


# pheatmap可视化SCENIC结果
library(pheatmap)
cellInfo <- readRDS("int/cellInfo.Rds")
celltype = subset(cellInfo,select = 'celltype')
AUCmatrix <- t(AUCmatrix)
BINmatrix <- t(BINmatrix)
#挑选部分感兴趣的regulons
my.regulons <- c('ETS1_2372g','ETV7_981g','IRF7_239g','XBP1_854g','ATF4_37g',
                 'KLF13_78g','ATF6_129g','CREB3L2_619g','TAGLN2_13g',
                 'STAT1_extended_1808g','CEBPB_extended_2290g','IRF5_extended_422g',
                 'SPI1_1606g','HMGA1_14g','SPIB_1866g','IRF8_348g','BCL11A_136g',
                 'EBF1_40g','MAF_45g','BATF_131g','FOXP3_55g','TBX21_388g',
                 'EOMES_extended_101g','TCF7_extended_31g','LEF1_extended_49g')
myAUCmatrix <- AUCmatrix[rownames(AUCmatrix)%in%my.regulons,]
myBINmatrix <- BINmatrix[rownames(BINmatrix)%in%my.regulons,]
#使用regulon原始AUC值绘制热图
pheatmap(myAUCmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myAUCmatrix_heatmap.png',
         width = 6, height = 5)
#使用regulon二进制AUC值绘制热图
pheatmap(myBINmatrix, show_colnames=F, annotation_col=celltype,
         filename = 'scenic_seurat/myBINmatrix_heatmap.png',
         color = colorRampPalette(colors = c("white","black"))(100),
         width = 6, height = 5)

高级分析三：细胞通讯分析

细胞通讯研究领域涵盖的内容很广，如上图所示包括通讯方式、功能、信号分子以及各种途径的机制。细胞之间通讯的介质有很多，例如钙离子、脂质、多肽、蛋白、外泌体以及电信号等。利用单细胞转录组数据分析的细胞通讯，仅限于蛋白质配体-受体复合物介导的细胞间通讯。**其分析的基础是基因表达数据和配体-受体数据库信息**，例如转录组数据表明 A、B 细胞分别表达了基因 α 和 β，通过数据库查询 α 和 β 是配体-受体关系，则认为 A-B 通过 α-β 途径进行了通讯。
CellPhoneDB(网页版，Python 版)：提供细胞通讯配体-受体数据库
iTALK(可视化效果最好，R)：集成了多种差异分析和可视化方法；将配体-受体注释为 4 大类：细胞因子、生长因子、免疫检查点和其他；配体-受体分析的中间数据和最终结果都可以轻松导出。
celltalker(R)：更倾向于分析样本之间具有差异的细胞通讯关系，认定细胞进行通讯的前提是配体和受体的表达值在通讯的细胞之间具有一致性。

高级分析四：scRNA 数据推断 CNV

inferCNV流程

原始数据进行标准化处理，去除测序深度不同造成的差异，结果是流程图中的第 1 个热图；
将基因表达值转换为 CNV 值，结果是流程图中的第 2 个热图；
对每个细胞的 CNV 值进行中心化处理，使细胞之间相同染色体区域的 CNV 值具有可比性，结果是流程图中的第 3 个热图；
用肿瘤细胞的 CNV 值减去对应区域正常细胞的 CNV 值，使热图展现的 CNV 结果更直观，结果是流程图中的第 4 个热图；
如果 infercnv::run 函数中的参数 denoise=TRUE，则使用算法进一步去除背景噪音凸显 CNV 区域，结果是流程图中的蓝框左图；
如果 infercnv::run 函数中的参数 HMM=TRUE，则使用隐马尔可夫模型（Hidden Markov Model, HMM）预测 CNV 区域，并用贝叶斯潜在混合模型（Bayesian Network Latent Mixture Model）对结果进行校正，结果是流程图中的蓝框下图。

inferCNV 运行需要三个输入文件，原始表达矩阵、细胞注释文件和基因定位文件
其中原始表达矩阵：行为基因，列为细胞；细胞注释文件即行为细胞名，第二列为细胞类型；基因定位文件行为基因名，第二列为所在染色体，第三列第四列为起始位置。

#创建inferCNV对象
infercnv_obj = CreateInfercnvObject(delim = '\t',
                  raw_counts_matrix = 'inferCNV/exprMatrix.txt',
                  annotations_file = 'inferCNV/cellAnnota.txt',
                  gene_order_file = 'inferCNV/geneLocate.txt',
                  ref_group_names = c("B_cell","Monocyte","NK_cell","T_cells"))
dir.create("inferCNV/gse149180")
#10x数据cutoff推荐使用0.1
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=0.1,
                             out_dir='inferCNV/gse149180/',
                             cluster_by_groups=TRUE,
                             denoise=TRUE,
                             HMM=TRUE)
# 输出文件在out_dir目录下
# 分为中间过程数据、初步分析结果、去噪分析结果、HMM预测后结果、最终分析结果5部分
# infercnv开头的为分析结果
# infercnv.png : 去噪之后的最终热图
# infercnv.references.txt : 正常细胞的CNV分值矩阵
# infercnv.observations.txt : 肿瘤细胞的CNV分值矩阵
# infercnv.observation_groupings.txt : 肿瘤细胞的聚类分组
# infercnv.observations_dendrogram.txt : 热图的newick格式文件

其中，HMM 模型结果：

0：完全缺失数的变异
0.5：缺失一个拷贝数的变异
1：正常
1.5：增加一个拷贝数的变异
2：增加两个拷贝数的变异
3：所有大于两个拷贝数的变异

高级分析五：GSEA 与 GSVA 分析

GSEA：Gene Set Enrichment Analysis
GSVA：Gene Set Variation Analysis
以基因集为单位研究数据。基于基因集的富集矩阵，我们可以做差异分析、表型相关性分析、生存分析，也可以基于基因集的功能来鉴定细胞类型。

高级分析六：整合 scATAC 数据

# scATAC数据预处理
library(Seurat)
library(ggplot2)
library(rtracklayer)
library(patchwork)
rm(list=ls())
dir.create("ATAC")
set.seed(811)
peaks <- Read10X_h5("Data/ATAC/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
# 初始文件是peak矩阵，行为染色体坐标，列为细胞barcode，即代表染色质可及性区域

# 假设基因上游0~2kb区域的reads数可以代表基因的活性，并把上面的peak矩阵转换为基因活性矩阵
activity.matrix <- CreateGeneActivityMatrix(peak.matrix = peaks,
                   annotation.file = "Data/ATAC/Homo_sapiens.GRCh37.82.gtf",
                   seq.levels = c(1:22, "X", "Y"), upstream = 2000, verbose = TRUE)

# 使用peak矩阵、活性矩阵和metadata创建seurat对象
pbmc.atac <- CreateSeuratObject(counts = peaks, assay = "ATAC", project = "10x_ATAC")
pbmc.atac[["ACTIVITY"]] <- CreateAssayObject(counts = activity.matrix)
meta <- read.table("Data/ATAC/atac_v1_pbmc_10k_singlecell.csv", sep = ",",
                   header = TRUE, row.names = 1, stringsAsFactors = FALSE)
meta <- meta[colnames(pbmc.atac), ]
pbmc.atac <- AddMetaData(pbmc.atac, metadata = meta)
pbmc.atac <- subset(pbmc.atac, subset = nCount_ATAC > 5000)
pbmc.atac$tech <- "atac"

# 使用peak矩阵对scATAC数据降维，并对活性矩阵进行标准化
#对活性矩阵执行标准化与中心化处理
DefaultAssay(pbmc.atac) <- "ACTIVITY"
pbmc.atac <- FindVariableFeatures(pbmc.atac)
pbmc.atac <- NormalizeData(pbmc.atac)
pbmc.atac <- ScaleData(pbmc.atac)
# 对peak矩阵降维
DefaultAssay(pbmc.atac) <- "ATAC"
VariableFeatures(pbmc.atac) <- names(which(Matrix::rowSums(pbmc.atac) > 100))
pbmc.atac <- RunLSI(pbmc.atac, n = 50, scale.max = NULL)
pbmc.atac <- RunUMAP(pbmc.atac, reduction = "lsi", dims = 1:50)

# 整合前scRNA与scATAC数据展示
pbmc.rna <- readRDS("Data/ATAC/pbmc_10k_v3.rds")
pbmc.rna$tech <- "rna"
#展示整合前的scATAC和scRNA数据
p1 <- DimPlot(pbmc.atac, reduction = "umap") + NoLegend() + ggtitle("scATAC-seq")
p2 <- DimPlot(pbmc.rna, group.by="celltype", label=TRUE, label.size=2,
              repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
plotc = p1 + p2
ggsave('ATAC/ATAC-RNA_beforeINT.png', plotc, width = 8, height = 4)

# 识别scATAC数据集的细胞类型(通过参考数据集进行)
#通过锚点建立数据集之间的配对关系
transfer.anchors <- FindTransferAnchors(reference=pbmc.rna, query = pbmc.atac,
                    features=VariableFeatures(object=pbmc.rna), reference.assay = "RNA",
                    query.assay = "ACTIVITY", reduction = "cca")
#数据迁移
celltype.predictions <- TransferData(anchorset = transfer.anchors, refdata = pbmc.rna$celltype,
                        weight.reduction = pbmc.atac[["lsi"]])
pbmc.atac <- AddMetaData(pbmc.atac, metadata = celltype.predictions)
#展示id转换后的scATAC和scRNA数据
pbmc.atac.filtered <- subset(pbmc.atac, subset = prediction.score.max > 0.5)
pbmc.atac.filtered$predicted.id <- factor(pbmc.atac.filtered$predicted.id, levels=levels(pbmc.rna))
p1 <- DimPlot(pbmc.atac.filtered, group.by="predicted.id", label=TRUE, label.size=2, repel=TRUE) +
              ggtitle("scATAC-seq cells") + NoLegend() + scale_colour_hue(drop = FALSE)
p2 <- DimPlot(pbmc.rna, group.by="celltype", label=TRUE, label.size=2, repel=TRUE) +
              ggtitle("scRNA-seq cells") + NoLegend()
plotc = p1 + p2
ggsave('ATAC/ATAC-RNA_transfer.png', plotc, width = 8, height = 4)

# 将scRNA细胞的转录谱通过锚点转移给scATAC细胞
genes.use <- VariableFeatures(pbmc.rna)
refdata <- GetAssayData(pbmc.rna, assay="RNA", slot="data")[genes.use, ]
imputation <- TransferData(anchorset=transfer.anchors, refdata=refdata, weight.reduction = pbmc.atac[["lsi"]])
pbmc.atac[["RNA"]] <- imputation
#合并两个数据集的seurat对象
coembed <- merge(x = pbmc.rna, y = pbmc.atac)
#合并后的数据集执行降维分析
coembed <- ScaleData(coembed, features = genes.use, do.scale = FALSE)
coembed <- RunPCA(coembed, features = genes.use, verbose = FALSE)
coembed <- RunUMAP(coembed, dims = 1:30)
coembed$celltype <- ifelse(!is.na(coembed$celltype), coembed$celltype, coembed$predicted.id)
p1 <- DimPlot(coembed, group.by = "tech")
p2 <- DimPlot(coembed, group.by="celltype", label=TRUE, label.size=2, repel=TRUE) + NoLegend()
plotc = p1 + p2
ggsave('ATAC/ATAC-RNA_coembed.png', plotc, width = 8, height = 4)

单细胞转录组还能分析些什么？