一般来说,公共数据集都会给出表达量矩阵和具体不同细胞亚群特异性基因,比如 GSE122083 数据集背后的文献,就给出来了这些分群: - 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),
但是,更大的可能性是大家需要在已经发表的单细胞数据集里面去可视化自己的基因表达量情况,结合自己的生物学背景去解释这些数据。同样的,不太可能每个人都学会代码,走我们分享过:人人都能学会的单细胞聚类分群注释 这样的教程。 这次仍然是有粉丝在在我们《生信技能树》公众号后台付费求助,希望可以复现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二维可视化。 img(2)在一篇文献中,作者使用另一种思路:利用k-means聚类,然后进行基于KNN(k-nearest neighbor)的可视化。 img下图是我根据文献流程绘制的结果,大致流程为 - k-means方法聚类(可进一步对cluster完成细胞类型注释);
imgKNN-graph是使用igraph包进行绘制,关于igraph包的相关介绍,会在笔记第二大点介绍。
1、具体绘图流程1.0 原始数据文献:Predicting bacterial infection outcomes using single cell RNA-sequencing analysis of human immune cells https:///10.1038/s41467-019-11257-y 测序数据:GSE122083的GSM3454528单细胞表达矩阵
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]
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)
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
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() }
img以上介绍了如何利用igraph包可视化基于KNN的单细胞聚类关系。从结果来看就是从另一种方式展现单细胞的分类结果,以及分类的准确性。
2、igraph包简介igraph是用于进行网络关系分析的开源工具,在R中有对应的R包 img2.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 , V and E 用于查看igraph对象属性。 第3、4行则主要具体展示了边的信息。
2.3 igraph对象可视化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./workshop-scRNAseq/oldlabs/igraph.html 2、igraph R package https:///r/ 3、igraph R manual pages https:///r/doc/aaa-igraph-package.html
|