手把手教你计算FRiP score的值

在Encode的ATAC文库质控标准中，认为一个高质量的文库其FRiP score值应该大于0.03，最低也要大于0.02。FRiP全称如下

Fraction of reads in peaks

表示的是位于peak区域的reads的比例，FRiP score是一个比值，其分子是位于peak区域的reads总数，分母是比对到参考基因组上的reads总数。该值最初在chip_seq文库质控中广泛使用

如上图所示，Encode汇总了不同转录因子的chip_seq数据，绘制了peak总数与FRiP score值的分布图。从图中可以看出，FRiP score值与peak 总数呈现正相关关系，而且不同转录因子对应的FRiP score值也不尽相同。

最初制定FRiP score阈值的时候，就是一个经验阈值，观察了上万例样本的FRip score值，发现绝大多数都位于0.01以上，所以采用0.01作为阈值。对于ATAC文库而言，也是同样的思路，由于ATAC的peak总数非常多，所以FRiP的阈值也比较大。需要指出的是，尽管FRiP score的阈值看上去是一个最低标准，但是由于不同组织，细胞类型的特异性，一刀切的标准是很难满足所有情况的。对于不符合FRiP score值的样本，应当结合TSS Enrichment score值等其他指标来进一步衡量其文库质量。

介绍完了FRip Score的概念和应用，我们来看下其计算过程。从概念出发，只需要peak区域内的read总数和mapping上参考基因组的reads总数即可。

在ATAC的peak caling中，使用了TagAlign这种bed文件来存储reads的比对信息，通过这个文件也可以非常快速的计算FRiP score, 步骤如下

1. 计算比对上参考基因组的reads总数

TagAlign格式中，每一行表示一个fragment的比对情况，要统计比对上的reads总数，直接统计行数即可，代码如下

wc -l sample.tn5.shift.tagalign

2. 计算peak区域的reads总数

计算peak区域的reads数目，实际上可以转换为peak区域与TagAlign这种bed文件取交集的操作，统计交集的行数即可，代码如下

bedtools intersect -a sample.tn5.shift.tagalign -b sample.narrowpeak -wa -u | wc -l

将两个数相除就得到了FRiP score, 需要注意的是，有些人会纠结fragment和read的概念，对于双端测序而言，一个fragment会产生两条reads, 上述的计算过程是针对fragment进行计数的，和reads数目的计算结果是有出入的。在我看来，FRiP score的核心思想是看peak区域的序列占所有比对上基因组序列的比例，用fragment也是对的，没必要精细到read, 而且用fragment可以同时使用单端和双端测序。

·end·