【原】肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)(三)：基因可視化

前言

Hello小伙伴們大家好，我是生信技能樹的小學(xué)徒”我才不吃蛋黃“。今天是肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)系列第三期。第二期我們對細(xì)胞亞群進(jìn)行了注釋（肺腺癌單細(xì)胞數(shù)據(jù)集GSE189357復(fù)現(xiàn)(二)：細(xì)胞注釋）。

本期，我們將使用多種可視化方法，繪制FeaturePlot，ggplot，DoHeatmap圖。

1.背景介紹

在胃癌系列推文中，我提到繪好看圖需要審美能力強(qiáng)。實(shí)際上，繪圖能力也能很好的反映我們的科研素養(yǎng)。繪圖能力與科研能力之間的關(guān)系非常緊密，因?yàn)橛行У目梢暬诳茖W(xué)研究中扮演著至關(guān)重要的角色。繪圖不僅是展示研究成果的手段，更是理解復(fù)雜數(shù)據(jù)、發(fā)現(xiàn)新知識以及向同行清晰傳達(dá)信息的核心工具。在單細(xì)胞數(shù)據(jù)可視化中，我們大部分情況下都是充當(dāng)“調(diào)包俠”，使用各種R包繪制各種圖形。

R語言是一種廣泛使用的統(tǒng)計(jì)分析和圖形表示語言，它在數(shù)據(jù)可視化方面有著強(qiáng)大的功能，但是也有其缺點(diǎn)。

R語言可視化的優(yōu)點(diǎn)有很多：

豐富的圖形庫：R擁有大量的圖形庫，如ggplot2、lattice、base plots等，這些庫提供了豐富的圖形類型和定制選項(xiàng)。

高度可定制：R的圖形可以高度定制，從顏色、形狀到布局，用戶可以精細(xì)控制圖形的每一個(gè)方面。

數(shù)據(jù)操作能力強(qiáng)：R是一種強(qiáng)大的數(shù)據(jù)分析工具，它可以在繪圖之前對數(shù)據(jù)進(jìn)行復(fù)雜的處理和分析。

集成統(tǒng)計(jì)分析：R語言在統(tǒng)計(jì)分析方面非常強(qiáng)大，可以直接在繪圖中集成統(tǒng)計(jì)測試和模型結(jié)果。

動態(tài)報(bào)告：R的可視化可以輕松地集成到動態(tài)報(bào)告中，如R Markdown和Shiny，這使得結(jié)果可以交互式地呈現(xiàn)。

社區(qū)支持：R有一個(gè)活躍的社區(qū)，用戶可以從中獲得大量的教程、論壇討論和包更新。

開源免費(fèi)：R是開源軟件，可以免費(fèi)使用和修改，這使得它在學(xué)術(shù)界和數(shù)據(jù)科學(xué)領(lǐng)域非常受歡迎。

跨平臺：R可以在多種操作系統(tǒng)上運(yùn)行，包括Windows、macOS和Linux。

R語言缺點(diǎn)：

學(xué)習(xí)曲線：對于初學(xué)者來說，R的學(xué)習(xí)曲線可能比較陡峭，尤其是對于那些沒有編程背景的用戶。

性能問題：雖然R在數(shù)據(jù)處理和可視化方面表現(xiàn)出色，但它在處理非常大的數(shù)據(jù)集時(shí)可能會變慢。

圖形更新：在R中更新圖形可能需要重新運(yùn)行整個(gè)腳本，這在迭代過程中可能會顯得繁瑣。

3D可視化有限：雖然R可以創(chuàng)建3D圖形，但與一些專門的軟件相比，它的3D可視化能力相對有限。

用戶界面：R的默認(rèn)用戶界面（RStudio除外）可能不如一些集成開發(fā)環(huán)境（IDE）友好。

依賴管理：R的包依賴管理有時(shí)可能會導(dǎo)致版本沖突，需要用戶手動解決。

專業(yè)圖形制作：雖然R可以制作高質(zhì)量的圖形，但與專業(yè)的圖形設(shè)計(jì)軟件相比，它在制作出版級別的圖形方面可能不夠直觀。

交互性：雖然R可以通過Shiny等工具創(chuàng)建交互式圖形，但這些交互性可能不如專門的交互式可視化工具（如Tableau）直觀和強(qiáng)大。

R繪圖最重要的一點(diǎn)是無法像PS一樣隨意改變圖形，只能在參數(shù)限定范圍內(nèi)修改調(diào)整。因此，要想把圖做的好看，就要選擇合適的繪圖函數(shù)，并不停的調(diào)參。所以說，理解函數(shù)是R繪圖的基礎(chǔ)。

2.可視化

首先加載R包，創(chuàng)建新的文件夾，讀取細(xì)胞注釋后的Seurat數(shù)據(jù)：

rm(list=ls())
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(SingleR)
library(celldex)
#BiocManager::install("celldex")
# library(devtools)
# install_github("arc85/singleseqgset")
library(singleseqgset)
library(devtools)
library(grid)
library(gridExtra)
getwd()
dir.create("4-plot")
setwd('4-plot/')
sce.all=readRDS( "../3-Celltype/sce_celltype.rds")
sce.all
#Idents(sce.all)

2.1 Featureplot可視化基因

#EPCAM,NKG7,LYZ,CD79A,CLDN5,DCN
marker <- c('EPCAM','NKG7','LYZ','CD79A','CLDN5','DCN')
gene = marker
FeaturePlot(sce.all,features = marker,cols = c("lightgrey" ,"#DE1F1F"),ncol=3,raster=FALSE)
ggsave('FeaturePlot_marker.pdf',width = 12,height = 8)

我們可以看到，上圖各個(gè)基因的顏色閾值范圍不一致，我們可以調(diào)參，使其保持一致：

p1 <- FeaturePlot(sce.all, features = marker, combine = FALSE,ncol=3,raster=FALSE )
#colours = c('lightgrey', "#DE1F1F")
fix.sc <- scale_color_gradientn( colours = c('lightgrey', "#DE1F1F"),  limits = c(0, 6))
#+NoLegend()+NoAxes()
p2 <- lapply(p1, function (x) x + fix.sc)
CombinePlots(p2)

批量畫基因,注意圖例的范圍不同：

FeaturePlot(sce.all, features =marker, 
            cols = c("lightgrey", 'red'),
            ncol = 3 ) & NoLegend() & NoAxes() & theme(
              panel.border = element_rect(color = "black", size = 1)
            )

Featureplot還可以把兩個(gè)基因畫在同一個(gè)圖中,看右上角可以發(fā)現(xiàn)黃色越深的地方兩個(gè)基因疊加越多：

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

結(jié)合ggplot函數(shù)，我們還可以把三個(gè)基因畫在同一個(gè)圖中：

提取tsne坐標(biāo)，并提取基因表達(dá)數(shù)據(jù)并與tsne坐標(biāo)合并：

tsne_df <- as.data.frame(sce.all@reductions$umap@cell.embeddings)
tsne_df$cluster <- as.factor(sce.all$celltype)
head(tsne_df)
gene_df <- as.data.frame(GetAssayData(object = sce.all, slot = "data")[c('S100A9','S100A8','CXCL8'), ])

ggplot繪制圖形：

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("umap_1", "umap_2", 'S100A9','S100A8','CXCL8'), aes(x = umap_1, y = umap_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()

上面的圖花里胡哨，目的是什么？個(gè)人認(rèn)為，可視化的主要目的是清晰的呈現(xiàn)數(shù)據(jù)的特征并揭示生物學(xué)意義。使用FeaturePlot將多個(gè)基因的表達(dá)繪制在同一個(gè)圖中，我們可以比較這不同基因在細(xì)胞中的共表達(dá)模式，以及探索它們在特定細(xì)胞類型或亞群中的潛在生物學(xué)相關(guān)性。

2.2 DoHeatmap繪制熱圖

DoHeatmap 在單細(xì)胞 RNA 測序數(shù)據(jù)分析中用于可視化細(xì)胞群體或類型的基因表達(dá)模式，其目的包括展示差異表達(dá)基因、識別細(xì)胞亞群、展示標(biāo)記基因、比較基因表達(dá)模式、分析細(xì)胞類型的群體特征以及探索潛在的功能和通路。這種可視化工具是解釋單細(xì)胞數(shù)據(jù)、理解細(xì)胞異質(zhì)性以及揭示潛在生物學(xué)機(jī)制的強(qiáng)大手段。

利用DoHeatmap函數(shù)繪制熱圖，可以展示不同細(xì)胞類型的top5 maker：

Idents(sce.all)
table(sce.all$celltype)
Idents(sce.all) = sce.all$celltype
sce1 = sce.all[, sce.all$celltype %in% c( 'B', 'Endothelial','Epithelial', 'Fibro','Myeloid' ,'T&NK' )]

if (!file.exists('sce.markers.csv')) {
  sce.markers <- FindAllMarkers(object = sce1, only.pos = TRUE, 
                                min.pct = 0.25, 
                                thresh.use = 0.25)
  write.csv(sce.markers,file='sce.markers.csv')
} else {
  
  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è)：

sce.Scale <- ScaleData(subset(sce1,downsample=100),
                       features = top5$gene )
DoHeatmap(sce.Scale,
          features = top5$gene ,
          # group.by = "celltype",
          assay = 'RNA', label = T)+
  scale_fill_gradientn(colors = c("white","grey","firebrick3"))

ggsave('markers_heatmap.pdf',width = 10,height = 7)

2.3 DotPlot繪制氣泡圖

top5_dotplot <- DotPlot(sce.all, features = top5$gene)+
  theme(axis.text.x = element_text(angle = 90, vjust = 1, hjust = 1))
top5_dotplot
ggsave('markers_top5_dotplot.pdf',width = 10,height = 7)

setwd('../')

3.參數(shù)解析

本期單細(xì)胞基因可視化用到的三種函數(shù)分別為FeaturePlot()、DoHeatmap()和DotPlot()。代碼參數(shù)及解析如下：

FeaturePlot()：

FeaturePlot(
object,   #Seurat對象，包含要進(jìn)行可視化的數(shù)據(jù)集。
features,  #字符向量，指定要在圖中展示的特征（例如基因或元數(shù)據(jù)列名）。
dims = c(1, 2),  #數(shù)字向量，長度為2，指定要用于繪圖的降維空間的維度（例如c(1, 2)表示第一和第二維度）。
cells = NULL,  #可選的，指定要在圖中展示的細(xì)胞的向量。默認(rèn)是所有細(xì)胞。
cols = if (blend) { c("lightgrey", "#ff0000", "#00ff00") } else {
c("lightgrey", "blue") },
pt.size = NULL,
order = FALSE,
min.cutoff = NA,
max.cutoff = NA,
reduction = NULL,
split.by = NULL,
shape.by = NULL,
slot = "data",
blend = FALSE,
blend.threshold = 0.5,
label = FALSE,
label.size = 4,
repel = FALSE,
ncol = NULL,
coord.fixed = FALSE,
by.col = TRUE,
sort.cell = NULL,
interactive = FALSE,
combine = TRUE
)

FeaturePlot函數(shù)參數(shù)解析：

object: Seurat對象，包含要進(jìn)行可視化的數(shù)據(jù)集。

features: 字符向量，指定要在圖中展示的特征（例如基因或元數(shù)據(jù)列名）。

dims: 數(shù)字向量，長度為2，指定要用于繪圖的降維空間的維度（例如c(1, 2)表示第一和第二維度）。

cells: 可選的，指定要在圖中展示的細(xì)胞的向量。默認(rèn)是所有細(xì)胞。

cols: 字符向量或顏色向量，用于定義繪圖中使用的顏色漸變。

pt.size: 點(diǎn)的大小。

order: 邏輯值，指定是否根據(jù)特征表達(dá)量的順序繪制細(xì)胞，可以幫助突出表達(dá)某特征的細(xì)胞。

min.cutoff, max.cutoff: 用于指定每個(gè)特征的表達(dá)值截止的向量。可以使用百分位數(shù)來指定截止點(diǎn)（例如，'q1', 'q99'表示1%和99%分位數(shù)）。

reduction: 指定要使用的降維技術(shù)，例如"umap"、"tsne"或"pca"。如果未指定，F(xiàn)eaturePlot會依次查找"umap"、"tsne"、"pca"中可用的結(jié)果。

split.by: 字符串，指定元數(shù)據(jù)中的一個(gè)變量，用于按該變量的不同類別拆分并分別繪制圖形。

shape.by: 可選的，允許根據(jù)某個(gè)細(xì)胞屬性改變點(diǎn)的形狀。

slot: 指定從哪個(gè)Seurat對象的槽中提取表達(dá)數(shù)據(jù)進(jìn)行可視化。

blend: 邏輯值，指定是否混合兩個(gè)特征的表達(dá)值來同時(shí)可視化它們。

blend.threshold: 設(shè)置用于混合特征表達(dá)的閾值，范圍從0到1。

label: 是否在圖上標(biāo)記群組。

label.size: 標(biāo)簽文字的大小。

repel: 邏輯值，指定是否使用標(biāo)簽排斥機(jī)制，以避免標(biāo)簽之間的重疊。

ncol: 數(shù)字，指定在使用split.by時(shí)，合并到一個(gè)圖中的列數(shù)。

coord.fixed: 邏輯值，指定是否使用固定的縱橫比繪制坐標(biāo)系。

by.col: 邏輯值，指定在分列展示時(shí)是否按列而非按行展示特征圖。

interactive: 邏輯值，指定是否啟用交互式FeaturePlot。

combine: 邏輯值，指定是否將多個(gè)特征圖合并為單個(gè)繪圖對象。如果為FALSE，則返回一個(gè)包含多個(gè)ggplot對象的列表。

DoHeatmap()

DoHeatmap(
object,
features = NULL,
cells = NULL,
group.by = "ident",
group.bar = TRUE,
group.colors = NULL,
disp.min = -2.5,
disp.max = NULL,
slot = "scale.data",
assay = NULL,
label = TRUE,
size = 5.5,
hjust = 0,
angle = 45,
raster = TRUE,
draw.lines = TRUE,
lines.width = NULL,
group.bar.height = 0.02,
combine = TRUE
)

DoHeatmap()函數(shù)參數(shù)解析：

object :  一個(gè)Seurat對象，包含要進(jìn)行可視化的數(shù)據(jù)集。

features : 要打印的特征向量，默認(rèn)為variableffeatures（object=object）

cells : 要繪制的細(xì)胞向量

group.by : 一種變量向量，用于按單元格分組；將“ident”傳遞給按單元格分組的標(biāo)識類

group.bar : 添加顯示單元格組狀態(tài)的顏色欄

group.colors : 要用于顏色欄的顏色

disp.min : 最小顯示值（以下所有值均被剪裁）

disp.max : 最大顯示值（以上所有值均被剪裁）；如果插槽為'比例數(shù)據(jù)，否則為6

slot : 要使用的數(shù)據(jù)槽，從'原始數(shù)據(jù)'、'數(shù)據(jù)'或'比例數(shù)據(jù)'

assay : 從中提取

label : 在顏色欄上方標(biāo)記單元格標(biāo)識

size : 顏色欄上方文本的大小

hjust : 顏色欄上方文本的水平對正

angle : 顏色欄上方文本的角度

raster : 如果為真，則使用//oraprdnt繪制，否則使用available。//oraprdnt在某些查看應(yīng)用程序（如預(yù)覽）上，由于光柵是如何插值的，因此可能會顯得模糊。如果遇到該問題，請將其設(shè)置為false（請注意，打印可能需要更長的時(shí)間來生成/渲染）。

draw.lines : 包括分隔組的白線

lines.width : 整數(shù)，用于調(diào)整分隔白色的寬度行。對應(yīng)每個(gè)組之間的“細(xì)胞”數(shù)量。

group.bar.height : 縮放顏色欄的高度

combine : 將繪圖合并到單個(gè)PatchworkedgPlot對象中。如果為FALSE，則返回ggplot對象的列表

DotPlot()

DotPlot(
object,
assay = NULL,
features,
cols = c("lightgrey", "blue"),
col.min = -2.5,
col.max = 2.5,
dot.min = 0,
dot.scale = 6,
idents = NULL,
group.by = NULL,
split.by = NULL,
cluster.idents = FALSE,
scale = TRUE,
scale.by = "radius",
scale.min = NA,
scale.max = NA
)

DotPlot()函數(shù)參數(shù)解析：

object : 一個(gè)Seurat對象，包含要進(jìn)行可視化的數(shù)據(jù)集。

assay : 要使用的分析名稱，默認(rèn)為活動分析

features : 特征的輸入向量，或特征向量的命名列表如果需要特征分組面板（復(fù)制舊SplitDotPlotGG的功能）

cols : 要打印的顏色：來自rcolorbrewer:的調(diào)色板的名稱：布魯爾.pal.info，一對定義漸變的顏色，或3+個(gè)定義多個(gè)漸變的顏色（如果拆分方式（已設(shè)置）

col.min : 最小縮放平均表達(dá)式閾值（everythingsmaller將設(shè)置為該值）

col.max : 最大縮放平均表達(dá)式閾值（所有較大的值都將設(shè)置為該值）

dot.min : 繪制最小點(diǎn)的單元格分?jǐn)?shù)（默認(rèn)值為0）。所有表達(dá)givengene的細(xì)胞組都沒有畫點(diǎn)。

dot.scale : 縮放點(diǎn)的大小，類似于cex

idents : 要包含在繪圖中的標(biāo)識類（默認(rèn)為all）

group.by : 將細(xì)胞分組的因子

split.by : 用于拆分組的因子（復(fù)制舊的SplitDotPlotGG的功能）；有關(guān)詳細(xì)信息，請參閱FetchData

cluster.idents : 是否根據(jù)給定的特征按層次聚類對標(biāo)識進(jìn)行排序，默認(rèn)值為FALSE

scale : 確定數(shù)據(jù)是否縮放，默認(rèn)為TRUE

scale.by : 按“大小”或“半徑”縮放點(diǎn)的大小

scale.min : 設(shè)置縮放下限，默認(rèn)使用NA

scale.max : 設(shè)置縮放上限，默認(rèn)使用NA

結(jié)語

本期，我們使用多種可視化方法，繪制DimPlot，F(xiàn)eaturePlot，ggplot，DoHeatmap圖，展示了不同細(xì)胞類型基因的表達(dá)。單細(xì)胞數(shù)據(jù)分析可視化的方法遠(yuǎn)不止以上幾種，新的繪圖R包層出不窮，各種好看的圖形讓我們目不暇接。在不停的學(xué)習(xí)代碼之余，我們還要清晰的認(rèn)識到，可視化的主要目的是清晰的呈現(xiàn)數(shù)據(jù)的特征并揭示生物學(xué)意義。所以，我們還是要不停的看教材、看文獻(xiàn)、學(xué)統(tǒng)計(jì)，不停的思考。任重而道遠(yuǎn)，我們?nèi)孕枧ΑＯ乱黄冢覀儗⒃诖嘶A(chǔ)上，繪制餅圖、堆積柱狀圖、箱線圖、氣泡圖等，比較不同分組之間細(xì)胞比例差異。干貨滿滿，歡迎大家持續(xù)追更，謝謝！