明码标价之单细胞转录组的质控降维聚类分群和生物学注释
单细胞转录组的质控降维聚类分群和生物学注释例子我们在《生信技能树》和《单细胞天地》都多次分享过：人人都能学会的单细胞聚类分群注释 。
一般来说，公共数据集都会给出表达量矩阵和具体不同细胞亚群特异性基因，比如 GSE122083 数据集背后的文献，就给出来了这些分群：
NK (NKG7 and GNLY), NKT (CD3D and NKG7), CD8 T cells (CD3D and CD8A), CD4 T cells (CD3D,LDHB and IL7R), B cells (MS4A1, CD79A and CD79B), monocytes (LYZ and CD14 and/or CD16), DC (LYZ and CCR7)
但是，更大的可能性是大家需要在已经发表的单细胞数据集里面去可视化自己的基因表达量情况，结合自己的生物学背景去解释这些数据。同样的，不太可能每个人都学会代码，走我们分享过：人人都能学会的单细胞聚类分群注释 这样的教程。
这次仍然是有粉丝在在我们《生信技能树》公众号后台付费求助，希望可以复现GSE122083 数据集的质控降维聚类分群和生物学注释结果，然后探索他感兴趣的8个基因的表达量情况。
就安排学徒读了一下这个文章：Predicting bacterial infection outcomes using single cell RNA-sequencing analysis of human immune cells. Nat Commun 2019 Jul 22; PMID: 31332193 ，这个《质控降维聚类分群和生物学注释》任务安排给了学徒，感谢学徒在这个春节假期还兢兢业业完成任务！
下面是学徒的探索
0、背景
（1）在Seurat等包中，在进行挑选高变基因，PCA分析后，多使用SNN(shared nearest neighbor)算法进行单细胞聚类，然后进行TSNE或者UMAP二维可视化。
（2）在一篇文献中，作者使用另一种思路：利用k-means聚类，然后进行基于KNN(k-nearest neighbor)的可视化。
下图是我根据文献流程绘制的结果，大致流程为
单细胞表达矩阵质控、过滤； 挑选Top5000高变基因； PCA主成分分析； k-means方法聚类（可进一步对cluster完成细胞类型注释）； KNN方法进行可视化；
KNN-graph是使用igraph包进行绘制，关于igraph包的相关介绍，会在笔记第二大点介绍。
1、具体绘图流程
1.0 原始数据
文献：Predicting bacterial infection outcomes using single cell RNA-sequencing analysis of human immune cells https://doi.org/10.1038/s41467-019-11257-y
测序数据：GSE122083的GSM3454528单细胞表达矩阵
1.1 表达矩阵质控（部分参考文献过滤标准）
（1）基因去重，取表达量高的
tmp1 <- read.table("GSM3454528_naive_cells.txt.gz",
header = T,#row.names = 1,
stringsAsFactors = F)
tmp1 <- tmp1[order(apply(tmp1[,-1], 1, sum),decreasing = T),]
tmp1 <- tmp1[!duplicated(tmp1$genes),]
tmp1 <- tmp1[order(tmp1$genes),]
rownames(tmp1) <- tmp1[,1]
tmp1 <- tmp1[,-1]
（2）按照文库因子标准化
lib.size=colSums(tmp1)/median(colSums(tmp1))
tmp1.new=tmp1
for(i in 1:length(lib.size)){
print(i)
tmp1.new[,i]=tmp1[,i]/lib.size[i]
}
tmp1=as.matrix(tmp1.new)
(3) log转换
tmp1=log2(tmp1+1)
dim(tmp1)
#[1] 18405 3515
#3515个细胞，18405个基因结果
1.2 挑选Top5000高变基因，进行主成分分析
构建Seurat对象，利用 FindVariableFeatures()函数处理
library("Seurat")
scRNA = CreateSeuratObject(counts=tmp1)
scRNA <- FindVariableFeatures(scRNA, selection.method = "vst", nfeatures = 5000)
hvg.gene=VariableFeatures(scRNA)
str(hvg.gene)
#chr [1:5000] "CCL5" "IGKC" "LYZ" "IGLC2" "HLA-DRA" "GNLY" "FTH1" "CD74" ...
得到新的表达矩阵
tmp1.hvg=tmp1[rownames(tmp1) %in% hvg.gene,]
dim(tmp1.hvg)
[1] 5000 3515
prcomp()主成分分析
pca <-prcomp(t(tmp1.hvg))
dim(pca$x)
#[1] 3515 3515
pca$x[1:4,1:4]
# PC1 PC2 PC3 PC4
#AAACCTGAGCTATGCT.1 3.1894482 -0.829263 3.1659335 -0.1827105
#AAACCTGCACTTCGAA.1 2.7927029 2.527055 0.8481548 -2.9785372
#AAACCTGGTTGGTGGA.1 3.1512823 1.051601 0.7916325 0.1704211
#AAACCTGTCCATGAAC.1 0.7509956 -8.282673 -4.5208889 -0.2187663
1.3 k-means cluster聚类
单细胞聚类分析时一般较少用到k-means方法。因为这种方法需要提前指定聚类数k。 如果使用这种方法，也很简单。根据主成分分析得到的结果，使用 kmeans()函数即可。
#这里使用前20个主成分，指定聚类数k=10
clust.kmeans <- kmeans(pca$x[,1:20], centers=10)
table(clust.kmeans$cluster)
# 1 2 3 4 5 6 7 8 9 10
#418 25 148 117 350 577 660 199 361 660
1.4 KNN可视化
主要根据每个细胞的主成分属性，找到与其相距“最近”的X个细胞。再利用 igraph包将这些关系可视化，并标准cluster信息其中主要涉及连个参数。一个是选用主成分数目；一个是相距最近的多少个细胞。 下面代码为采用前50个主成分，以及top20最近细胞进行操作
#得到所有细胞两两间的距离矩阵
dist<-as.matrix(dist(pca$x[,1:20]))
dist[1:3,1:3]
# AAACCTGAGCTATGCT.1 AAACCTGCACTTCGAA.1 AACCTGGTTGGTGGA.1
#AAACCTGAGCTATGCT.1 0.000000 6.929171 6.658774
#AAACCTGCACTTCGAA.1 6.929171 0.000000 5.526362
#AAACCTGGTTGGTGGA.1 6.658774 5.526362 0.000000
# 定义具有两列的空矩阵
edges <- mat.or.vec(0,2)
# for循环为每个细胞寻找最近的20个细胞
for (i in 1:nrow(dist)){
# find closes neighbours(matches即表示最近细胞的编号)
matches <- setdiff(order(dist[i,],decreasing = F)[1:21],i) #去除细胞自己与自己的距离
# add edges
edges <- rbind(edges,cbind(i,matches))
}
head(edges, 50)
# 创建igraph对象
library(igraph)
graph <- graph_from_edgelist(edges,directed=F)
graph
如下图，该igraph对象有3515个节点(细胞)，70300条边(最近距离关系)
#颜色标记细胞分类
cols<-rainbow(10)
names(cols) <- unique(clust.kmeans$cluster)
col.clust <- cols[clust.kmeans$cluster]
#由于窗口绘图，所以保存为图片再查看
png("test1.png")
set.seed(1)
plot(graph,vertex.size=1,vertex.label=NA,vertex.frame.color=NA,vertex.color=col.clust,
edge.width=0.5,main="50PCs; k=20")
legend("topright",names(cols),col=cols,
pch=16,cex=0.5,bty='n')
dev.off()
plot会自动调用plot.igraph进行绘制。值得注意的是其中的layout参数默认为layout_nicely，即自动根据节点间关系绘制最适宜的排版，但每次绘图结果会略有差异，可设置set.seed()保证结果重现性。
编写成函数function，其中参数k设定最近细胞数
make.knn.graph<-function(D,k){
# calculate euclidean distances between cells
dist<-as.matrix(dist(D))
# make a list of edges to k nearest neighbors for each cell
edges <- mat.or.vec(0,2)
for (i in 1:nrow(dist)){
# find closes neighbours
matches <- setdiff(order(dist[i,],decreasing = F)[1:(k+1)],i)
# add edges in both directions
edges <- rbind(edges,cbind(rep(i,k),matches))
}
# create a graph from the edgelist
graph <- graph_from_edgelist(edges,directed=F)
#取消节点的边框颜色
V(graph)$frame.color <- NA
# make a layout for visualizing in 2D
set.seed(1)
#指定layout_with_fr类型布局风格
g.layout<-layout_with_fr(graph)
return(list(graph=graph,layout=g.layout))
}
函数中使用了layout_with_fr排版方式，并设置了随机种子，所以绘图结果稳定。但随着其它参数的改变，出图也会发生较大的改变。如下图分别绘制了最近5、10、30、50细胞的KNN图
ks <- c(5,10,30,50)
for (k in ks){
g.pca20 <- make.knn.graph(pca$x[,1:20],k)
# plot all 4:
print(k)
png(paste0("k",k,".png"))
plot.igraph(g.pca20$graph,layout=g.pca20$layout,vertex.color=col.clust,
vertex.size=1,vertex.label=NA,main=paste0("K",k,"--","50 PCs"))
legend("topright",names(cols),col=cols,
pch=16,cex=0.5,bty='n')
dev.off()
}
以上介绍了如何利用igraph包可视化基于KNN的单细胞聚类关系。从结果来看就是从另一种方式展现单细胞的分类结果，以及分类的准确性。
2、igraph包简介
igraph是用于进行网络关系分析的开源工具，在R中有对应的R包
2.0 igraph的两要素vertices and edges与igraph对象
vertices节点，表示同一类别的具体实例，例如人名、地名、细胞、基因..... edges边，连接两端的节点，以表示这两个节点存在联系（可进一步设置权重weight、方向direct）。 igraph对象是igraph包分析的中心，其储存着vertices与edges信息，以及对应的属性。
2.1 创建igraph对象
library(igraph)
#手动创建
g1 <- graph_from_literal( Alice-Bob-Cecil-Alice, Daniel-Cecil-Eugene,
Cecil-Gordon )
如上表示A与B，B与C，C与A，D与C，C与E，C与G存在关系
但一般都不用
graph_from_literal创建，更方便的是graph_from_edgelist,graph_from_data_frame与graph_from_adjacency_matrix三种函数更符合我们分析的需求。例如上面KNN绘图中使用的是graph_from_edgelist函数。具体可参看帮助文档。
2.2 了解igraph对象
g1
如下图igraph对象简介信息分为4行3类
第1行重点关注那两个数字，前者表示节点(实例)总数，后者表示边(关系)总数。
第2行表示属性attribute，分为三大类vertices(v) or edges(e) or graph(g)。例如本例中就只有节点vertices的names属性，为character类型。还可以进一步设置节点属性(
set_vertex_attr())比如每个人的年龄，性别；边属性(set_edge_attr)如关系等级，互相打电话数等。函数vertex_attr,VandE用于查看igraph对象属性。第3、4行则主要具体展示了边的信息。
2.3 igraph对象可视化
直接使用 plot函数即可
plot(g1)
可视化过程涉及许多参数的设置，例如节点的大小/颜色/标签...，边的宽度/曲直/颜色....，还有最重要的是igraph提供多种多样的整体布局方式layout。如不特殊设置，均按照默认值。 关于igraph绘图参数，可从3种途径进行设置。下图代码展示的是在绘图时，进行设置。
#设置随机种子，保证结果稳定
set.seed(1)
plot(g, layout=layout_with_gem, vertex.size=4,
vertex.label.dist=0.5, vertex.color="red",edge.width=2)
此外
tkplot()函数可绘制交互式结果，rglplot()函数可绘制3D结果
以上简单介绍了创建igraph对象--了解igraph对象--igraph可视化三个步骤，实际可进行多种多样的复杂分析(该包的函数帮助文档有400+页)。
到时可结合需求，再了解下这个包，例如上面提到的KNN展示。
参考链接
1、Create a single cell Graph https://nbisweden.github.io/workshop-scRNAseq/oldlabs/igraph.html
2、igraph R package https://igraph.org/r/
3、igraph R manual pages https://igraph.org/r/doc/aaa-igraph-package.html
网页工具完成不了这样个性化数据分析哦
目前没有成熟的网页工具支持这样的分析。其实呢，如果你有时间请务必学习编程基础，自由自在的探索海量的公共数据，辅助你的科研，那么：
如果你没有时间从头开始学编程，也可以委托专业的团队付费拿到同样的数据分析， 普通的单细胞公共数据集的《质控降维聚类分群和生物学注释》，仅需800元，如果是文章里面的《利用igraph包可视化基于KNN的单细胞聚类关系》，价格翻倍，但都是可以拿到全部的数据和代码哦！
需要自己读文献筛选合适的数据集 提供半个小时左右的一对一讲解单细胞数据处理背景知识。
如果需要委托，直接在我们《生信技能树》公众号留言即可，我们会安排合适的生信工程师对接具体的项目。