單細胞分析實錄(7): 差別表達分析/細胞類型註釋

前面已經講解了：
單細胞分析實錄(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

咱們接下來要把這個結果整合到test.seu@meta.data中，而後畫tsne/umap展現一下

test.seu@meta.data$CB=rownames(test.seu@meta.data)
test.seu@meta.data=merge(test.seu@meta.data,result_main_hpca,by="CB")
rownames(test.seu@meta.data)=test.seu@meta.data$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))
)

test.seu@meta.data$celltype <- plyr::mapvalues(x = as.integer(as.character(test.seu@meta.data$seurat_clusters)), from = current.cluster.ids, to = new.cluster.ids)

新的test.seu@meta.data是這樣的：

統計一下各類細胞的數目

> table(test.seu@meta.data$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(test.seu@meta.data%>%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)

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