有时很难想出一个明显的过滤截止点。在这种情况下，自适应阈值可以帮助我们识别在我们用于QC的任何变量中与中位数绝对偏差（MAD）相差超过3个中位数的点。
01、单细胞数据集构建和质量控制
一旦基因表达被量化，它被总结为一个表达矩阵，每一行对应一个基因(或转录本)，每一列对应一个细胞。下一步，应检查矩阵以去除质量差的细胞。如果不能在这一阶段去除低质量的细胞，可能会增加技术噪声，从而有可能模糊下游分析中感兴趣的生物信号。
由于目前没有执行scRNA-seq的标准方法，因此本文将介绍的各种质量控制QC的期望值可能因实验而异。因此，为了执行QC，我们将寻找相对于数据集其余部分的异常值，而不是与独立的质量标准进行比较。因此，在比较使用不同测序协议的数据集的质量时，应该分情况讨论。
02、Tung数据集
我们将使用scater包，以及AnnotationDbiorg.Hs.eg.db将ENSEMBLID转换为基因名称（符号）。
library(scater)library(SingleCellExperiment)library(AnnotationDbi)library(org.Hs.eg.db)library(EnsDb.Hsapiens.v86)
接下来，我们将读取矩阵和每个细胞的注释。后者转换为因子。
molecules<-read.delim("data/tung/molecules.txt",row.names=1)annotation<-read.delim("data/tung/annotation.txt",stringsAsFactors=T)
快速浏览一下数据集：
head(molecules[,1:3])##NA19098.r1.A01NA19098.r1.A02NA19098.r1.A03##ENSG00000237683000##ENSG00000187634000##ENSG00000188976361##ENSG00000187961000##ENSG00000187583000##ENSG00000187642000head(annotation)##individualreplicatewellbatchsample_id##1NA19098r1A01NA19098.r1NA19098.r1.A01##2NA19098r1A02NA19098.r1NA19098.r1.A02##3NA19098r1A03NA19098.r1NA19098.r1.A03##4NA19098r1A04NA19098.r1NA19098.r1.A04##5NA19098r1A05NA19098.r1NA19098.r1.A05##6NA19098r1A06NA19098.r1NA19098.r1.A06
在这里，我们设置altExp包含ERCC，从主对象中删除ERCC特征：
umi<-SingleCellExperiment(assays=list(counts=as.matrix(molecules)),colData=annotation)altExp(umi,"ERCC")<-umi[grep("^ERCC-",rownames(umi)),]umi<-umi[grep("^ERCC-",rownames(umi),invert=T),]
现在，让我们将ENSEMBLID映射到基因符号。从命令中table，我们可以看到大多数基因都被注释了;但是，846返回了“NA”。默认情况下，mapIds每个ID还原一个符号;可以使用参数multiVals更改此行为。
删除所有没有找到符号的基因：
umi<-umi[!is.na(rowData(umi)$SYMBOL),]
检查一下是否可以在新注释的符号中找到线粒体蛋白。
grep("^MT-",rowData(umi)$SYMBOL,value=T)##namedcharacter(0)
奇怪的是，这什么也没返回。查找核糖体蛋白（以RPL或RPS开头）的类似命令按预期工作：
grep("^RP[LS]",rowData(umi)$SYMBOL,value=T)
快速搜索线粒体蛋白ATP8（也称为MT-ATP8）显示该名称不包含“MT-”。但是，正确的特征（ENSEMBLIDENSG00000228253）存在于我们的注释中。
grep("ATP8",rowData(umi)$SYMBOL,value=T)##ENSG00000143515ENSG00000132932ENSG00000104043ENSG00000081923ENSG00000130270##"ATP8B2""ATP8A2""ATP8B4""ATP8B1""ATP8B3"##ENSG00000124406ENSG00000228253##"ATP8A1""ATP8"
大多数现代注释，例如使用的CellRanger注释，将具有以MT-开头的线粒体基因名称。出于某种原因，我们发现的那个没有。注释问题通常很常见，应始终仔细考虑。在我们的例子中，我们也找不到基因的位置，因为染色体不受支持org.Hs.eg.db——这个数据库中没有基因组位置列：
columns(org.Hs.eg.db)##[1]"ACCNUM""ALIAS""ENSEMBL""ENSEMBLPROT""ENSEMBLTRANS"##[6]"ENTREZID""ENZYME""EVIDENCE""EVIDENCEALL""GENENAME"##[11]"GENETYPE""GO""GOALL""IPI""MAP"##[16]"OMIM""ONTOLOGY""ONTOLOGYALL""PATH""PFAM"##[21]"PMID""PROSITE""REFSEQ""SYMBOL""UCSCKG"##[26]"UNIPROT"
让我们尝试一个不同的、更详细的数据库-EnsDb.Hsapiens.v86。利用这一资源，我们可以找到位于线粒体中的13个蛋白质编码基因：
ensdb_genes<-genes(EnsDb.Hsapiens.v86)MT_names<-ensdb_genes[seqnames(ensdb_genes)=="MT"]$gene_idis_mito<-rownames(umi)%in%MT_namestable(is_mito)##is_mito##FALSETRUE##1806513
03、基本的质量控制
现在，我们可以使用将上面计算的指标添加到每个细胞和每个基因元数据的函数：
umi<-addPerCellQC(umi,subsets=list(Mito=is_mito))umi<-addPerFeatureQC(umi)
手动过滤可以使用我们选择任何截止值。为了找到一个好的值，最好看一下分布：
hist(umi$total,breaks=100)abline(v=25000,col="red")
hist(umi_cell$detected,breaks=100)abline(v=7000,col="red")
有时很难想出一个明显的过滤截止点。在这种情况下，自适应阈值可以帮助我们识别在我们用于QC的任何变量中与中位数绝对偏差（MAD）相差超过3个中位数的点。注意指定偏差的正确方向：事实上，检测到的基因数量少，但MT基因百分比高，是低质量细胞的标志：
让我们添加另一个元数据列，该列将保留有关单元格是否被丢弃的信息：umi$discard<-reasons$discard
不感兴趣
看过了
取消
人点赞
人收藏
打赏
我有话说
0/500
同步到新浪微博
您的申请提交成功
您已认证成功，可享专属会员优惠，买1年送3个月！开通会员，资料、课程、直播、报告等海量内容免费看！