测序饱和度【cellranger】

前言

  单细胞转录组测序技术逐渐的，就像父辈人年轻时，彰显身份的“大哥大”，变成了我们这一代人，人手一部的现代手机，成了很多文章中的“头牌”或是“锦上添花”的部分，就像普通的转录组测序的发展趋势。
  有时候，为了能在领域内”活下去“，我们不得不拼命的追赶潮流😑

什么是测序饱和度

  常用的单细胞转录组定量工具：cellranger count，其Web Summary中会有一个叫测序饱和度(Sequencing Saturation) 的概念（图示红框）

引用10X官网的解释，测序饱和度就是：

具体地说，它表示非唯一的可靠比对、有效的细胞条形码和有效的UMI Reads 所占的比例（与现有的细胞Barcode、UMI和基因组合匹配）。

  它是一个比例，反映当前测序量与文库复杂度相关性的指标，其大小主要取决于测序深度和文库复杂度。如何考量？
  一般来说，测序Reads越多，被检测到是独特转录本的频次就越大。被检测到的基因数会随着测序深度的增加而增加，但当测序深度达到一定程度，被检测到基因数量的增加程度逐渐变缓，直到不再随测序深度增加而增加，即达到饱和。测序饱和度可以一定程度上反应测序的深度是否足够。

  不同类型的细胞通常有不同数量、不同类型的转录本，因此不同类型细胞建成的文库复杂度存在差异。复杂度高的文库中，检测一个新转录本所需的额外Reads更多，即需要测更多的Reads增加饱和度。
  那么，多少的测序饱和度是合适的、足够的？ 我不能给出确切的数值，只能说，要根据实验材料和需求等，”因地制宜“，综合考量测序饱和度、基因中个位数等信息，比如，拟南芥的单细胞转录组测序，一个样本测了100G，拿到了10k左右细胞，基因中位数有3.5k多，测序饱和度有30%，我会认为饱和度可能够了，这个数据很好，因为拟南芥大概有2.5w个基因，同样的，如果换成了水稻的单细胞转录组，一个样本同样测了100G，10k左右细胞，基因中位数有3.5k左右，测序饱和度有30%，我会认为，这个数据还不错，因为水稻的基因大概有3.5w，比拟南芥要多出很多，这个时候，我可能会考虑，是不是要再加测50G(提高测序深度)来提高测序饱和度，获得更多信息。
  当我拿到一份数据的时候，往往会先分析下去，然后再回头考量数据的好坏。

如何计算

最后，聊一聊对于我来说，可能不那么重要的计算方式：
$$SequencingSaturation=1-\frac{n_{dedupedReads}}{n_{Reads}}$$
$n_{dedupedReads}$：可靠比对的Reads中，唯一的（有效Barcode、UMI、基因）组合数量
$n_{Reads}$：所有的可靠比对、有效Barcode、UMI的总Reads
$n_{Reads}=UniqueConfidentlyMappedReads+(DuplicateReads-UnmappedDuplicates)$

一个例子

10X的示例数据： https://support.10xgenomics.com/single-cell-gene-expression/datasets/3.0.0/pbmc_1k_v3

其 Web Summary 中显示测序饱和度为：0.7085123，计算数据及过程大概为：

unique_confidently_mapped_reads = 10,196,940
duplicate_reads = 24,785,461

x = 1 - (unique_confidently_mapped_reads/(unique_confidently_mapped_reads + duplicate_reads))
x = 1 - (10,196,940/(10,196,940 + 24,785,461))
x = 1 - (10,196,940/34,982,401)
x = 1 - 0.29148771
x = 0.70851229

那这些数据怎么来的？
  首先，unique_confidently_mapped_reads，就是前面翻译成唯一可靠比对Reads的东西……在BAM文件中会被 xf tag 标记值为25（1 + 8 + 16），具体的含义可以参考10X output，可以使用samtools查看：

samtools view pbmc_possorted_genome_bam.bam | grep 'xf:i:25' | wc -l 

  然后，duplicate_reads会被SAM/BAM flag 标记为 0X400，可以通过samtools flagstat获取：

samtools flagstat pbmc_possorted_genome_bam.bam
76920923 + 0 in total (QC-passed reads + QC-failed reads)
10319036 + 0 secondary
0 + 0 supplementary
24785461 + 0 duplicates
73840063 + 0 mapped (95.99% : N/A)
...

  之后，计算unmapped_duplicates
如 10X output 中最后一个表格所述，Feature Barcode Read都是Unmapped。在这些Reads中，有一些会被标记为duplicate，此时SAM flag为1028。可是使用samtools view获取：

samtools PBMC_possorted_genome_bam.bam | grep 'fb:Z:' | cut -f2 | grep '1028' | wc -l

之后，就是把数字带进公式，计算出测序饱和度……

参考：

1. https://kb.10xgenomics.com/hc/en-us/articles/115003646912-How-is-sequencing-saturation-calculated-
2. https://kb.10xgenomics.com/hc/en-us/articles/115005062366