前面已經講解了：
單細胞分析實錄(1): 認識Cell Hashing
單細胞分析實錄(2): 使用Cell Ranger得到表達矩陣
單細胞分析實錄(3): Cell Hashing資料拆分
單細胞分析實錄(4): doublet檢測
單細胞分析實錄(5): Seurat標準流程
單細胞分析實錄(6): 去除批次效應/整合資料

這一節我們進入細胞型別註釋的學習，總的來說，兩條路線，手動註釋和軟體註釋。我在實際處理的過程中，主要還是手動註釋，軟體的註釋結果只作為一個參考，來輔助證實手動註釋的結果是準確無誤的。
相信看過幾篇單細胞文獻的小夥伴基本都會知道幾種常見細胞大類的marker，我們在註釋的時候用這些marker就可以基本確定細胞大類了。另外，以SingleR為例，對於細胞大類的註釋還算準確，當然也沒有到很準確的程度，我試過更換SingleR裡面不同的參考集，最後得到的大類註釋結果一致性不到80%。對於細胞小類的註釋，軟體就更加無法勝任了，我幾乎沒有看過高分的文獻會用SingleR的小類註釋結果。當然不排除隨著單細胞研究的普及和深入，以後能有更準確的軟體出現。
接下來，我以上一節的資料為例，走一遍我的分析流程。

之前的資料呈現出了16個cluster，至於resolution引數的選擇，我的原則是能在降維圖上分開的細胞亞群，能有它自己的cluster label，並且成團較好，較緊緻的一群細胞沒有必要再強行分群（比如上圖的第4群）。
這16個cluster不一定都會用到，因為doublet、細胞數太少等原因，可能還得去掉。
我推薦用常見的marker先區分一下，大概能知道cluster對應什麼型別。這裡用到的marker都是在文獻裡面經常出現的，我自己也整理了一份更全一點的細胞型別marker清單，有需要的可以在公眾號後臺說明。

celltype_marker=c(
  "EPCAM",#上皮細胞 epithelial
  "PECAM1",#內皮細胞 endothelial
  "COL3A1",#成纖維細胞 fibroblasts
  "CD163","AIF1",#髓系細胞 myeloid
  "CD79A",#B細胞
  "JCHAIN",#漿細胞 plasma cell
  "CD3D","CD8A","CD4",#T細胞
  "GNLY","NKG7",#NK細胞
  "PTPRC"#免疫細胞
)

VlnPlot(test.seu,features = celltype_marker,pt.size = 0,ncol = 2)
ggsave(filename = "marker.png",device = "png",width = 44,height = 33,units = "cm")
 

結合對應的關係，初步確認細胞型別如下，需要說明的是，NK細胞和T細胞不是很能區分開，一些文獻也是直接把這兩個當做一個大群來做的，此處我根據GNLY基因確定第5個cluster為NK細胞，是否正確我們後面再看。

#免疫細胞
0: B_cell
4: Plasma_cell
1 2: T_cell
5: NK_cell
13: Unknown

#非免疫細胞
3 6 7 8 10 11 12: Epithelial
14: Endothelial
9: Fibroblasts

#doublet
15: Doublet (Myeloid+CD4)

以上根據經典的marker初步確定了細胞大類，接著我們找差異基因，看看找出來的每一個cluster的差異基因是不是和前面鑑定的型別一致。

1. 找差異基因

markers <- FindAllMarkers(test.seu, logfc.threshold = 0.25, min.pct = 0.1, 
                          only.pos = TRUE, test.use = "wilcox")
markers_df = markers %>% group_by(cluster) %>% top_n(n = 500, wt = avg_logFC)
write.table(markers_df,file="test.seu_0.5_logfc0.25_markers500.txt",quote=F,sep="\t",row.names=F,col.names=T)

FindAllMarkers()這個函式會將cluster中的某一類和剩下的那些類比較，來找差異基因。與其對應的有個FindMarkers()函式，可以指定哪兩個cluster對比。logfc.threshold表示logfc的閾值，這裡有兩個地方需要注意：一是Seurat裡面的logfc計算公式很特別，並不是我們平常在bulk裡面那樣算均值，相除，求log，但其實也不要糾結怎麼算的，只需知道這是表示倍數的一個指標即可；二是如果想畫火山圖，這個閾值可以設為0，不然最後畫出來的火山圖中間會缺一段，我們完全可以把全部基因先拿出來，畫火山圖，再根據p value, logFC這些閾值自己過濾。

min.pct表示基因在多少細胞中表達的閾值，only.pos = TRUE表示只求高表達的基因，test.use表示檢測差異基因所用的方法。
第2行程式碼，根據avg_logFC這一列把前500個差異基因取出來，小於500個時，有多少保留多少。
事實上，不必在意我們保留多少個差異基因，實際用到的，也就前面十幾（幾十）個基因。
最終得到的差異基因list如下：

> head(markers_df)
# A tibble: 6 x 7
# Groups:   cluster [1]
  p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene    
  <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
1     0      3.58 0.999 0.114         0 0       HLA-DRA 
2     0      2.72 1     0.331         0 0       CD74    
3     0      2.59 0.983 0.169         0 0       HLA-DPA1
4     0      2.50 0.985 0.148         0 0       HLA-DPB1

然後，我們可以根據每一個cluster排在前面的基因驗證前面鑑定的結果是否正確。我把剛剛得到的的表格看一下，基本符合前面的鑑定。

2. SingleR註釋

我剛開始用SingleR是在19年的年底，現在使用方法可能有一些變化了。因為Single在註釋細胞的過程中，會用到一些參考資料集，我們可以把這些資料集儲存下來，下一次使用就可以直接載入，省去了重新下載參考集的時間。

library(celldex)
hpca.se <- HumanPrimaryCellAtlasData()
hpca.se
## class: SummarizedExperiment 
## dim: 19363 713 
## metadata(0):
## assays(1): logcounts
## rownames(19363): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
## rowData names(0):
## colnames(713): GSM112490 GSM112491 ... GSM92233 GSM92234
## colData names(3): label.main label.fine label.ont

save(hpca.se,file="/ref/singler/HPCA/hpca.se.RData")

注意參考集裡面是logcounts矩陣，後面對於單細胞資料集也要做類似的處理。
匯入UMI count矩陣

library(SingleR)
library(scater)
library(SummarizedExperiment)
test.count=as.data.frame(test.seu[["RNA"]]@counts)

匯入參考集，保證兩個資料集的基因相同，然後log處理

load(file="hpca.se.RData")
common_hpca <- intersect(rownames(test.count), rownames(hpca.se))
hpca.se <- hpca.se[common_hpca,]
test.count_forhpca <- test.count[common_hpca,]
test.count_forhpca.se <- SummarizedExperiment(assays=list(counts=test.count_forhpca))
test.count_forhpca.se <- logNormCounts(test.count_forhpca.se)

接下來是註釋步驟，在這一步裡，我只用了main大類註釋，還有一個fine小類註釋，我就不演示了，因為我覺得小類註釋不太準。

###main
pred.main.hpca <- SingleR(test = test.count_forhpca.se, ref = hpca.se, labels = hpca.se$label.main)
result_main_hpca <- as.data.frame(pred.main.hpca$labels)
result_main_hpca$CB <- rownames(pred.main.hpca)
colnames(result_main_hpca) <- c('HPCA_Main', 'CB')
write.table(result_main_hpca, file = "HPCA_Main.txt", sep = '\t', row.names = FALSE, col.names = TRUE, quote = FALSE) #儲存下來，方便以後呼叫

得到的結果是這樣的，每個CB都有一個label

> head(result_main_hpca)
         HPCA_Main                 CB
1           B_cell A_AAACCCAAGGGTCACA
2          T_cells A_AAACCCAAGTATAACG
3 Epithelial_cells A_AAACCCAGTCTCTCAC
4           B_cell A_AAACCCAGTGAGTCAG
5          T_cells A_AAACCCAGTGGCACTC
6          T_cells A_AAACGAAAGCCAGAGT

我們接下來要把這個結果整合到[email protected]中，然後畫tsne/umap展示一下

[email protected]$CB=rownames([email protected])
[email protected]=merge([email protected],result_main_hpca,by="CB")
rownames([email protected])[email protected]$CB

p5 <- DimPlot(test.seu, reduction = "tsne", group.by = "HPCA_Main", pt.size=0.5)+theme(
  axis.line = element_blank(),
  axis.ticks = element_blank(),axis.text = element_blank()
)
p6 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident",   pt.size=0.5, label = TRUE,repel = TRUE)+theme(
  axis.line = element_blank(),
  axis.ticks = element_blank(),axis.text = element_blank()
)
fig_tsne <- plot_grid(p6, p5, labels = c('ident','HPCA_Main'),rel_widths = c(2,3))
ggsave(filename = "tsne4.pdf", plot = fig_tsne, device = 'pdf', width = 36, height = 12, units = 'cm')

整體看上去也是符合最初的鑑定的

到這兒，算是把大類鑑定了，最後我們把註釋資訊新增上去

B_cell=c(0)
Plasma_cell=c(4)
T_cell=c(1,2)
NK_cell=c(5)
Unknown=c(13)
Epithelial=c(3, 6, 7, 8, 10, 11, 12)
Endothelial=c(14)
Fibroblasts=c(9)
Doublet=c(15)

current.cluster.ids <- c(B_cell,
                         Plasma_cell,
                         T_cell,
                         NK_cell,
                         Unknown,
                         Epithelial,
                         Endothelial,
                         Fibroblasts,
                         Doublet)

new.cluster.ids <- c(rep("B_cell",length(B_cell)),
                     rep("Plasma_cell",length(Plasma_cell)),
                     rep("T_cell",length(T_cell)),
                     rep("NK_cell",length(NK_cell)),
                     rep("Unknown",length(Unknown)),
                     rep("Epithelial",length(Epithelial)),
                     rep("Endothelial",length(Endothelial)),
                     rep("Fibroblasts",length(Fibroblasts)),
                     rep("Doublet",length(Doublet))
)

[email protected]$celltype <- plyr::mapvalues(x = as.integer(as.character([email protected]$seurat_clusters)), from = current.cluster.ids, to = new.cluster.ids)

新的[email protected]是這樣的：

統計一下各種細胞的數目

> table([email protected]$celltype)

     B_cell     Doublet Endothelial  Epithelial Fibroblasts     NK_cell Plasma_cell 
       1188          22          27        2023         205         441         638 
     T_cell     Unknown 
       2167          35

我們把Doublet和Unknown去掉，畫最後的tsne圖

plotCB=as.data.frame([email protected]%>%filter(seurat_clusters!="13" &
                                                seurat_clusters!="15"))[,"CB"]
DimPlot(test.seu, reduction = "tsne", group.by = "celltype", pt.size=0.5,cells = plotCB)
ggsave(filename = "tsne5.pdf", width = 15, height = 12, units = 'cm')
saveRDS(test.seu,file = "test.seu.rds") #儲存test.seu物件，下次可以直接呼叫，不用重複以上步驟

最後分享一個找差異基因的小技巧，有時候，我們希望對自己定義的類群找差異基因，可以用下面這種方法

Idents(test.seu)="celltype"
markers2 <- FindAllMarkers(test.seu, logfc.threshold = 0.25, min.pct = 0.1, 
                          only.pos = TRUE)

因水平有限，有錯誤的地方，歡迎批評指正！