全网最详细单细胞保姆级分析教程

各位读者,好久不见,我又归来了,之后的一段时候我将以Rstudio分析单细胞的RNA-seq流程为主,希望各位读者朋友多多支持!

1. pbmc单样本分析

1.包的加载

library(multtest)
library(dplyr)
library(Seurat)
library(patchwork)
library(R.utils)

2. 清除环境变量

rm(list = ls))

3. 下载示例数据并加压

download.file(url='https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz','~/Desktop/pbmc/mypbmc.gz')
untar(gunzip("mypbmc.gz"))

4.读取数据

pbmc.data <- Read10X(data.dir = '~/Desktop/pbmc/filtered_gene_bc_matrices/hg19')

请添加图片描述

5.创建seurat对象

pbmc <- CreateSeuratObject(counts = pbmc.data,project = 'pbmc3k',min.cells = 3,min.features = 200)

请添加图片描述

6. 获取 RNA 测序数据的计数矩阵

counts_matrix <- GetAssayData(pbmc, assay = "RNA", slot = "counts")

7. 质控及数据可视化

添加线粒体基因的百分比,用作评估细胞质量的指标之一

pbmc[['percemt.mt']] <- PercentageFeatureSet(pbmc,pattern = 'MT-')
head(pbmc@meta.data)
#                  orig.ident nCount_RNA nFeature_RNA percemt.mt
# AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
# AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
# AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
# AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
# AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898
# AAACGCACTGGTAC-1     pbmc3k       2163          781  1.6643551

8. 数据可视化

8.1 小提琴图

VlnPlot(pbmc,features = c('nCount_RNA','nFeature_RNA','percemt.mt'),pt.size = 0)

请添加图片描述

FeatureScatter(pbmc,feature1 = 'nCount_RNA',feature2 = 'nFeature_RNA')
FeatureScatter(pbmc,feature1 = 'nCount_RNA',feature2 = 'percemt.mt')

请添加图片描述

9. 数据归一化

pbmc <- NormalizeData(pbmc,normalization.method = 'LogNormalize',scale.factor = 10000)

数据归一化是数据分析中的常见步骤，旨在消除数据集中不同特征之间的量纲差异，并使它们在相同的尺度上进行比较。通过归一化，可以提高后续分析的准确性和可靠性。在这个例子中，使用 LogNormalize 方法对数据进行对数转换，并使用缩放因子 10000 进行调整。这有助于将数据转换为更适合分析的形式。

10. 寻找特征基因

pbmc <- FindVariableFeatures(pbmc,selection.method = 'vst',nfeatures = 2000)
# 指定使用的特征选择方法为 vst（方差稳定变换）。vst 是一种常用的方法，用于识别在不同样本或条件下具有较高方差的特征。选取2000个基因# VariableFeatures(pbmc) -- 可变特征列表
top10 <- head(VariableFeatures(pbmc),10)
# [1] "PPBP"   "S100A9" "IGLL5"  "LYZ"    "GNLY"   "FTL"    "PF4"    "FTH1"   "FCER1A" "GNG11" plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1,points = top10,repel = TRUE,xnudge = 0,ynudge = 0)

请添加图片描述

11. 数据缩放 – 标准化

数据缩放是数据分析中的常见步骤，旨在将数据的特征值调整到相似的范围，以避免某些特征对分析结果的过度影响。

pbmc <- ScaleData(pbmc,features = rownames(pbmc))

12. 主成分分析(PCA)

主成分分析是一种常用的数据降维和可视化技术，它可以将高维数据投影到低维空间中，同时保留数据的主要特征和变异性。

用途

数据可视化：通过绘制PCA的前几个主成分，可以观察数据的分布和聚类情况。
特征选择：根据主成分的重要性，可以选择对数据解释贡献较大的特征。
数据降维：将高维数据投影到低维空间中，减少数据的复杂性，便于后续分析。

pbmc <- RunPCA(pbmc,features = VariableFeatures(object = pbmc))
print(pbmc[['pca']],nfeatures = 5)
PC_ 1 
Positive:  MALAT1, LTB, IL32, CD2, ACAP1 
Negative:  CST3, TYROBP, LST1, AIF1, FTL 
PC_ 2 
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DRA 
Negative:  NKG7, PRF1, CST7, GZMA, GZMB 
PC_ 3 
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1 
Negative:  PPBP, PF4, SDPR, SPARC, GNG11 
PC_ 4 
Positive:  HLA-DQA1, CD79A, CD79B, HIST1H2AC, HLA-DQB1 
Negative:  VIM, S100A8, S100A6, S100A4, S100A9 
PC_ 5 
Positive:  GZMB, FGFBP2, NKG7, GNLY, PRF1 
Negative:  LTB, VIM, AQP3, PPA1, MAL

12.1 维度载荷图

维度载荷表示每个特征在主成分中的贡献程度。通过可视化维度载荷，可以了解哪些特征对主成分的形成起到了重要作用，以及不同特征之间的相关性。

VizDimLoadings(pbmc,dims = 1:3,reduction = 'pca',ncol = 3)

请添加图片描述

12.2 散点图

DimPlot(pbmc,reduction = 'pca')

请添加图片描述

12.3 热图

热图中的颜色表示细胞在每个主成分上的表达水平，颜色越深表示表达越高，颜色越浅表示表达越低。通过观察热图，可以了解不同细胞在主成分上的分布情况，以及主成分与细胞特征之间的关系。

DimHeatmap(pbmc,dims = 1:15,cells = 500,balanced = TRUE)
# 设置为 TRUE 时，热图将根据细胞在主成分上的表达值进行平衡，以突出不同细胞之间的差异。

请添加图片描述

13. 数据评估

13.1 重抽样分析

JackStraw 函数用于评估主成分分析（PCA）中特征的显著性。通过重抽样，可以生成多个随机数据集，并计算每个特征在这些随机数据集中的统计量（如 p 值）,并将结果储存在pbmc对象中。

pbmc <- JackStraw(pbmc,num.replicate = 100)
# num.replicate = 100：指定重抽样的重复次数。在这里，设置为 100 次。

重抽样分析的结果可以用于以下方面：

特征选择：根据重抽样分析得到的 p 值，可以确定哪些特征在 PCA中具有显著的贡献。
模型评估：通过比较原始数据和重抽样数据的结果，可以评估PCA模型的稳定性和可靠性。

13.2 打分

ScoreJackStraw 函数用于计算每个主成分的 JackStraw 得分。JackStraw 得分是一种用于评估主成分中特征的显著性的指标。

pbmc <- ScoreJackStraw(pbmc,dims = 1:15)

得分结果可以用于以下方面：

特征选择：根据得分的大小，可以确定哪些主成分中的特征具有较高的显著性，从而可以选择这些特征进行进一步的分析或建模。
主成分解释：得分可以帮助解释每个主成分的生物学意义，通过查看得分较高的特征，可以了解主成分所代表的生物学过程或细胞状态。

请添加图片描述

JackStrawPlot 函数会绘制每个主成分的 JackStraw 得分分布情况。通过观察这些分布，可以评估每个主成分中特征的显著性。通常，显著的特征会在 JackStraw 得分分布中表现出较高的峰值或明显的分离。这意味着这些特征在主成分中具有较强的贡献，并且可能与数据的潜在结构或生物学意义相关。

散点图

#散点图
ElbowPlot(pbmc)

请添加图片描述

14. 聚类分析

14.1 FindNeighbors

FindNeighbors 函数会根据指定的主成分维度计算数据点之间的距离，并找到每个数据点的邻居。邻居的定义可以根据距离阈值或其他相似性度量来确定。通过运行这段代码，FindNeighbors 函数将在 pbmc 对象中查找邻居，并将结果存储在对象中。

pbmc <- FindNeighbors(pbmc,dims = 1:10)

邻居信息可以用于后续的分析:

聚类分析：可以使用邻居信息来进行聚类，将相似的数据点分组在一起。
可视化：可以绘制邻居关系图，以了解数据点之间的连接和分布情况。
差异表达分析：可以比较不同邻居组之间的基因表达差异。

14.2 FindClusters – 分群

FindClusters 函数会根据数据点之间的相似性将其分组到不同的聚类中。具体的聚类算法和相似性度量方法可能因函数的实现而有所不同。

pbmc <- FindClusters(pbmc,resolution = 0.5)
# resolution指定聚类的分辨率参数。分辨率值越高，得到的聚类数量就越多；分辨率值越低，得到的聚类数量就越少。

聚类结果可以用于进一步的分析，例如：

可视化：可以使用各种可视化方法（如 t-SNE 或 UMAP）来展示聚类结果，以便更好地理解数据的结构和分组情况。
差异表达分析：可以比较不同聚类之间的基因表达差异，以发现与聚类相关的生物学特征。
功能富集分析：可以对聚类中的基因进行功能富集分析，以了解聚类所代表的生物学过程或通路。

15. 数据可视化

15.1 umap降维

pbmc <- RunUMAP(pbmc,dims = 1:10)
DimPlot(pbmc,reduction = 'umap')

请添加图片描述

15.2 tsne降维

pbmc <- RunTSNE(pbmc,dims = 1:10)
DimPlot(pbmc,reduction = 'tsne')

请添加图片描述

16. 寻找Marker并对Cluster进行命名

pbmc.markers <- FindAllMarkers(pbmc,only.pos = TRUE,min.pct = 0.25,logfc.threshold = 0.25)
# only.pos = TRUE：表示只寻找上调的标记基因（即表达增加的基因）。如果将其设置为 FALSE，则会同时寻找上调和下调的标记基因
# min.pct = 0.25：指定标记基因在至少多少百分比的细胞中表达。这里设置为 0.25，表示标记基因必须在至少 25%的细胞中表达
# logfc.threshold = 0.25：指定 log2 倍变化阈值。只有当基因的表达在不同组之间的差异达到或超过这个阈值时，才会被认为是标记基因# 这段代码使用了 dplyr 包中的函数对 pbmc.markers 数据框进行分组，并在每个组中选择前 n 个具有最大 avg_log2FC 值的行
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2,wt = avg_log2FC)
# A tibble: 18 × 7
# Groups:   cluster [9]p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene         <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>        1 1.10e- 81       2.35 0.432 0.111 1.51e- 77 0       CCR7         2 2.53e- 48       2.12 0.336 0.108 3.46e- 44 0       PRKCQ-AS1    3 5.80e-139       7.15 0.297 0.005 7.96e-135 1       FOLR3        4 1.39e-123       6.77 0.275 0.006 1.90e-119 1       S100A12      5 4.78e- 57       2.06 0.42  0.113 6.56e- 53 2       AQP3         6 4.21e- 46       2.14 0.289 0.064 5.77e- 42 2       CD40LG       7 2.40e-276       7.26 0.562 0.009 3.29e-272 3       LINC00926    8 6.96e-238       6.96 0.484 0.007 9.55e-234 3       VPREB3       9 4.29e-201       4.92 0.599 0.049 5.88e-197 4       GZMK         
10 2.45e- 81       3.61 0.381 0.059 3.35e- 77 4       GZMH         
11 2.77e-222       5.49 0.516 0.009 3.79e-218 5       CDKN1C       
12 1.77e-175       5.95 0.377 0.005 2.43e-171 5       CKB          
13 1.09e-259       5.88 0.966 0.072 1.49e-255 6       GZMB         
14 1.05e-190       6.25 0.493 0.013 1.44e-186 6       AKR1C3       
15 1.43e-267       8.04 0.833 0.009 1.96e-263 7       FCER1A       
16 3.99e-213       8.15 0.472 0.002 5.48e-209 7       SERPINF1     
17 0              14.2  0.533 0     0         8       LY6G6F       
18 6.60e-197      13.8  0.333 0     9.04e-193 8       RP11-879F14.2

FindAllMarkers 函数将返回一个包含标记基因信息的数据框。数据框的每一行代表一个标记基因，包含以下列：

gene：基因名称。
p_val：统计学显著性检验的 p 值。
avg_log2FC：平均 log2 倍变化。
pct.1 和 pct.2：标记基因在两个组中的表达百分比。
cluster：标记基因所属的聚类（如果进行了聚类分析）。

请添加图片描述

# 观测基因在不同类群中的分布情况
VlnPlot(pbmc,features = c("NKG7","PF4"),slot = "counts",log = TRUE)

请添加图片描述

# 特征基因的表达量
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ"))

请添加图片描述

top_10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10,wt = avg_log2FC)
DoHeatmap(pbmc,features = top_10$gene)

请添加图片描述