單細胞轉錄組的質控降維聚類分群和生物學注釋例子我們在《生信技能樹》和《單細胞天地》都多次分享過：人人都能學會的單細胞聚類分群注釋 。
一般來說，公共數據集都會給出表達量矩陣和具體不同細胞亞群特異性基因，比如 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完成細胞類型注釋）；img
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 表達矩陣質控（部分參考文獻過濾標準）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包
img
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, 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.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的單細胞聚類關係》，價格翻倍，但都是可以拿到全部的數據和代碼哦！
提供半個小時左右的一對一講解單細胞數據處理背景知識。如果需要委託，直接在我們《生信技能樹》公眾號留言即可，我們會安排合適的生信工程師對接具體的項目。