scRNA表达矩阵的基因层面检查
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！
豆豆写于19.9.22
有了单细胞数据的表达矩阵，要如何对基因进行检查和过滤呢？怎么知道这个表达矩阵质量好不好？
1 检查高表达基因
目的是看看实验设计如何，是否在文库制备、比对环节出了问题
先看下高表达基因（例如top50）的ID。通常来讲，其中应该存在许多连续表达的转录本，比如核糖体或线粒体相关。如果出现了与实验预期不相符的基因，就应该关注一下。例如存在许多的ERCC转录本，就说明在文库制备中添加了太多的spike-in RNA，如果核糖体基因缺失或出现了一些假基因，那么可能在比对过程中选择了不太好的比对结果。
# 这个代码可以借鉴，调整字体大小的参数
fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotHighestExprs(sce, n=50) + fontsize
2 过滤低丰度基因
定义
低丰度基因就是表达量为0或接近0，在下游分析中基本起不到什么统计作用的基因。在假设检验中，它们不能提供足够的证据去推翻原假设，但同时它们的存在会增加多重检验校正的计算量。另外还有一个弊端，可能用英文的解释比较好理解：The discreteness of the counts may also interfere with statistical procedures, e.g., by compromising the accuracy of continuous approximations。因此下游分析之前，一般是要先去除低丰度基因。
作者认为“最佳的过滤方法”是：并非简单采用一次过滤应用于整个分析，而是每一步根据需要去过滤。他的理由是
A more aggressive filter is usually required to remove discreteness (e.g., for normalization) compared to that required for removing underpowered tests. For hypothesis testing, the filter statistic should also be independent of the test statistic under the null hypothesis (Bourgon, Gentleman, and Huber 2010)
利用平均表达量进行过滤
有许多指标可以定义低丰度，其中最常用的就是平均表达量，利用calcAverage函数，同时它还会根据各个文库大小进行一个校正（默认设置use_size_factors=TRUE）
# 这里就只是简单统计一个平均值
ave.counts <- calcAverage(sce, use_size_factors=FALSE)
# 一个很有用的操作，就是如何在坐标轴上显示下标，看xlab如何设置的就知道了
hist(log10(ave.counts), breaks=100, main="", col="grey80",
xlab=expression(Log[10]~"average count"))
图中会观察到一个“山峰”在一片“高原”上拔地而起，这个“高原”就是大量的低丰度基因，我们过滤的话一般去掉平均表达量在1以下的，也就是图中横坐标0左边的区域
demo.keep <- ave.counts >= 1
filtered.sce <- sce[demo.keep,]
summary(demo.keep)
## Mode FALSE TRUE
## logical 33490 13206
# 可以看到去掉了3w多基因，可以考虑换一种方式，想办法多留一些基因
利用每个基因在多少细胞中表达进行过滤
这个指标和上面平均表达量是紧密相关的，因为一个基因在许多细胞中都表达，也意味着它的平均表达量不会太低，就像下面这张图中表示的一样。
num.cells <- nexprs(sce, byrow=TRUE)
smoothScatter(log10(ave.counts), num.cells, ylab="Number of cells",
xlab=expression(Log[10]~"average count"))
如果发现某个基因在很少细胞中表达，一般对这个基因是不感兴趣的，因为它很有可能是由于扩增偏差产生的“假表达量”。不过更严谨考虑，这个基因也可能就只是存在于某个罕见的细胞群体，而这个群体本身数量就很少。
这里选择宽松的过滤条件：一个基因至少在一个细胞中有表达，这个基因就留下
to.keep <- num.cells > 0
sce <- sce[to.keep,]
summary(to.keep)
## Mode FALSE TRUE
## logical 22833 23863
# 这样剩下2w多个基因
点击底部的“阅读原文”，获得更好的阅读体验哦😻
初学生信，很荣幸带你迈出第一步。
我们是生信星球，一个不拽术语、通俗易懂的生信知识平台。由于是2018年新号，竟然没有留言功能。需要帮助或提出意见请后台留言、联系微信或发送邮件到jieandze1314@gmail.com，每一条都会看到的哦~