8成生信审稿人必问！单细胞第一张可视化图，你漏了吗？
平均阅读时长为 2分钟
单细胞基础分析第二课：文章复现
大家好，我是风。欢迎来到风风的从零开始单细胞系列一。通过前面两次推文，我们已经能够从公共数据库下载单细胞数据，然后使用R语言读取构建对象，包括scater对象和Seurat对象。如果你的进度够快，这时候应该已经拟定好一个题目，并且把数据下载完成准备开始自己的课题分析了。公共数据挖掘的关键在于找到一个好的切入点，比如我们前面文献讲解的时候提到的解决免疫治疗的难题，如果是小众领域也可以直接鉴定细胞亚型，总之一个字——“快”。我个人观点是：第一篇文章不需要特别复杂，先快速发表并且熟悉整个流程，然后再通过不断学习完善知识体系发表更高分的文章。
上一期那么重要居然没啥人看，我差点被主编大大抓去烤了！！！（单细胞都有两个“对象”，你一个都没有）所以我决定今天放几张花里胡哨的图出来，我们来聊聊数据质控（QC）。一旦我们得到了表达矩阵，每一行对应一个基因（或转录本），每一列对应一个细胞，接下来我们就应该对矩阵进行质控，让数据可信度更高，以得到更可靠的结果。举个例子：未质控的数据可能包含低质量的细胞，那么可能导致下游分析的时候掩盖了一些特定的生物学结果。需要说明的是，目前并没有完全统一的单细胞数据质控标准。所以在做项目的时候，需要根据大家具体项目来进行质控，比如对于线粒体基因，人的线粒体基因集就不适用于植物的数据质控。此外，QC阈值的设置也没有特定的标准，根据自己的矩阵选择合适的阈值即可。
我们将使用2017年发表在“Scientific Reports”上的一篇文章：“Batch effects and the effective design of single-cell gene expression studies”介绍的方法对数据进行质控，作者来自芝加哥大学。
我们使用作者的数据来进行QC，数据我放在了后台，可以回复关键词下载，也可以自己阅读文章下载对应的数据，我的电脑是8+512，跑起来需要一点点时间，虽然不长，但是需要大家有点耐性。接下来我们先读取数据：
options(stringsAsFactors = FALSE)library(SingleCellExperiment)library(scater)# Loading required package: ggplot2molecules <- read.table("molecules.txt", sep = "\t") # 加载矩阵head(molecules[1:3]) # 查看矩阵数据# NA19098.r1.A01 NA19098.r1.A02 NA19098.r1.A03# ENSG00000237683 0 0 0# ENSG00000187634 0 0 0# ENSG00000188976 3 6 1# ENSG00000187961 0 0 0# ENSG00000187583 0 0 0# ENSG00000187642 0 0 0anno <- read.table("annotation.txt", sep = "\t", header = TRUE) # 加载注释文件head(anno) # 查看注释数据# individual replicate well batch sample_id# 1 NA19098 r1 A01 NA19098.r1 NA19098.r1.A01# 2 NA19098 r1 A02 NA19098.r1 NA19098.r1.A02# 3 NA19098 r1 A03 NA19098.r1 NA19098.r1.A03# 4 NA19098 r1 A04 NA19098.r1 NA19098.r1.A04# 5 NA19098 r1 A05 NA19098.r1 NA19098.r1.A05# 6 NA19098 r1 A06 NA19098.r1 NA19098.r1.A06 使用singlecelexperimental (SCE)和scater软件包对数据进行标准化，也就是我们上一讲说的构建scater对象counts <- as.matrix(molecules)umi <- SingleCellExperiment(assays = list(counts = counts), colData = anno) 去除不表达的基因（也就是总表达量为0的基因）keep_feature <- rowSums(counts(umi) > 0) > 0umi <- umi[keep_feature, ] 使用isSpike获取特征基因（互补修复基因（ERCC）和线粒体基因）isSpike(umi, "ERCC") <- grepl("^ERCC-", rownames(umi))# Warning: 'isSpike<-' is deprecated.# Use 'isSpike<-' instead.# See help("Deprecated")# Warning: 'spikeNames' is deprecated.# See help("Deprecated")# Warning: 'isSpike' is deprecated.# See help("Deprecated")isSpike(umi, "MT") <- rownames(umi) %in% c("ENSG00000198899", "ENSG00000198727", "ENSG00000198888", "ENSG00000198886", "ENSG00000212907", "ENSG00000198786", "ENSG00000198695", "ENSG00000198712", "ENSG00000198804", "ENSG00000198763", "ENSG00000228253", "ENSG00000198938", "ENSG00000198840")# Warning: 'isSpike<-' is deprecated.# Use 'isSpike<-' instead.# See help("Deprecated")# Warning: 'spikeNames' is deprecated.# See help("Deprecated")# Warning: 'isSpike' is deprecated.# See help("Deprecated")# Warning: 'isSpike' is deprecated.# See help("Deprecated") 计算QC metrics，对细胞和基因进行数据质量控制umi <- calculateQCMetrics( umi, feature_controls = list( ERCC = isSpike(umi, "ERCC"), MT = isSpike(umi, "MT")))# Warning: 'calculateQCMetrics' is deprecated.# Use 'perCellQCMetrics' or 'perFeatureQCMetrics' instead.## Warning: 'isSpike' is deprecated.# See help("Deprecated")## Warning: 'isSpike' is deprecated.# See help("Deprecated")# Warning in calculateQCMetrics(umi, feature_controls = list(ERCC = isSpike(umi, :# spike-in set 'ERCC' overwritten by feature_controls set of the same name做了QC自然也可以对结果进行可视化，结果也可以放在补充材料：
# QC结果可视化plotHighestExprs(umi, exprs_values = "counts",colour = "batch")
# plotHighestExprs函数用来可视化高表达基因（T藕片0）的表达,图中行表示每个基因，圆圈代表这个基因在所有细胞中表达量的中位数。 也可以使用plotExprsFreqVsMean函数可视化，用于观察异常值plotExprsFreqVsMean(umi)## Warning: 'plotExprsFreqVsMean' is deprecated.# Use 'plotRowData' instead.# See help("Deprecated")# `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
# plotScater函数查看Top基因对每个细胞文库的贡献值plotScater(umi,nfeatures = 300, exprs_values = "counts",colour = "batch")
可视化的函数和方法有很多，我们后面论文实操的时候再来一一展示和讲解，现在继续质控，接下来是进行细胞质控，由于分子很少的孔可能已经破裂并且无法捕获细胞，因此我们需要移除这类细胞：
# 通过total_counts对细胞进行过滤，这里选择50000total_counts <- umi$total_countsrange(total_counts)## [1] 6028 234862hist(umi$total_counts, breaks = 100)abline(v = 50000, col = "blue")
# 基于基因总数的过滤hist(umi$total_features_by_counts, breaks = 100)abline(v = 7000, col = "red")
# 另一个衡量细胞质量的指标是ERCC spike-in rna和内源性rna之间的比率，这个比值可以用来估计捕获细胞中RNA的总量。一般具有高水平spike-in RNA的细胞起始RNA量较低，这可能是由于细胞死亡或应激导致RNA降解。plotColData( umi, x = "total_features_by_counts", y = "pct_counts_MT", # 线粒体基因 colour = "batch")
plotColData( umi, x = "total_features_by_counts", y = "pct_counts_ERCC", # ERCC基因 colour = "batch")
由图片可以看到大部分来自NA19098.r2批次的细胞的ERCC/Endo比率比较高当然我们也可以使用scater对一组QC指标进行PCA分析，然后使用自动离群值进行过滤： 进行PCA分析umi <- runColDataPCA( umi, outliers = TRUE)# Warning: 'variables' must now be explicitly specifiedreducedDimNames(umi)# [1] "PCA_coldata"outlier)# # FALSE TRUE # 815 49plotReducedDim( umi, dimred = "PCA_coldata", size_by = "total_features_by_counts", shape_by = "replicate", colour_by = "outlier")
最后我们再对基因进行过滤，保留至少存在两个细胞中的基因。注意！必须在细胞过滤后对基因进行过滤，因为一些基因可能只在质量差的细胞中检测到：
keep_feature <- nexprs( umi, byrow = TRUE, detection_limit = 1) >= 2use <- keep_featuretable(keep_feature)# keep_feature# FALSE TRUE # 4512 14214 当然我们也可以通过一些个性化的条件如线粒体基因等进行基因的过滤 对基因的表达量进行均一化处理umi <- normalize(umi)# Warning: 'normalizeSCE' is deprecated.# Use 'logNormCounts' instead.# See help("Deprecated")# Warning in .local(object, ...): using library sizes as size factors# Warning: 'centreSizeFactors' is deprecated.# See help("Deprecated")# Warning in .get_all_sf_sets(object): spike-in set 'ERCC' should have its own# size factors# Warning in .get_all_sf_sets(object): spike-in set 'MT' should have its own size# factorsplotExplanatoryVariables(umi)# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'is_cell_control' with fewer than 2 unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_100_features_feature_control' with# fewer than 2 unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_200_features_feature_control' with# fewer than 2 unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_500_features_feature_control' with# fewer than 2 unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_100_features_ERCC' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_200_features_ERCC' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_500_features_ERCC' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_50_features_MT' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_100_features_MT' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_200_features_MT' with fewer than 2# unique levels# Warning in .get_variance_explained(assay(x, exprs_values), variables =# variables, : ignoring 'pct_counts_in_top_500_features_MT' with fewer than 2# unique levels
也可以通过使用细胞文库的size factors进行缩放，使得均一化后的数值范围与count一致sizeFactors(umi) <- librarySizeFactors(umi)summary(sizeFactors(umi))# Min. 1st Qu. Median Mean 3rd Qu. Max. # 0.08848 0.74696 1.00179 1.00000 1.21631 3.44738umi <- normalize(umi)# Warning: 'normalizeSCE' is deprecated.# Use 'logNormCounts' instead.# See help("Deprecated")# Warning: 'centreSizeFactors' is deprecated.# See help("Deprecated")# Warning in .get_all_sf_sets(object): spike-in set 'ERCC' should have its own# size factors# Warning in .get_all_sf_sets(object): spike-in set 'MT' should have its own size# factors好了，今天的内容就到这里。与上一讲构建scater对象和Seurat对象一致，数据质控也有多种做法，这里我们只介绍了其中一种。相比使用Seurat质控，今天介绍的这种方法虽然看起来复杂，但是有助于我们理解整个质控的流程。如果你想学习Seurat进行质控，那也不要着急，后面我们会有专门几次推文来介绍Seurat的流程化分析，理解了今天的内容，再看Seurat就会觉得很简单了。质控的过程的阈值需要大家根据实际分析的结果进行调整，有些参数也需要不断尝试，这些就需要不断试错啦！
好啦，今天的内容就到这里，今天是第二part的第三推，我相信有的人已经掉队了U•ェ•*U，能坚持的都是好样的，因为换做我我也坚持不下来/(ㄒoㄒ)/~~。
后台回复“feng42”获取数据和代码，我们下期见啦！
单细胞分析专栏传送门
碎碎念专栏传送门（完结）
风之美图系列传送门（完结）
—
END—
撰文丨风 风
排版丨四金兄
值班 | 风间琉璃
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明：以上内容为用户推荐收藏至CareerEngine平台，其内容（含文字、图片、视频、音频等）及知识版权均属用户或用户转发自的第三方网站，如涉嫌侵权，请通知[email protected]进行信息删除。如需查看信息来源，请点击“查看原文”。如需洽谈其它事宜，请联系[email protected]。
版权声明：以上内容为用户推荐收藏至CareerEngine平台，其内容（含文字、图片、视频、音频等）及知识版权均属用户或用户转发自的第三方网站，如涉嫌侵权，请通知[email protected]进行信息删除。如需查看信息来源，请点击“查看原文”。如需洽谈其它事宜，请联系[email protected]。