前言

Hello小伙伴們大家好，我是生信技能樹的小學(xué)徒”我才不吃蛋黃“。今天是胃癌單細(xì)胞數(shù)據(jù)集GSE163558復(fù)現(xiàn)系列第四期。第三期我們選擇0.5分辨率，對細(xì)胞進(jìn)行了分群注釋。本期，我們將在第三期基礎(chǔ)上使用多種方法可視化細(xì)胞和基因。

1.背景介紹

繪圖質(zhì)量一定程度上會(huì)直接影響文章的發(fā)表。站在審稿人的角度，好看的圖會(huì)令人賞心悅目，不好看的圖會(huì)成為巨大的減分項(xiàng)。好看的圖往往邏輯清晰、布局合理、配色協(xié)調(diào)。

邏輯清晰與否主要體現(xiàn)在機(jī)制/流程圖上。這張圖是全文的重中之重，是作者科研能力、繪圖水平和文章質(zhì)量的集中體現(xiàn)。布局是否合理大圖和小圖都會(huì)涉及。小圖需要我們?nèi)ゲ粩嗟恼{(diào)整坐標(biāo)軸、標(biāo)簽、文字的位置、大小和粗細(xì)，還有圖片的透明度、分組的順序等等等等。配色協(xié)調(diào)考驗(yàn)的是繪圖人的審美水平。筆者是直男審美，在學(xué)習(xí)R繪圖的時(shí)候，數(shù)據(jù)和代碼的問題往往能解決，但就是畫的一手丑圖。

特別是配色一言難盡，只能找愛逛街的師妹幫忙配色，或者去找高分文章的配色然后直接copy。由于是審美不在線，因此本系列所有的顏色代碼，我就老老實(shí)實(shí)的copy。

看看第三期的圖，配色是不是看著很舒服：

2.可視化

在分群注釋之后，我們可以使用DimPlot、FeaturePlot、DoHeatmap、DotPlot等多種函數(shù)對細(xì)胞或基因可視化。

在繪圖之前，我們首先創(chuàng)建新的工作目錄，并讀取第三期亞群注釋后的數(shù)據(jù)：

dir.create("4-plot")
setwd('4-plot/')
sce.all=readRDS( "../3-Celltype/sce_celltype.rds")
sce.all

使用DimPlot函數(shù)展示T細(xì)胞（"0","1"）在tSNE圖中的位置：

Idents(sce.all)
DimPlot(sce.all,label=T,
        reduction = "tsne",pt.size = 0.2,
        cells.highlight=list(
          T=WhichCells(sce.all, idents = c("0","1"))
        ),
        cols.highlight = c("red"),cols = "grey")

使用FeaturePlot函數(shù)展示T細(xì)胞maker基因('CD3D', 'CD3E', 'CD8A', 'CD4','CD2')的表達(dá)情況：

T_marker <- c('CD3D', 'CD3E', 'CD8A', 'CD4','CD2')
FeaturePlot(sce.all,features = T_marker,reduction = "tsne",cols = c("lightgrey" ,"#DE1F1F"),ncol=2,raster=FALSE)

使用DimPlot函數(shù)展示上皮細(xì)胞("5","11","12")在tSNE圖中的位置：

DimPlot(sce.all,label=T,
        reduction = "tsne",pt.size = 0.2,
        cells.highlight=list(
          Epithelial=WhichCells(sce.all, idents = c("5","11","12"))
        ),
          cols.highlight = c("red"),cols = "grey")

使用FeaturePlot函數(shù)展示上皮細(xì)胞maker基因('EPCAM','KRT19','CLDN4')的表達(dá)情況：

Epithelial_marker <- c('EPCAM','KRT19','CLDN4')
FeaturePlot(sce.all,features = Epithelial_marker,reduction = "tsne",cols = c("lightgrey" ,"#DE1F1F"),ncol=2,raster=FALSE)

FeaturePlot函數(shù)展示maker基因的表達(dá)也可以幫助我們檢查注釋準(zhǔn)確與否，因此，我們可以批量繪制多個(gè)maker基因：

genes_to_check = c('EPCAM','KRT19','CLDN4',  #上皮
                   'PECAM1' , 'CLO1A2', 'VWF',  #基質(zhì)
                   'CD3D', 'CD3E', 'CD8A', 'CD4','CD2', #T
                   'CDH5', 'PECAM1', 'VWF',  #內(nèi)皮
                   'LUM' , 'FGF7', 'MME',  #成纖維
                   'AIF1', 'C1QC','C1QB','LYZ',  #巨噬
                   'MKI67', 'STMN1', 'PCNA',  #增殖
                   'CPA3' ,'CST3', 'KIT', 'TPSAB1','TPSB2',#肥大
                   'GOS2', 'S100A9','S100A8','CXCL8', #中性粒細(xì)胞
                   'KLRD1', 'GNLY', 'KLRF1','AREG', 'XCL2','HSPA6', #NK
                   'MS4A1','CD19', 'CD79A','IGHG1','MZB1', 'SDC1',  #B
                   'CSF1R', 'CSF3R', 'CD68') #髓系

##批量畫基因
FeaturePlot(sce.all, features =genes_to_check, 
            cols = c("lightgrey", 'red'),
            reduction = "tsne",
            ncol = 8 ) & NoLegend() & NoAxes() & theme(
              panel.border = element_rect(color = "black", size = 1)
            )

FeaturePlot除了可以展示單個(gè)基因，還可以把多個(gè)基因畫在同一個(gè)圖中。

Featureplot把兩個(gè)基因畫在同一個(gè)圖中,看右上角可以發(fā)現(xiàn)黃色越深的地方兩個(gè)基因疊加越多。這樣可以看兩個(gè)基因的共表達(dá)，是不是和免疫熒光圖比較類似：

FeaturePlot(sce.all, features = c('S100A9','S100A8'),
            reduction = "tsne",
            cols = c("lightgrey", "green", "orange"),
            blend=T,blend.threshold=0)

Featureplot還可以把三個(gè)基因畫在同一個(gè)圖中：

# 提取tsne坐標(biāo)
tsne_df <- as.data.frame(sce.all@reductions$tsne@cell.embeddings)
tsne_df$cluster <- as.factor(sce.all$celltype)
head(tsne_df)
# 提取基因表達(dá)數(shù)據(jù)并與tsne坐標(biāo)合并
gene_df <- as.data.frame(GetAssayData(object = sce.all, slot = "data")[c('S100A9','S100A8','CXCL8'), ])

library(ggnewscale)
merged_df <- merge(t(gene_df), tsne_df, by = 0, all = TRUE)
head(merged_df)
colnames(merged_df)
ggplot(merged_df, vars = c("tSNE_1", "tSNE_2", 'S100A9','S100A8','CXCL8'), aes(x = tSNE_1, y = tSNE_2, colour = S100A9)) +
  geom_point(size=0.3, alpha=1) +
  scale_colour_gradientn(colours = c("lightgrey", "green"), limits = c(0, 0.3), oob = scales::squish) +
  new_scale_color() +
  geom_point(aes(colour = S100A8), size=0.3, alpha=0.7) +
  scale_colour_gradientn(colours = c("lightgrey", "blue"), limits = c(0.1, 0.2), oob = scales::squish) +
  new_scale_color() +
  geom_point(aes(colour = CXCL8), size=0.3, alpha=0.1) +
  scale_colour_gradientn(colours = c("lightgrey", "red"), limits = c(0, 0.3), oob = scales::squish)+
  theme_classic()

此外，我們還可以使用DoHeatmap函數(shù)繪制熱圖，展示每群細(xì)胞top基因：

Idents(sce.all)
Idents(sce.all) = sce.all$celltype
sce.markers <- FindAllMarkers(object = sce.all, only.pos = TRUE, 
                              min.pct = 0.25, 
                              thresh.use = 0.25)
write.csv(sce.markers,file='sce.markers.csv')
#sce.markers = read.csv('sce.markers.csv',row.names = 1)
library(dplyr) 
top5 <- sce.markers %>% group_by(cluster) %>% top_n(5, avg_log2FC)
#為了防止數(shù)據(jù)量太大不好出圖，這里在每個(gè)亞群提取出來100個(gè)
DoHeatmap(subset(sce.all,downsample=100),top5$gene,size=3)+scale_fill_gradientn(colors=c("#94C4E1","white","red"))

使用DotPlot函數(shù)繪制氣泡圖：

top5_dotplot <- DotPlot(sce.all, features = top5$gene)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1))
top5_dotplot

結(jié)語

本期，我們使用了DimPlot、FeaturePlot、DoHeatmap、DotPlot等多種函數(shù)對細(xì)胞群和基因進(jìn)行了可視化。下一期，我們將在此基礎(chǔ)上，繪制餅圖、堆積柱狀圖、箱線圖、氣泡圖等，比較不同分組之間細(xì)胞比例差異。干貨滿滿，歡迎大家持續(xù)追更，謝謝！