单细胞01.就这样愉快地开始吧！做个test！附代码~
关于生信学习，我们的建议是纯生信 → R 语言 → 单细胞、空间组及多组学三步走策略。费曼学习法可能是生信入门的好方法，但是真正提高还是要靠在课题进展中不断地实践！本次单细胞 21天教程以《癌生物学》《肿瘤免疫12讲2.0》为基础，并结合经典前沿综述和单细胞文章为基础，手把手展示单细胞数据分析流程。今天是环境配置。
关于单细胞教程，大家从公众号、知乎、简书、小红书或CSDN上检索，一定能搜到很多。这里，我们借助DeepSeek和接触到的单细胞教程，做个简要归纳，供果友们参考。在尝试多种提问方法对DeepSeek之后，我们发现有几个网站必定会被推荐：①10x Genomics；②Seurat；③Scanpy；④Monocle。
10x Genomics，网址: https://www.10xgenomics.com/。10x Genomics 提供单细胞测序技术和数据分析工具，该网站提供了丰富的资源，包括数据分析软件（如Cell Ranger）和教程。适合人群：希望通过真实数据练习的用户。
Seurat，网址: https://satijalab.org/seurat/。Seurat 是一个广泛使用的 R 包，可用于单细胞RNA测序数据的分析和可视化，网站有详细文档、教程和示例代码。官方教程包含从数据预处理、降维、聚类到差异表达分析的完整流程。适合人群：初学者到中级用户。
Scanpy，网址: https://scanpy.readthedocs.io/。Scanpy 是一个基于Python的单细胞数据分析工具，适用于大规模单细胞RNA测序数据的处理和分析。网站提供详细的文档和教程，有从数据加载、预处理、可视化到高级分析的完整指南。适合人群：Python 用户或偏好 Python 的分析者。
Monocle3，网址: https://cole-trapnell-lab.github.io/monocle3/。Monocle 是一个用于单细胞RNA测序数据分析的 R 包，特别适用于细胞轨迹分析和伪时间排序。网站提供详细文档和教程。适合人群：希望系统学习单细胞组学的用户。其他当然还有很多教程、书籍、资料，有资源型网站，也有数据分析型网站。但是总体来说，都不如上述流程权威或者广泛，我们就不展开了，大家如果感兴趣，可以自行检索。由于单细胞测序数据以10x Genomics为主，而我们针对的主要是初学者，因此，我们的分享是以Seurat包展开的，数据集也是以10x Genomics的为主。
国内也有很多单细胞学习教程（数据分析型教程）。生信宝典是中国中医科学院陈同教授专注于生物信息学和单细胞数据分析，定期分享教程和工具使用指南的公众号。里面内容非常丰富，涵盖 R 语言，Linux，Python和文献解读等，单细胞教程和文献也非常全面。
博客：https://blog.sciencenet.cn/blog-118204-1248788.html
CSDN：https://blog.csdn.net/qazplm12_3?type=blog
知乎：https://www.zhihu.com/people/chen586
单细胞天地专注于单细胞测序数据分析，提供从基础到进阶的教程，是生信技能树的姊妹公众号，涵盖广泛的生物信息学内容，包括单细胞数据分析的教程和实战案例。有些文章在两个公众号上都会发布，应该是国内做单细胞教程最早、最全面的团队了。国内学生信的小伙伴，估计都看过他们的推文，或者学习使用过其中的代码。
博客：http://www.bio-info-trainee.com/
简书： https://www.jianshu.com/u/06ae70ef31bc （周运来老师）
知乎：https://www.zhihu.com/people/biotrainee
生信百晓生：https://www.xiaolvji.com/u/taoxiangjiang
随着单细胞测序的推广，2020年后的单细胞文献解读、教程和代码分享雨后春笋般涌现。虽然有些测序公司或者科研团队也推出了不少教程，但是个人感觉还是在读硕博研究生写的教程更好，因为是刚需嘛，而且是精力、创作能力做好的时候，实战性很强。当然也有很多生信教培课程，内容很棒，就是价格贵。感兴趣的话，大家自行了解。
因此，现在是资源丰富、甚至过剩的时代，用马斯克的话说，现在网络上有足够的资源让你学会任何你想学习的东西，只要你愿意并且付出努力。
代码
####----01.单细胞相关的R包安装（test一下）----####
install.packages("devtools", dependencies=T)
install.packages("BiocManager", dependencies=T)
install.packages("tidyverse", dependencies=T)
install.packages('Seurat', dependencies=T)
BiocManager::install(c("SingleR","monocle", "DESeq2"),ask = F,update = F)
BiocManager::install(c("clusterProfiler","DOSE","pheatmap"),ask = F,update = F)
BiocManager::install(c("org.Hs.eg.db","org.Mm.eg.db","org.Rn.eg.db"),ask = F,update = F)
BiocManager::install("MAST") #或者从Github安装
library(Seurat)
remotes::install_github("bnprks/BPCells", quiet = TRUE)
## 加载 R 包
# install.packages(c("Seurat","tidyverse","ggplot2"))
library(Seurat)
library(tidyverse)
library(ggplot2)
## 数据读取
counts = Read10X(data.dir = "F:/GSE176029/WL02/") ## 读取10x数据
seurat = CreateSeuratObject(counts, project="WL02",
min.cells = 3,
min.features = 500) ## 创建seurat数据
{
library(Matrix) ## 非10x数据
Counts = readMM("F:/GSE176029/WL02/matrix.mtx.gz")
barcodes = read.table("F:/GSE176029/WL02/barcodes.tsv.gz", stringsAsFactors = F)[,1]
features= read.csv("F:/GSE176029/WL02/features.tsv.gz", stringsAsFactors = F, sep = "\t", header = F)
rownames(Counts) = make.unique(features[,2])
colnames(Counts) = barcodes
}
## 质量控制
seurat[["percent.mt"]] = PercentageFeatureSet(seurat, pattern = "^MT-")
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
VlnPlot(seurat, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3,
pt.size = 0)
library(patchwork)
plot1 = FeatureScatter(seurat, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 = FeatureScatter(seurat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
seurat = subset(seurat, subset = nFeature_RNA > 500 & nFeature_RNA < 7500 & percent.mt < 20)
## 数据标准化
seurat = NormalizeData(seurat)
## 识别高变基因
seurat = FindVariableFeatures(seurat, nfeatures = 3000)
# 可视化结果
top_features = head(VariableFeatures(seurat), 20)
plot1 = VariableFeaturePlot(seurat)
plot2 = LabelPoints(plot = plot1, points = top_features, repel = TRUE)
plot1 + plot2
## 数据缩放
seurat = ScaleData(seurat)
seurat = ScaleData(seurat, vars.to.regress = c("nFeature_RNA", "percent.mt"))
##对数归一化法可能会引入零膨胀伪影，使用SCTransform
seurat = SCTransform(seurat, variable.features.n = 3000)
seurat = SCTransform(seurat,
vars.to.regress = c("nFeature_RNA", "percent.mt"),
variable.features.n = 3000)
## 线性降维
seurat = RunPCA(seurat, npcs = 50)
ElbowPlot(seurat, ndims = ncol(Embeddings(seurat, "pca")))
PCHeatmap(seurat, dims = 1:20, cells = 500, balanced = TRUE, ncol = 4)
## 非线性降维
seurat = RunTSNE(seurat, dims = 1:20)
seurat = RunUMAP(seurat, dims = 1:20)
plot1 = TSNEPlot(seurat)
plot2 = UMAPPlot(seurat)
plot1 + plot2
plot1 = FeaturePlot(seurat, c("MKI67","NES","TFAP2A"),
ncol=3, reduction = "tsne")
plot2 = FeaturePlot(seurat, c("MKI67","NES","TFAP2A"),
ncol=3, reduction = "umap")
plot1 /plot2
## 细胞聚类
seurat = FindNeighbors(seurat, dims = 1:20)
seurat = FindClusters(seurat, resolution = 1)
# 聚类数据的可视化
plot1 = DimPlot(seurat, reduction = "tsne", label = TRUE)
plot2 = DimPlot(seurat, reduction = "umap", label = TRUE)
plot1 + plot2
## 细胞类群注释
ct_markers <- c("PTPRC","CD3D",
"CD4","CD8A",
"CD38","MS4A1","CD79A",
"FCGR3A","NKG7",
"CD14","CD68",
"IL2","CD69",
"MKI67")
DoHeatmap(seurat, features = ct_markers)
## 鉴定各个亚群的Marker基因
cl_markers = FindAllMarkers(seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(1.2))
library(dplyr)
cl_markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
library(devtools)
install_github('immunogenomics/presto')
library(presto)
cl_markers_presto = wilcoxauc(seurat)
cl_markers_presto %>%
filter(logFC > log(1.2) & pct_in > 20 & padj < 0.05) %>%
group_by(group) %>%
arrange(desc(logFC), .by_group=T) %>%
top_n(n = 2, wt = logFC) %>%
print(n = 40, width = Inf)
# 热图展示
top10_cl_markers = cl_markers %>%
group_by(cluster) %>%
top_n(n = 10, wt = avg_log2FC)
DoHeatmap(seurat, features = top10_cl_markers$gene)
plot1 = FeaturePlot(seurat, c("PTPRC","CD8A"), ncol = 1)
plot2 = VlnPlot(seurat, features = c("PTPRC","CD8A"), pt.size = 0) + NoLegend()
plot1 + plot2 + plot_layout(widths = c(1, 2))
最后编辑于 03-30 · 浏览 707
