在RNA-Seq流程中,第一步是载入counts数据,我们可以使用DGEList工具或者自己组建表达矩阵。接下来就应该是筛选过低表达的基因,是必须要做的一项工作。
Filtering to remove low counts
Genes that have very low counts across all the libraries should be removed prior to downstream analysis. This is justified on both biological and statistical grounds. From biological point of view, a gene must be expressed at some minimal level before it is likely to be translated into a protein or to be considered biologically important. From a statistical point of view, genes with consistently low counts are very unlikely be assessed as significantly DE because low counts do not provide enough statistical evidence for a reliable judgement to be made. Such genes can therefore be removed from the analysis without any loss of information.
无论从生物学角度还是统计学角度,都应该去除低表达基因;生物学角度看,一个基因必须有一定的表达量,才能够转录、翻译蛋白质并有生物学意义;从统计角度看,一个过低表达的基因不能够获得显著性差异,不具有统计学意义
我们假设使用一个样本数20-25M的样本库举例 —— 好吧如果较真的话是 GSE60450
。此样本分组情况是6个组,每组两个重复。
group
B.lactating B.pregnant B.virgin L.lactating L.pregnant L.virgin
2 2 2 2 2 2
使用cpm手动筛选
As a rule of thumb, we require that a gene have a count of at least 10–15 in at least some libraries before it is considered to be expressed in the study. We could explicitly select for genes that have at least a couple of counts of 10 or more, but it is slightly better to base the filtering on count-per-million (CPM) values so as to avoid favoring genes that are expressed in larger libraries over those expressed in smaller libraries. For the current analysis, we keep genes that have CPM values above 0.5 in at least two libraries:
首先要遵循的一条原则是:一个基因的count数,至少在某些样本中应达到10-15。我们可以人工规定,只计算10以上的。但是,用cpm代替rowcounts要好那么一丢丢。cpm众所周知是消除测序深度差异。
使用cpm>0.5进行筛选,得到了1w+结果
> keep <- rowSums(cpm(y) > 0.5) >= 2
> table(keep)
keep
FALSE TRUE
10704 15653
为什么选择cpm>0.5这个数值呢?因为cpm0.5 ≈ 10/L,L是 minimum library size in millions, 而这个实验的L是20-25million。那么10/20=0.5。注意这个0.5并不是很重要,我们只是简单估计。这个数值对后续分析的影响也不会太大。
>=2 , 表示对每个基因(每行),其在至少两个libraries 中cpm要>0.5。选择2因为每组都有2两个重复。换句话说就是,至少在每组的两个重复样本中cpm>0.5, 我们就保留,否则就去掉。
使用counts直接筛选
上面使用cpm的筛选可以最大限度的保留符合条件的基因。如果我们赶时间的话,也可以直接使用counts筛选:keep <-rowSums(y$counts) > 50
使用edgeR的filterByExpr方法筛选
keep <- filterByExpr(exprSet,group = group)
这个方法比较方便,而且实际筛选出来的基因数量较少,我试着我试着去看edgeR的源码,居然没有搜到这个方法,只有介绍。
应用筛选结果
y <- y[keep, , keep.lib.sizes=FALSE]
keep.lib.sizes=FALSE
参数在过滤结束后重新计算样本库size,一般推荐是这样做的。尽管做不做对后续影响都不大...
小结
Whatever the filtering rule, it should be independent of the information in the targets file. It should not make any reference to which RNA libraries belong to which group, because doing so would bias the subsequent differential expression analysis.
无论何种筛选方法,都应该独立于目标文件进行。避免指定哪个RNA属于哪个group。否则会对下游分析产生影响