前情提要
之前在整理单细胞数据质控的时候，整理了ncount_RNA 和nFeature_RNA辅助过滤以及根据线粒体基因进行过滤
在根据线粒体基因进行过滤的推文中，简单提到了计算并查看一下单细胞数据中表达量高的TOP50基因
有朋友对于计算代码有点疑惑，那今天一起来看看如何计算获取Top基因并可视化。
根据每个基因的相对表达量获取Top基因
这里的示例数据不是pbmc的数据，而是GSE161045这个单细胞数据集，那首先我们需要读取数据并创建seurat对象
rm(list=ls())
options(stringsAsFactors = F)
source('scRNA_scripts/lib.R')
#读取数据
dir='GSE161045_RAW/outputs/'
samples=list.files( dir )
samples
sceList = lapply(samples,function(pro){
# pro=samples[1]
print(pro)
tmp = Read10X(file.path(dir,pro ))
if(length(tmp)==2){
ct = tmp[[1]]
}else{ct = tmp}
sce =CreateSeuratObject(counts = ct ,
project = pro ,
min.cells = 5,
min.features = 300 )
return(sce)
})
do.call(rbind,lapply(sceList, dim))
sce.all=merge(x=sceList[[1]],
y=sceList[ -1 ],
add.cell.ids = samples )
# Alternate accessor function with the same result
LayerData(sce.all, assay = "RNA", layer = "counts")
sce.all <- JoinLayers(sce.all)
在生信技能树整理分享的qc.R的脚本中，就有通过计算每个基因的中位数表达量并进行排序，从而识别在细胞中平均表达量最高的50个基因的代码。为了查看top基因，我们可以将计算的代码从打包好的函数中提取出来，然后直接运行即可。
#查看top50基因及可视化
C=subset(sce.all,downsample=100)@assays$RNA$counts
dim(C)
dim(sce.all)
C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100
most_expressed <- order(apply(C, 1, median), decreasing = T)[50:1]
sorted_gene <- rownames(C)[most_expressed]
先借助AI帮我们解读一下这段代码，可以看到不是完全准确，在进行计算平均表达量的时候运算的先后级稍微有点调换
依次理解这些代码：
C=subset(sce.all,downsample=100)@assays$RNA$counts
使用subset 函数用于从seurat对象中提取子集，并使用downsample=100 参数指定了每个细胞的基因计数进行随机抽样，随后取出原始的表达counts计数矩阵
运行完得到一个包含29769基因，1200细胞（12个样品，每个样品随机取100个细胞）的矩阵
C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100
这个计算过程可能理解起来比较复杂，那我们可以依次拆分开每一个计算步骤，然后看看每个步骤运行得到的结果
a = Matrix::t(C)
b = Matrix::colSums(C)
d = a/b
e = Matrix::t(d) * 100
首先 Matrix::t(C)就是对我们提取的计数矩阵C进行转置，使得基因成为行，细胞成为列。
Matrix::colSums(C)计算矩阵C的每列的和，即每个细胞的基因计数总和
Matrix::t(C)/Matrix::colSums(C)——将转置后的矩阵的每个元素除以其对应细胞的总和，得到每个基因在每个细胞中的相对表达量
Matrix::t(Matrix::t(C)/Matrix::colSums(C))再将我们的计算结果的基因放回行，细胞放回列最后，将结果乘以100，将比例转换为百分比，方便进行后续的可视化
计算之后，我们就得到了每个基因在每个细胞中的相对表达量
most_expressed <- order(apply(C, 1, median), decreasing = T)[50:1]
首先使用apply(C, 1, median)，计算每个基因在所有细胞中的中位数表达量，然后使用order函数对中位数进行排序, decreasing = T表示降序排列，然后[50:1]选择表达量中位数最高的50个基因
sorted_gene <- rownames(C)[most_expressed]——查看这些基因
使用箱线图可视化
boxplot(as.matrix(Matrix::t(C[most_expressed, ])),
cex = 0.2, las = 2,
ylim=c(0,8),
xlab = "% total count per cell",
col = (scales::hue_pal())(50)[50:1],
horizontal = TRUE)
从矩阵 C 中选取 most_expressed 索引对应的行（即Top50基因），再使用 Matrix::t 函数将结果转置，使得基因成为列，细胞成为行。使用as.matrix将结果转换为矩阵格式用于绘图
从箱线图中我们可以看到，top基因中每个基因在每个细胞中的相对表达量情况。
可以简单查看一下某个基因，比如"MALAT1"基因在每个细胞中的相对表达量，在中位值2.206539左右
其他计算获取Top基因的方式
1. 计算每个基因在所有细胞中的平均相对表达量
#使用rowSums(C)/colSums(C)得到每个基因在所有细胞中的相对表达量的平均值
C=subset(sce.all,downsample=100)@assays$RNA$counts
f = as.matrix(rowSums(C)/colSums(C)) *100
g <- order(apply(f, 1, median), decreasing = T)[50:1]
对比两种计算方式的交集基因，共计33个
list <- as.vector(rownames(C)[g])
list1 <- as.vector(rownames(C)[most_expressed])
intersection <- Reduce(intersect, list(list1, list3))
可见仍然是"MALAT1"这个参与基因表达调控的lncRNA的基因表达量较高
C=Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100计算的是每个基因在每个细胞中的相对表达量。
as.matrix(rowSums(C)/colSums(C)) *100计算的是每个基因在所有细胞中的平均相对表达量，这种方法得到的是一个基因的总体表达水平的度量，而不是在单个细胞中的相对表达量。
2. 根据基因的阳性百分比的排序
测试的话可以简单设定一个阈值，判断是否为阳性——这边假定counts是大于0就代表阳性
data_matrix <- subset(sce.all,downsample=100)@assays$RNA$counts
# 计算每个基因的阳性细胞数
positive_cells <- apply(data_matrix, 1, function(x) sum(x > 0))
# 计算阳性百分比
positive_percentage <- positive_cells / ncol(data_matrix)
# 根据阳性百分比排序
gene_order <- order(positive_percentage, decreasing = TRUE)[50:1]
# 输出排序后的基因
sorted_gene_ids <- rownames(data_matrix)[gene_order]
也可以对比一下获得的基因情况，有37个交集基因
list1 <- as.vector(rownames(C)[most_expressed])
list2 <- as.vector(sorted_gene_ids)
intersection <- Reduce(intersect, list(list1, list2))
由于我们设定了阈值，人为将基因区分为阳性和阴性两个维度，所以再使用箱线图查看每个基因的表达量情况就没有意义。这边就不可视化了。
小结
在跑通流程之后理解每一步的代码就非常有意义，需要知道通过函数及一些计算方法得到了什么样的结果数据，基于这些结果数据可视化的图可以怎样理解。
如果你需要系统性学习生信分析，可以了解一下我们的生物信息学马拉松授课，你的生物信息学入门课。单细胞转录组数据分析，最好是有自己的计算机资源哦，我们的2024的共享服务器交个朋友福利价仍然是800。
