R包封面

該R包由國(guó)外友人Enblacar完成，目前處于預(yù)印本階段阻问，旨在提供一種簡(jiǎn)化的方式梧税，為已知的單細(xì)胞可視化生成可發(fā)布的圖形沦疾。與“審美愉悅”一詞一樣主觀称近，這是一組跨不同情節(jié)類型實(shí)施的主題修改。該軟件包也可作為個(gè)人項(xiàng)目哮塞，具有未來的增長(zhǎng)前景刨秆。

可以使用以下命令安裝此軟件包：

if(!requireNamespace("devtools", quietly = T)){
  install.packages("devtools") # If not installed.
}
devtools::install_github("enblacar/SCpubr")

同時(shí)，為了使該包正常執(zhí)行忆畅，需要安裝下列包：

• colortools

• dplyr

• enrichR

• forcats

• ggbeeswarm

• ggplot2

• ggpubr

• ggrastr

• ggrepel

• Matrix

• Nebulosa

• patchwork

• pbapply

• purrr

• rlang

• scales

• Seurat

• stringr

• svglite

• tidyr

• viridis

可以使用以下命令安裝所有軟件包：

cran_packages <- c("colortools",
                   "dplyr",
                   "enrichR",
                   "forcats",
                   "ggbeeswarm",
                   "ggplot2",
                   "ggpubr",
                   "ggrepel",
                   "Matrix",
                   "patchwork",
                   "purrr",
                   "rlang",
                   "scales",
                   "Seurat",
                   "stringr",
                   "svglite",
                   "tidyr",
                   "viridis")
install.packages(cran_packages)
bioconductor_packages <- c("Nebulosa")
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(bioconductor_packages)

seura對(duì)象的準(zhǔn)備質(zhì)控與降維

#引用R包
library(dplyr)
library(Seurat)
library(patchwork)
library(ggplot2)
library(SingleR)
library(CCA)
library(clustree)
library(cowplot)
library(monocle)
library(tidyverse)
library(SCpubr)

#設(shè)置目錄并讀取data文件夾下的10X文件
data_dir <- paste0(getwd(),"/data") #路徑必須中文
samples=list.files(data_dir)
dir=file.path(data_dir,samples)
afdata <- Read10X(data.dir = dir)
# 簡(jiǎn)單創(chuàng)建一個(gè)seurat對(duì)象“sample”衡未，每個(gè)feature至少在3細(xì)胞中表達(dá)同時(shí)每個(gè)細(xì)胞中至少200個(gè)feature被檢測(cè)到
sample <- Seurat::CreateSeuratObject(counts = afdata,
                                     project = "SeuratObject",
                                     min.cells = 3,
                                     min.features = 200)
# 計(jì)算一下線粒體與核糖體基因的百分比，在sample@meta.data添加列名"percent.mt"與"percent.rb"
sample <- Seurat::PercentageFeatureSet(sample, pattern = "^MT-", col.name = "percent.mt")
sample <- Seurat::PercentageFeatureSet(sample, pattern = "^RP", col.name = "percent.rb")
# 質(zhì)控
mask1 <- sample$nCount_RNA >= 1000
mask2 <- sample$nFeature_RNA >= 500
mask3 <- sample$percent.mt <= 20
mask <- mask1 & mask2 & mask3
sample <- sample[, mask]
# 標(biāo)準(zhǔn)化
sample <- Seurat::SCTransform(sample)
# 簡(jiǎn)單降個(gè)維
sample <- Seurat::RunPCA(sample)
sample <- Seurat::RunUMAP(sample, dims = 1:30)
sample <- Seurat::RunTSNE(sample, dims = 1:30)
# 分一下cluster
sample <- Seurat::FindNeighbors(sample, dims = 1:30)
sample <- Seurat::FindClusters(sample, resolution = 1.2)
#簡(jiǎn)單注釋一下
refdata=celldex::HumanPrimaryCellAtlasData()
afdata <- GetAssayData(sample, slot="data")
cellpred <- SingleR(test = afdata,
                    ref = refdata,
                    labels = refdata$label.main,
                    method = "cluster",
                    clusters = sample@meta.data$seurat_clusters)
metadata <- cellpred@metadata
celltype = data.frame(ClusterID = rownames(cellpred),
                      celltype = cellpred$labels,
                      stringsAsFactors = F)
#########細(xì)胞注釋后結(jié)果可視化
newLabels=celltype$celltype
names(newLabels)=levels(sample)
sample=RenameIdents(sample, newLabels)

三種算法簡(jiǎn)單可視化一下：

#簡(jiǎn)單可視化一下
p1<-DimPlot(sample, reduction = "pca")
p2<-DimPlot(sample, reduction = "tsne")
p3<-DimPlot(sample, reduction = "umap")
p1+p2+p

image.png

增加一列celltype到sample@metada里：

#增加一列為celltype的對(duì)象
kk=as.data.frame(sample@active.ident)
colnames(kk)="celltype"
identical(colnames(sample),row.names(kk))
sample$celltype <-kk$celltype
identical(sample@meta.data[,"celltype"],kk[,"celltype"])
afmeta <- sample@meta.data

metada表頭

使用SCpubr的do_DimPlot函數(shù)降維：

p4<-do_DimPlot(sample = sample,
                   plot.title = "af object",
                   reduction = "pca",
                   dims = c(2, 1)) #選擇展示的主成分家凯，這邊是PC2與PC1
p5<-do_DimPlot(sample = sample,
                   plot.title = "af object",
                   reduction = "tsne",
                   dims = c(2, 1))
p6<-do_DimPlot(sample = sample,
                   plot.title = "af object",
                   reduction = "umap",
                   dims = c(2, 1))
p4+p5+p6

image.png

圖片調(diào)整：

#更改圖例位置置頂top并修改列數(shù)為2（放在左邊則為left）
p4<-do_DimPlot(sample = sample,
                       plot.title = "My object", #標(biāo)題
                       reduction = "pca",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1))
p5<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "tsne",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1))
p6<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "umap",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1))
p4+p5+p6

image.png

換個(gè)展示形式（細(xì)胞反倒沒那么好看了）：

#使用標(biāo)簽而不是圖例
p4<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "pca",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE)
p5<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "tsne",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE)
p6<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "umap",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE)
p4+p5+p6

image.png

換個(gè)顏色吧（配色鬼才）：

# colors.use換個(gè)顏色吧缓醋，紅黃藍(lán)綠橙紫（配色鬼才）
colors <- c("Chondrocytes" = "red",
            "Endothelial_cells" = "yellow",
            "Tissue_stem_cells" = "blue",
            "Monocyte" = "green",
            "T_cells" = "Orange",
            "NK_cell" = "Purple"
            )
p4<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "pca",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE,
                       colors.use = colors)
p5<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "tsne",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE,
                       colors.use = colors)
p6<-do_DimPlot(sample = sample,
                       plot.title = "My object",
                       reduction = "umap",
                       legend.position = "top",
                       legend.ncol = 2,
                       dims = c(2, 1),
                       label = TRUE,
                       legend = FALSE,
                       colors.use = colors)
p4+p5+p6

image.png

看看多個(gè)特征（如基因，線粒體含量等）的分布：

#單個(gè)特征基因的展示
do_FeaturePlot(sample,
      features = "CD14",
        plot.title = "CD14",
        ncol = 1,
        dims = c(2, 1))
#多個(gè)特征基因查詢
do_FeaturePlot(sample,
           features = c("percent.mt","percent.rb", "CD14"),
           plot.title = "A collection of features",
           ncol = 3,
           dims = c(2, 1))

單個(gè)基因

多個(gè)基因或特征

**弱化某些細(xì)胞亞群的展示绊诲，即標(biāo)成灰色送粱，只看我們關(guān)注的細(xì)胞群之間的差異

#最后一行輸入想弱化展示的細(xì)胞名
do_FeaturePlot(sample,
                       features = "CD14",
                       plot.title = "CD14",
                       ncol = 3,
                       dims = c(2, 1),
                       idents.highlight = levels(sample)[!(levels(sample) %in% c("Monocyte","Tissue_stem_cells","T_cells","NK_cell" ))])

image.png

換個(gè)配色：

do_FeaturePlot(sample,
                       features = "percent.mt",
                       plot.title = "percent.mt",
                       ncol = 1,
                       dims = c(2, 1),
                       viridis_color_map = "A")

SCpubr包 ****viridis_color_map****參數(shù)支持8種連續(xù)變量配色：

• A - 巖漿顏色圖。
• B - 地獄色圖掂之。
• C - 等離子顏色圖抗俄。
• D - viridis 顏色圖脆丁。
• E - cividis 彩色地圖。
• F——火箭彩圖动雹。
• G - mako 彩色地圖槽卫。

• H - 渦輪彩色圖。
  image.png

小提琴圖

do_VlnPlot(sample = sample,
                   features = "ACTB",
                   boxplot_width = 0.5) #箱子的寬度

image.png

換個(gè)配色（配色鬼才）

# colors.use換個(gè)顏色吧胰蝠，紅黃藍(lán)綠橙紫（配色鬼才）
colors <- c("Chondrocytes" = "red",
            "Endothelial_cells" = "yellow",
            "Tissue_stem_cells" = "blue",
            "Monocyte" = "green",
            "T_cells" = "Orange",
            "NK_cell" = "Purple"
)
do_VlnPlot(sample = sample,
                   features = "ACTB",
                   boxplot_width = 0.5,
                   colors.use = colors)

image.png

蜂群圖

#蜂群圖,看看某個(gè)基因的含量吧
do_BeeSwarmPlot(sample = sample,
               feature_to_rank = "ACTB",
               group.by = "celltype", #按什么分組比較
                continuous_feature = T)

image.png

同樣可以使用viridis_color_map參數(shù)更換魔鬼配色

#魔鬼配色
do_BeeSwarmPlot(sample = sample,
                feature_to_rank = "ACTB",
                group.by = "celltype",
                continuous_feature = T, viridis_color_map = "A")

image.png

**畫個(gè)點(diǎn)圖吧家人們

#點(diǎn)圖歼培，將需要展示的基因復(fù)制為
genesgenes <- c("IL7R", "CCR7", "CD14", "LYZ", "S100A4", "MS4A1", "CD8A", "FCGR3A", "MS4A7", "GNLY", "NKG7", "FCER1A", "CST3", "PPBP")
do_DotPlot(sample = sample,
           features = genes,
           flip = T) #翻轉(zhuǎn)坐標(biāo)軸

image.png

當(dāng)然也可以分區(qū)塊進(jìn)行細(xì)胞注釋

#當(dāng)然也可以分區(qū)塊,這個(gè)時(shí)候genes就得是列表了
genes <- list("Chondrocytes" = c("IL7R", "CCR7"),
             "Endothelial_cells" = c("CD14", "LYZ"),
             "Monocyte" = c("CD8A"),
            "T_cells" = c("FCGR3A", "MS4A7"),
             "NK_cell" = c("GNLY", "NKG7"))
             
do_DotPlot(sample = sample,features = genes)

image.png

換個(gè)配色

#換個(gè)配色
do_DotPlot(sample = sample,
                   features = genes,
                   colors.use = c("blue","red"))

image.png

柱狀圖展示一下各細(xì)胞群的cluster含量

#柱狀圖
do_BarPlot(sample = sample,
           features = "celltype",
           group.by = "seurat_clusters",
           legend = T,
           horizontal = T,
           add.summary_labels = T,
           add.subgroup_labels = T,
           repel.summary_labels = T,
           repel.subgroup_labels = T)

image.png

只有一個(gè)樣本展示不了樣本中細(xì)胞亞群的比例了，這里放一下官方文檔的示例：

image.png

單細(xì)胞的基因集富集展示：

#基因集富集
GS <- list("MAPK single" = c("MAPK1","MAPK2","MAPK4","IL1B","AKT","IL7R", "CCR7"),
              "P53 single" = c("BAX","AKR", "BCL2","ACTB","TP53"),
              "JAK-STAT single" = c("MAPK2","MAPK4","IL1B","AKT","IL7R", "CCR7"),
              "cell cycle" = c("FCGR3A", "MS4A7","GNLY", "NKG7","CD8A"),
              "Death" = c("BAX","AKR", "BCL2","MAPK1","MAPK2"))
do_EnrichmentHeatmap(sample = sample,
                     list_genes = GS,
                     group.by = "celltype",
                     transpose = T, #是否裝置
                     column_names_rot = 90,row_names_rot = 0, #行列的寬度
                     cell_size = 5        #格子大小
                     )

image.png

關(guān)注”生信堿移“公眾號(hào)后臺(tái)回復(fù)：newafdata姊氓，即可獲得示例文件與代碼
說這么多了丐怯，快去實(shí)操一下吧