超棒的单细胞测序学习的宝藏公众号

最新推荐文章于 2024-09-14 21:53:48 发布
Valentine_faye
最新推荐文章于 2024-09-14 21:53:48 发布
阅读量443
点赞数 7
文章标签：学习 r语言
本文链接：https://blog.csdn.net/Valentine_faye/article/details/136473688
版权
指路生信小博士
单细胞直播一理解seurat数据结构与pbmc处理流程
是从数据结构带着大家理解，视频中讲到的一些坑，我有踩到
.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2",
             "/home/data/t040413/R/yll/usr/local/lib/R/site-library",  
             "/home/data/refdir/Rlib/", "/usr/local/lib/R/library"))

sessionInfo()
getwd()

##dir.create('~/')
setwd('~/single_cell_learning')

# 数据下载
#https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
#tar -zxvf pbmc3k_filtered_gene_bc_matrices.tar.gz

library(dplyr)
library(Seurat)
library(patchwork)



1#1 导入 PBMC 数据集----
pbmc.data <- Read10X(data.dir = "~/single_cell_learning/filtered_gene_bc_matrices/hg19/")
# 初始化 Seurat 对象，使用原始数据（non-normalized data）
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc


# pbmc=readRDS("pbmc3k_final.rds")
# pbmc@assays
# pbmc@meta.data %>%head()
# pbmc@active.ident
# table(pbmc@active.ident)
# 
# # 向量 数据框 列表 增删改查
# #SummarizedExperiment
# #表达矩阵
# 
# pbmc@assays$RNA
# pbmc@assays$RNA@counts %>%head() %>%.[,1:3]
# 
# pbmc@meta.data %>%head()
# 
# pbmc@meta.data$group=ifelse(grepl(pattern = "^AAAC",x = colnames(pbmc)   ) ,'con','exp')
# 
# pbmc@meta.data %>%tail()

# 让我们检查前三十个细胞中的一些基因
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]

## 3 x 30 sparse Matrix of class "dgCMatrix"
##                                                                    
## CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
## TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
## MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

#矩阵中的.值代表 0（未检测到分子）。由于 scRNA-seq 矩阵中的大多数值都是 0，因此 Seurat 尽可能地使用稀疏矩阵表示。这会显着节省数据的内存和速度。


2#2 将 QC 指标可视化为小提琴图----
# [[ 运算符可以向对象 metadata 添加列。这是存储 QC 统计数据的好地方
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

pbmc[["ribosomal.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^RP[SL]")
head(pbmc@meta.data)



VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)


# FeatureScatter 通常用于可视化特征与特征之间的关系，但也可以用于由对象计算的任何内容，即对象 metadata 中的列、PC 分数等。

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2


# 显示前 5 个细胞的 QC 指标
head(pbmc@meta.data, 5)

#nFeature_RNA is the number of genes detected in each cell.
# nCount_RNA is the total number of molecules detected within a cell. 
2.1 #2.1 qc---------
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)


3 #3 数据归一化----
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)



4# 4 识别高变特征（特征选择）----
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 确定前 10 个变异最大的基因
top10 <- head(VariableFeatures(pbmc), 10)

# 绘制带标签和不带标签的变异特征
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2


5 # 5 标准化数据------
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

5.1
pbmc <- ScaleData(pbmc)

5.2
head(pbmc@meta.data)
pbmc <- ScaleData(pbmc, vars.to.regress = c('ribosomal.mt', 'percent.mt'))
##这里是不考虑核糖体和线粒体百分比这两个要素来进行标准化

#SCTransform()##是之后推荐的新的标准化处理方式

##归一化，标准化做完之后，就可以做PCA降维
6#6 线性降维----
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))


# 通过几种不同的方式检查和可视化 PCA 结果
#我们选择多少个PCA去做下游呢，其实可以看elbowplot的结果
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

7# 7 确定数据集的维度------

# 注意：对于大数据集，此过程可能需要很长时间，为了方便起见，请注释掉。更多的近似技术，例如在 ElbowPlot() 中实现的技术，可用于减少计算时间
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

ElbowPlot(pbmc)##看图多久变平缓，判断降维聚类成多少个维度合适，选择标准变化偏差不大的点。这里觉得>=10还不错


8 #8 细胞聚类\分群-------

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
##resolution越大，清晰度越大，亚群越多。如果聚类的个数太多，可以降低resolution

# 查看前 5 个细胞的 cluster IDs
head(Idents(pbmc), 5)

9#9 运行非线性降维 (UMAP/tSNE)----
# 如果你还没有安装 UMAP，你可以通过 reticulate::py_install(packages = 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
##这一步小博士遇到了matrix函数在多个包里冲突的error，我这次没有遇到，之后要是遇到了就有个心理准备了
pbmc <- RunTSNE(pbmc, dims = 1:10)##UMAP, tSNE是对PCA的降维
# 请注意，您可以设置 `label = TRUE` 或使用 LabelClusters 函数来帮助标记单个类群
DimPlot(pbmc, reduction = "umap")

#saveRDS(pbmc, file = "pbmc_tutorial.rds")
#pbmc=readRDS("./pbmc3k_final.rds")
#关于保存方式，如果用saveRDS保存，则用readRDS读取，如果用save保存，就用load载入


#
10# 10寻找marker基因 cluster 2 的所有 markers----
##这个时候已经命名好了，想找谁的marker基因，就把ident.1 = '1'中的'1'改成它的名字
table(Idents(pbmc))
#看一看都取了什么名字
cluster2.markers <- FindMarkers(pbmc, ident.1 = 'B', min.pct = 0.25)

head(cluster2.markers, n = 5)
##cd79a, cd79b就是B细胞的marker

allmarkers <- FindAllMarkers(pbmc, min.pct = 1,only.pos = TRUE)
#找所有细胞类型的markers,
##这里找出来的markers是有正有负的，如果只想要正的markers，加参数'only.pos = TRUE'即可,pos是positive的意思，正的嘛
#min.pct = 1这个参数意思是在多少个细胞里面有表达，min.pct = 0.25，时运行速度会很慢，改成=1之后会加快
head(allmarkers)


10.1
# 10.1 寻找不同亚群之间的差异（基因）  区分 cluster 5 与 cluster 0 和 3 的所有 markers---------
cluster5.markers <- FindMarkers(pbmc, ident.1 = 'Memory CD4 T', ident.2 = c('Naive CD4 T', 'B'), min.pct = 0.25)
##想比较哪一群细胞就把5，0，3更改成想要的就行，比如
head(cluster5.markers, n = 5)

# 与所有剩余细胞相比，找到每个 cluster 的 markers，仅报告阳性的那些
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)

10.2 # 10.2 Seurat 有几个差异表达测试，可以使用 test.use 参数设置（有关详细信息，请参阅我们的 DE vignette）。例如，ROC 测试返回任何单个标记的“分类能力”（范围从 0-random 到 1-perfect）。----------
cluster0.markers <- FindMarkers(pbmc, ident.1 = 'B', logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster0.markers)

10.3 #10.3 热图----
pbmc.markers %>%
  group_by(cluster) %>%
  top_n(n = 10, wt = avg_log2FC) -> top10
top10
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

table(Idents(pbmc))

10.4#10.4 细胞命名----------
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
###上面代码的细胞命名方式，作者不喜欢，推荐下面的

table(pbmc@active.ident)
###看看怎么命名的
##以下是作者推荐的简单的改名方式，很方便！！----
pbmc = RenameIdents(pbmc,'Naive CD4 T' = 'NAIVE T')

saveRDS(pbmc, file = "~/gzh/featureplot_dotplot_vlnplot/pbmc3k_final.rds")

getwd()

str(pbmc)

#@layer
#数据结构对于单细胞学习的重要性

pbmc=RenameIdents(pbmc,'Naive CD4 T'='NAIVE T')
###作者总结了画图的帖子，可以去他的微信公众号找一找，生信小博士，非常感谢~