第一次投sci拒稿率（8成生信审稿人必问）

第一次投sci拒稿率（8成生信审稿人必问）我们使用作者的数据来进行QC，数据我放在了后台，可以回复关键词下载，也可以自己阅读文章下载对应的数据，我的电脑是8 512，跑起来需要一点点时间，虽然不长，但是需要大家有点耐性。接下来我们先读取数据：读取数据上一期那么重要居然没啥人看，我差点被主编大大抓去烤了！！！（单细胞都有两个“对象”，你一个都没有）所以我决定今天放几张花里胡哨的图出来，我们来聊聊数据质控（QC）。一旦我们得到了表达矩阵，每一行对应一个基因（或转录本），每一列对应一个细胞，接下来我们就应该对矩阵进行质控，让数据可信度更高，以得到更可靠的结果。举个例子：未质控的数据可能包含低质量的细胞，那么可能导致下游分析的时候掩盖了一些特定的生物学结果。需要说明的是，目前并没有完全统一的单细胞数据质控标准。所以在做项目的时候，需要根据大家具体项目来进行质控，比如对于线粒体基因，人的线粒体基因集就不适用于植物的数据质控。此外，QC阈值的

单细胞基础分析第二课：文章复现

更多生信干货教程

私信我回复“数据库”即可免费领取哦

大家好，我是风。欢迎来到风风的从零开始单细胞系列。通过前面两次推文，我们已经能够从公共数据库下载单细胞数据，然后使用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: ggplot2 molecules <- 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 0 anno <- 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) > 0 umi<-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. ##Seehelp("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对细胞进行过滤，这里选择50000 total_counts <- umi$total_counts range(total_counts) ## [1] 6028 234862 hist(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 specified reducedDimNames(umi) ## [1] "PCA_coldata" table(umi$outlier) ## ## FALSE TRUE ## 815 49 plotReducedDim( 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 ) >= 2 rowData(umi)$use <- keep_feature table(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 ## factors plotExplanatoryVariables(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.44738 umi <- 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ㄒ)/~~。

—END—