利用igraph包可视化基于KNN的单细胞聚类关系
分享是一种态度
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展示。
References
[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
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦！
长按扫码可关注