单细胞数据挖掘

目录

1. 实验原理
   1. 数据集
   2. 质控
   3. 挑选高变异基因
   4. 降维与聚类
   5. 细胞类型注释
   6. 轨迹分析
2. R语言实现
   1. 提取数据
   2. 质控
   3. 挑选高变异基因
   4. 降维与聚类，主要是找marker基因
   5. QC之二
   6. 细胞类型注释
   7. 轨迹分析
3. 尾声
   1. 多样本整合
   2. infeCNV
   3. copyKAT
   4. 细胞通讯

案例里用了《Glioblastoma cell differentiation trajectory predicts the immunotherapy response and overall survival of patients》，主要是使用integrated scRNA-seq（综合（？）单细胞RNA测序）和bulk RNA-seq（整体（？）RNA测序）来找出GBM分化相关基因然后预测患者总生存期之类的。

换一个类似的做实验，也不懂各种数据有什么差别，就选《Cell differentiation trajectory predicts patient potential immunotherapy response and prognosis in gastric cancer》这篇文章，主要是使用integrated scRNA-seq和bulk RNA-seq来找出GC分化相关基因然后预测患者总生存期之类的。（干的事都差不多，真是原原又神神啊）
这个用的数据集是GSE112302和GSE84437，和案例类似的部分用的是GSE112302。

实验原理

主要就是干了质控、基因可视化、降维与聚类、周期判断、double droplet、细胞类型注释、拟时序分析。

数据集

这堆数据集就是一轴是细胞，一轴是基因种类，每个点代表某个细胞/样本的某个基因表达量。

质控

质控主要就是过滤：

• 低质量的细胞(比如基因少、线粒体基因多)
• 双细胞（两个细胞被捕获到同一个液滴中）
• 高ERCC基因（External RNA Controls Consortium），“是一组合成的RNA分子，它们的浓度已知，并且在测序实验中添加到样本中，以便可以用来校正技术变异和进行数据标准化”，在所有细胞中的浓度应该是一致的，所以如果ERCC基因百分比高，那这个细胞的RNA可能较少。

Seurat包可以傻瓜实现。

案例里还提到了double droplet和周期判断。
double droplet就是上面说的双细胞。
而周期判断就是说，一些细胞处于分裂期，它们会高表达某些基因，可能会影响聚类，所以要去掉。
scran包可以傻瓜实现。

挑选高变异基因

highly Variable gene是在一组样本或细胞中，表达水平变化大的基因，它们可能是驱动细胞异质性的关键。 Seurat包同样可以傻瓜实现。

降维与聚类

降维机器学习里也有，“将数据从原始的高维空间转换到低维空间，同时尽可能保留数据的变异性”，总的来说目的就是去掉冗杂信息。一个方法就是PCA（主成分分析），能使得降维后数据的方差尽可能大（保留变异性），而和原本数据的均方误差尽可能小(减小和原数据的差距)，这两个其实殊途同归。更多可见另一篇博客。
PCA对数据的尺度（可以用标准差度量）非常敏感，如果不同的特征有不同的尺度（0-1vs0-1000），那么PCA的结果将主要受到尺度大的特征的影响，所以PCA之前可以先缩放特征到同一尺度。

聚类就是把数据分成几类，使得同一个类内数据的相似性尽可能大，同时不同类数据的差异性也尽可能地大。聚类是一种无监督学习，也就是机器自己把数据根据特征分类，而不需要标注的指导。更多可见另一篇博客。

案例里PCA之后聚类然后又用t-SNE降维（也可以用UMAP），后者主要是为了降至2或3维方便可视化。

Seurat包还是可以傻瓜实现。that's right the square hole

细胞类型注释

字面意思，就比如把一些细胞标注成B细胞之类的甚至更精细一些，幼稚B细胞之类的。这靠的是一些特异性高的marker基因。 可以用SingleR包。

轨迹分析

就是“构建细胞的变化轨迹，重塑细胞随时间的变化过程”，目的是绘制祖细胞到分化细胞的进展情况。
可以用Monocle包。

---

本来这点东西一天下班时间就能搞完的，现在下来两三天了，不光是因为只有业余时间，还真是怠惰啊。

---

R语言实现

代码详见single-cell-data-mining。

提取数据

怎么两个数据集形式差那么多，你们能不能统一一下（骂骂咧咧）。

先看GSE84465学习一下。
案例里的GSE84465集合到了一个.csv里，列名是某样本的plate_id和Well，行名就是基因。
然后还有一个SraRunTable.txt，刚刚还奇怪，上面的文件都没有分组信息，这个文件里就有，就是每个样本的一些元信息，比如组织类型、患者编号、对应的样本编号等。这个就是SRA里的元数据，进SRA Run Selector里就能找到。
所以这个SraRunTable.txt就是用来给上面那些scRNA-seq数据分组。

我们要做的就是从SRA生成GSE对应的列名，并使用SRA的信息筛选出肿瘤组织，借助对应的列名筛选出GSE中的肿瘤组织。
还有SRA记录着分组信息，也要一起筛选出来。

所以说拿到GSE112302也是一样的，就是GSE112302没有给全部合在一起的数据，需要我们自己合起来（附加列，也就是样本），然后用SRA筛选出肿瘤组织。

质控

创建Seurat对象，简单筛选细胞

要先创建Seurat对象。CreateSeuratObject()函数的文档如下：

• counts：这是一个矩阵类对象，包含未经标准化的数据，其中细胞作为列，特征作为行，或者是一个Assay派生对象。
• assay：初始count矩阵的名称。
• meta.data：要添加到Seurat对象的额外细胞级元数据。应该是一个数据框，其中行是细胞名称，列是额外的元数据字段。元数据中的行名称需要与计数矩阵的列名称匹配。
• project：Seurat对象的项目名称。
• ...：传递给其他方法的参数。
• min.cells：至少在这么多细胞中检测到的特征。也会对计数矩阵进行取子集。要重新引入排除的特征，使用较低的截止值创建一个新对象。
• min.features：至少在这么多特征中检测到的细胞。

（对了R语言的函数是这样调用的：function_name(arg1 = value1, arg2 = value2, ...)）

其实只要用到counts、meta.data、min.cells、min.features就行了。 这里counts就是count矩阵，meta.data就是分组信息，用data.frame()上面说了是这样的：...是列（额外的元数据字段，这里就是病人编号，要用来分组嘛），row.names是行（样本编号）,要和counts的列名对应。
然后后两个前者用来筛选表达的细胞的数量低的基因，后者用来筛选基因少的细胞，顺便还能过滤空细胞。

筛选线粒体基因和ERCC基因过多细胞

然后还要筛选掉线粒体基因过多的细胞：

scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-") #先全部算了比较好，PercentageFeatureSet()计算满足给定模式的特征在每个细胞中的百分比，后面那个参数就是正则（真是先进啊（大嘘））
scRNA <- subset(scRNA, subset = percent.mt < 5) #案例取5%，subset()就是Array.filter()（确信）

还有ERCC基因一样的，选40%。

筛选双细胞

筛选双细胞应当在一个样本内进行，先一套下来到聚类，然后

sweep.res.list <- paramSweep(seurat_object, PCs = pc.num, sct = F)
sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <- find.pK(sweep.stats)
pK_bcmvn <- bcmvn$pK[which.max(bcmvn$BCmetric)] %>% as.character() %>% as.numeric()
homotypic.prop <- modelHomotypic(seurat_object$seurat_clusters)
nExp_poi <- round(doublet.rate*ncol(seurat_object))
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
seurat_object <- doubletFinder(seurat_object, PCs = pc.num, pN = 0.25, pK = pK_bcmvn, nExp = nExp_poi.adj, reuse.pANN = F, sct = F)

周期分析

这个要先标准化...别急...看seurat#7942。 我说怎么把这两个放在后面。

seurat_object <- CellCycleScoring(seurat_object, s.features= gene_s, g2m.features= gene_g2m)

gene_s和gene_g2m可以下载也可以自己提供，注意物种区别。

质控绘图

最后就是绘图。
看他们都是 - 图A：观察不同组cell的counts、feature分布 - 图B：nCount_RNA与对应的nFeature_RNA关系

图A用小提琴图，复习一下，VlnPlot()：

• object：Seurat对象。
• features：要绘制的特征。
• group.by：用于分组的元数据。
• pt.size：点的大小。
• cols：颜色。可以用rainbow(n)生成。

然后可以用+连接修改，theme(legend.position = "none")隐藏图例，labs(tag = "A")添加全局标签A。

图B用FeatureScatter()，绘制两个特征之间的关系。

• object：Seurat对象。
• feature1：x轴特征。
• feature2：y轴特征。
• group.by：用于分组的元数据。
• 剩下通用的都差不多。

案例还对ERCC基因比例和nCount用了FeatureScatter()，这个图一看也验证了上面说的，ERCC基因百分比高，nCount就低。

挑选高变异基因

用FindVariableFeatures()：

• object：Seurat对象。
• selection.method：选择方法。可以是"vst"（variance stabilizing transformation）或"mean.var.plot"（均值方差图）。
• nfeatures：选择为顶部变异特征的特征数量；只在selection.method设置为'dispersion'或'vst'时使用。
• 还有一些其它的，还有使用vst法时的参数。

然后用VariableFeaturePlot()绘制变异特征图。
LabelPoints()可以标记点，head()就可以取前几个。

降维与聚类，主要是找marker基因

PCA

先NormalizeData()，对每个细胞之间的某个基因表达量进行标准化，可以消除不同细胞之间由于测序深度不同而产生的差异，所以归一化后，不同的基因可能仍然有不同的表达量级别。
然后ScaleData()，对每个基因的表达量进行缩放并中心化，可以消除不同基因之间由于表达量级别不同而产生的差异，并且PCA需要基于中心化的数据。

然后就是简单的scRNA <- RunPCA(scRNA, features = VariableFeatures(scRNA))，绘图可以用DimPlot()。

最后挑选主成分，jackStraw()评估每个基因的PCA得分的显著性（对于每个主成分，它会对数据进行PCA，然后对PCA的结果进行随机重排，然后再次进行PCA，然后计算原始PCA（实际）和重排PCA（随机下理论）的差异，这个差异越大，说明这个主成分越重要。），ScoreJackStraw用于分析jackStraw函数的结果，JackStrawPlot()绘制jackStraw图。

这个图，empirical是实际的PCA得分，theoretical是随机重排的PCA得分（这两个对应彩色点），虚线是"均匀分布的 p 值的累积分布函数"，如果主成分轴显着，则二者应显示出差异。不是很清楚，反正选不一样的就完了。
这里结果都挺好的，盗用一下某人的例子，像11、12这种就不应该选了： 《JackStraw Plot (QQ-plot) を理解する》

或者可以画ElbowPlot，也可以选择累积解释方差≥90%的主成分数。

聚类还有一堆东西

先FindNeighbors()计算最近邻，然后用FindClusters()聚类。

这里结果可能不太一样（甚至教程的复现和原文不一样），有可能是质控的参数不一样，也有可能作者给出的数据经过质控。

FindAllMarkers()找marker基因，真慢啊，有必要换RRO或者python。
完了再筛选出P<0.05的，还有在不同类型细胞中表达量差异大的基因。

QC之二

刚才说的要先聚类或标准化的筛选双细胞和周期判断。

emptyDrops()检测空液滴，而检测双细胞可以用doubletFinder()，案例里的doubletCluster()其实寻找的是表达介于两个Cluster之间的Cluster，但他解释的是“有无double droplet聚在一起的类”（而且好像已经废弃了）。
CellCycleScoring()评估细胞周期。

细胞类型注释

按下不表。

轨迹分析

选择特定基因代表细胞不同阶段的特征，然后用Monocle包，用的2，可以去看他们的文档和论文。

exprs <- as(as.matrix(GetAssayData(seu_obj_malOS[["RNA"]], layer = "counts")), 'sparseMatrix')
featureData <- data.frame(gene_short_name = row.names(seu_obj_malOS),row.names = row.names(seu_obj_malOS))
fd <- new('AnnotatedDataFrame', data = featureData)
phenoData <- seu_obj_malOS@meta.data
pd <- new('AnnotatedDataFrame', data = phenoData)
cds <- newCellDataSet(exprs,
                      phenoData = pd, 
                      featureData = fd,
                      lowerDetectionLimit = 0.5,
                      expressionFamily = negbinomial.size()) # 其实创建cds对象也可以使用importCDS，但是好像seurat V5没法兼容

cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)

cds <- detectGenes(cds, min_expr = 0.1)
expressed_genes <- row.names(subset(fData(cds),
                                    num_cells_expressed >= 10))

diff <- differentialGeneTest(cds[expressed_genes,],
                             fullModelFormulaStr="~seurat_clusters",
                             cores=10)

ordergene <- rownames(subset(diff, qval < 0.01)) 
cds <- setOrderingFilter(cds, ordergene)

cds <- reduceDimension(cds, max_components = 2,
                       method = 'DDRTree')
                       
cds <- orderCells(cds) # 可能igraph版本不兼容，要换成2.0.3

plot_cell_trajectory(cds,color_by="Pseudotime", cell_size=1,show_backbone=TRUE)

尾声

多样本整合

用来去掉样本之间的批次效应，可以使用harmony包的RunHarmony()（很强劲），也可以使用IntegrateLayers()。

infeCNV

要先准备一份文件。
inferCNVpy快得多，可以试一下。

celltype <- cbind(colnames(seu_obj_inferCNV),as.character(seu_obj_inferCNV@meta.data$celltype))
write.table(celltype,file = "celltype_seu_obj_inferCNV.txt",
            sep = '\t',row.names=F,col.names=F,quote=F)

mtx <- as.data.frame(GetAssayData(seu_obj_inferCNV, slot='counts',assay='RNA'))

geneInfor=annoGene(rownames(mtx),"SYMBOL",'human')
geneInfor_fromEN = annoGene(rownames(mtx),"ENSEMBL",'human')
geneInfor<-rbind(geneInfor,geneInfor_fromEN)

ensembl_to_symbol <- setNames(geneInfor$SYMBOL, geneInfor$ENSEMBL)
new_rownames <- rownames(mtx)
for (i in seq_along(new_rownames)) {
  if (new_rownames[i] %in% names(ensembl_to_symbol)) {
    new_name <- ensembl_to_symbol[new_rownames[i]]
    if (!(new_name %in% new_rownames)) {
      new_rownames[i] <- new_name
    }
  }
}
rownames(mtx) <- new_rownames

geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]

geneInfor$chr_numeric <- as.numeric(gsub("chr", "", geneInfor$chr))
geneInfor = geneInfor[order(geneInfor$chr_numeric, geneInfor$start), ]
geneInfor$chr_numeric <- NULL
rownames(geneInfor)<-NULL

mtx=mtx[rownames(mtx) %in% geneInfor[,1],]
mtx=mtx[match( geneInfor[,1], rownames(mtx) ),]

geneFile='geneFile.txt'
write.table(geneInfor,file = geneFile,sep = '\t',quote = F,col.names = F,row.names = F)

infercnv_obj = CreateInfercnvObject(
  raw_counts_matrix = mtx,
  annotations_file = 'celltype_seu_obj_inferCNV.txt',
  gene_order_file= geneFile,
  ref_group_names= c("T/NK cells","B cells","Myeloid cells")) # 可以使用免疫细胞作为参考

infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=0.1, 
                             # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir = 'infercnv_large/',# dir is auto-created for storing outputs
                             cluster_by_groups=F, 
                             HMM=F, plot_steps=F)

copyKAT

细胞通讯

可以去看他们的论文。

data.input <- seu_obj[["RNA"]]$data
meta <- seu_obj@meta.data[,c("sampleID","celltype")]
colnames(meta) <-  c("samples","celltypes")

meta$celltypes <- gsub(" cells", "", meta$celltypes)
table(meta$celltypes)

ordered_indices <- order(meta$celltypes)
meta <- meta[ordered_indices, ]
data.input <- data.input[, ordered_indices]
identical(rownames(meta),colnames(data.input))

cellchat <- createCellChat(object = data.input, 
                           meta = meta, 
                           group.by = "celltypes")
CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)
dplyr::glimpse(CellChatDB$interaction)
CellChatDB.use <- CellChatDB

cellchat@DB <- CellChatDB.use
cellchat <- subsetData(cellchat) 

options(future.globals.maxSize= 10000*1024^2)
future::plan("multisession", workers = 10)
cellchat <- identifyOverExpressedGenes(cellchat)

cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat, type = "triMean",raw.use = TRUE)

cellchat <- filterCommunication(cellchat, min.cells = 10)

cellchat <- computeCommunProbPathway(cellchat)

cellchat <- aggregateNet(cellchat)

groupSize <- as.numeric(table(cellchat@idents))

netVisual_circle(cellchat@net$count, vertex.weight = groupSize,
                 weight.scale = T, label.edge= F, title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, 
                 weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")

有时间再用Python实现一下。

おわり