上期我在 cellranger定量结果详解 记录了我手动计算Sequencing Saturation 时遇到的问题
根据文章作者给的公式计算，如果grep 'xf:i:25'就拿到了唯一reads，那么sequencing saturation为50%
如果根据后面使用
samtools flagstat得到的汇总信息，duplicates得到的为非唯一reads，那么sequencing saturation为25%都不为33.7%，所以这里我还是没太弄清楚，也希望有知道的老师可以在评论区告诉我
评论区有老师给出了回答：
正好我这几天在学习lncRNA系统复习了一下转录组上游fastqc报告
里面有个Sequence Duplication Levels部分也让我一开始有点迷惑
这期推文就联系起来，系统地讲讲duplicated deduplicated 这都是啥
参考这篇文章 一文带你读懂 10X Cellranger Count 网页结果解析 我已经通过samtools view pbmc_1k_v3_possorted_genome_bam.bam | grep 'xf:i:25' | wc -l得到了n_deduped_reads也就是只被检测到一次的reads数（12435096）
通过samtools flagstat pbmc_1k_v3_possorted_genome_bam.bam得到了total reads和duplicates：
这里就需要重点提一下我当时理解错误的地方：
根据10X官网给的建议，饱和度=1-（unique_confidently_mapped_reads）/ (unique_confidently_reads+duplicate_reads)
我当时咋算的？
在数据准确有效的情况下，每检测到一种独特的reads，该项目的reads类型计数增加1，N_reads表示该项目共检测到了N种独特的reads。n_deduped_reads表示在N种独特的reads之中，有n种reads仅被检测到了1次。测序饱和度是指至少被检测到2次的reads占比，也就是1 - (n_deduped_reads / N_reads)
1减去唯一reads占比就是饱和度：
1 - (12435096 / 24451006)
[1] 0.491428
直接重复reads占比就是饱和度：
6328497 / 24451006
[1] 0.2588236
其实这两种方法都没有错，错的是分母：N_reads表示该项目共检测到了N种独特的reads
公式计算中所有的reads都应该是独特的reads来进行计算，至于是唯一还是重复都只是这个独特read的属性
如（A、A、B、B、B、C）中独特reads（A[重复]、B[重复]、C[唯一]），计算饱和度就应该是：2/3，而不是5/6
所以这也是我们在解读samtools flagstat的汇总信息时需要注意的地方：
Total指的是所有reads数，并不是所有独特reads数，作为参与计算饱和度的所有独特reads数
评论区老师提到官网的公式更加清晰：unique_confidently_reads+duplicate_reads
而一文带你读懂 10X Cellranger Count 网页结果解析公式关于N_reads需要理解他前面说的“N_reads表示该项目共检测到了N种独特的reads”
duplicates则是所有独特重复reads数，并不是所有重复reads
这两个观测并不在一个level上
统计序列完全一致的reads的频率，横轴表示重复水平，纵轴表示重复⽔平序列列占所有序列的百分比
duplicate是全部序列的duplicate的情况吗？还是随机筛选了一部分？为什什么要这样做？
是选择的每一个⽂文件里前100,000条序列作为样本进行的计算，因为样本本身很⼤，前100,000已经能够代表样 本的重复性
存在两条线，一蓝一红
蓝线“% Total Sequence”很好理解，就是测到的所有reads，红线“% Deduplicated sequences”就需要联系我们上面谈到的单细胞cellranger输出报告中的饱和度来理解，纳入计算的是独特的reads
根据上面草图中举的例子可以看到，在Sequence Duplication Level为1时，因为这个时候所有reads都是唯一的，唯一reads总数就等于独特reads总数，但分母肯定是“% Deduplicated sequences”要小，因为只是所有独特的reads（每种reads只有一个做代表），这也就解释了为什么图中红线在Sequence Duplication Level为1时要比蓝线高
但是随着Sequence Duplication Level的增加，“% Total Sequence”比“% Deduplicated sequences”保留了更多非独特reads，这也就解释了为什么后面会出现蓝线高于红线
这里的蓝线后面出现小峰是因为横坐标把多个levels合在一起了没所以扩大了视觉上的效果
综上，经过上面两个案例的学习，我们似乎可以得到这样的一个tip：在转录组上游的统计信息中，total字样一般指所有reads，而unique、deduplicated、duplicated字样往往只统计独特的reads，这两个并不在同一level
这其实类似R语言中的unique函数， 如
虽然使用unique函数，这里独特reads变量名使用spec是为了区别于前文里的unique reads：只被检测到一次的reads，在这个小例子中unique_reads就只有“C”
上面两个案例都是从测序reads饱和角度（唯一/sequence duplication level-1 vs 重复/多个levels）出发，对于不同level的数据，曾老师给我分享了信息熵的4个量化指标的R代码实现 ，大家也可以学习看看
联系客服
