ixxmu
commented
2 months ago 单细胞分析全流程：数据下载、多样本整合、QC及DecontX by 生信菜鸟团
从我们生信技能树历年的几千个马拉松授课学员里面募集了一些优秀的创作者，某种意义来说是传承了我们生信技能树的知识整理和分享的思想！
今天的是学员一点一滴整理的授课知识点笔记哦，还有互动练习题哈，欢迎大家点击文末的阅读原文去关注我们学员的公众号哦！
准备开一个系列教程，回顾一下自己学过的单细胞知识，希望对大家有所帮助~
数据下载
首先要说一下如何找到自己想要的数据集。我一般会通过文献，在pubmed上按疾病+scRNA的策略进行检索，能够得到相关文献。绝大多数的文章都会包含下图中的信息，告诉你文章用到的数据上传到了哪个数据库，这里作者是上传到了cellxgene数据库上面。如果是上传到GEO，那么往往会提供GSE号。
我们最常用的数据库还是GEO，下载数据很简单，搜索在前面得到的GSE号就可以进入到数据集下载页面。之前的推文中已经介绍过标准10x三个文件还有矩阵形式的单细胞数据的读取方式，今天我选择了一个h5格式的数据进行分析，掌握这三种格式的读取方式基本能满足90%以上数据集的需求。这个数据集包含了5个样本，点击下图中的http开始进行下载。完成后解压到文件夹中。开始分析，具体说明在代码部分进行了注释。
读取数据
#读取h5格式数据需要用到hdf5r包
options(stringsAsFactors = F)
BiocManager::install("hdf5r")
library(Seurat)
library(hdf5r)
library(tidyverse)
library(data.table)
dir = "E:/scRNA/单细胞转录组/output/"#设置数据存放路径
setwd(dir)#工作路径设置
fs <- list.files(pattern = ".h5")#读取该目录下所有h5格式文件
sceList = lapply(fs, function(x){
# x=fs[1]
print(x)
a=Read10X_h5( x )
a[1:4,1:4]
library(stringr)
(p=str_split(x,'_',simplify = T)[,2])#把文件名按_进行拆分，保留第二部分，自行按需修改
sce <- CreateSeuratObject( a ,project = p )
sce
})
setwd('../')
最终我们得到的是一个包含多个seurat对象的list
质控
for(i in 1:length(sceList)){
sc <- sceList[[i]]
sc <- sceList[[i]]
# 计算线粒体比例
sc[["mt_percent"]] <- PercentageFeatureSet(sc, pattern = "^mt-")
# 将sc赋值给scRNAlist[[i]]
sceList[[i]] <- sc
# 删除sc
rm(sc)
}
sceList <- lapply(X = sceList, FUN = function(x){
x <- subset(x,
subset = nFeature_RNA > 300 & nFeature_RNA < 7000 &
mt_percent < 10 )})
sceList <- merge(sceList [[1]],y=sceList [c(2:length(sceList))])
table(sceList$orig.ident)
saveRDS(sceList,file = "sce.RDS")
这里结合了原文的一些质控条件只根据线粒体基因比例进行了过滤，这里的10并不是固定的，还可以根据红细胞基因比例、细胞周期等进行过滤，无非是把pattern改一改，方法大同小异，在此不再过多介绍。
DecontX去除RNA污染
这一步可做可不做，正好前两天生信技能树上写到了这个包，这里就尝试一下。作用是什么呢？我们来看下面两张图，很明显一张图中亚群之间黏连细胞很多而另一张则分群十分清晰，这种黏连散在的细胞一般认为是污染严重的低质量细胞，这时候我们可以通过DecontX去除这些细胞，使分群的界限更加清晰。
rm(list = ls())
sce <- readRDS("sce.RDS")
library(devtools)
devtools::install_github("campbio/decontX")
library(decontX)
#提取counts矩阵
view(sce)
counts <- sce@assays$RNA@counts
decontX_results <- decontX(counts)
sce$contamination <- decontX_results$contamination
###根据contamination值进行过滤，一般是保留<0.2的细胞,这里过滤掉了7000多个细胞,其实有点多，可能会丢失过多的信息。
sce_filt <- sce[,sce$contamination<0.2]
saveRDS(sce,file = "sce.RDS")
saveRDS(sce_filt,file = "sce_filt.RDS")
整合数据：harmony
原文用的是CCA，但我个人更喜欢harmony，速度快，计算资源消耗小
rm(list = ls())
sce <- readRDS("../sce_filt.RDS")
dir.create("2-harmony")
setwd("./2-harmony")
sce <- NormalizeData(sce,normalization.method = "LogNormalize",
scale.factor = 1e4)
sce <- FindVariableFeatures(sce)#找高变，可以自行修改参数，写完整的话就是FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
sce <- ScaleData(sce)#这一步的数据其实仅为PCA降维服务
sce <- RunPCA(sce,features = VariableFeatures(object = sce))
ElbowPlot(scRNA1, ndims=50, reduction="pca")#查看多少PC数比较合理，一般来讲可以选择10-50之间。PC数越大保留的信息越完整，但随之而来的是技术噪音的增大。后续我们选择了15个PC进行分析
library(harmony)
obj <- RunHarmony(sce,"orig.ident")
names(obj@reductions)
obj <- RunUMAP(obj,dims = 1:15,
reduction = "harmony")
DimPlot(obj, reduction = "umap",label = T)
sce <- obj
sce <- FindNeighbors(sce,reduction = "harmony",
dims = 1:15)
#在不同分辨率下进行聚类
for (res in c(0.01, 0.05, 0.1, 0.2 ,0.3, 0.5, 0.8, 1)) {
sce = FindClusters(sce,
resolution = res, algorithm = 1)
}
rm(obj)
head(colnames(sce))
table(sce$orig.ident)
apply(sce@meta.data[,grep("RNA_snn",colnames(sce@meta.data))],2,table)#把不同resolution的结果添加到metadata中。
DimPlot(sce,reduction = "umap",label = F,group.by = "orig.ident")
比较一下DecontX前后的分群状态，这里呢效果不是非常好，但可以很明显的感受到DecontX之后的分群要更加分明。
下一期我们来讲讲细胞注释那些事。
参考文章：使用DecontX预测和去除单细胞转录组的环境游离RNA污染 · Issue #3858 · ixxmu/mp_duty (github.com)
文末友情宣传
强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：
生物信息学马拉松授课（买一得五） ，你的生物信息学入门课 时隔5年，我们的生信技能树VIP学徒继续招生啦 2024的共享服务器交个朋友福利价仍然是800 千呼万唤始出来的独享生物信息学云服务器 永久免费的生信进阶培训（线下）
https://mp.weixin.qq.com/s/UiO7AQczrcMdKENCWkLRCQ