数据模拟软件 VarBen

随着NGS在分子诊断中变异(体系和胚系)检测的应用普及。临床NGS检测产品的检测范围越来越大,带来的一个优势是可以识别更多基因中的变异;但是同时也使性能验证面临巨大的挑战。临床方法学验证过程中很难获得真实样本,能同时具有大量关注且能够被检测到频率也合适的变异。所以大多数实验室都是通过对具有已知变异的标准品或细胞系进行测序分析,并依靠测序指标(即具有足够覆盖率的目标部分)来推断其他区域的性能。

在这个大背景下,很多实验室会开始考虑使用模拟数据进行性能评估,使用数据模拟的方法在临床实验室中具有非常高的实用性,可作为临床检测验证的辅助工具。通过模拟数据可以便捷的构建出真实样本中难以获取的变异数据,从而使实验室能够更高效经济同时准确地测试生物信息学流程的性能(不同频率梯度、流程灵敏度,针对不同长度插入/缺失的性能边界)。从而进行更全面的方法学验证和生信流程验证,同时模拟数据也已经成为AMP的建议方法,同时针对不同的检测目的提供了不同的模拟方案,可以参考文章。同时文章中也提供了非常多的数据模拟软件和方案,可以参考模拟软件参考列表

而今天我们介绍的是VarBen,一款由中检院开发的可以模拟4中变异(SNV、InDel、CNV、SV)的软件。也许整体性能和使用便捷性上不是最好的,但是中检院背书,在国内资质申报等注册相关渠道来说,无疑是最容易被接受的。

安装

软件是基于Python2开发的,借助python重点pysam(>=0.9.4)模块进行对Bam文件进行处理,对比对到参考基因组上的数据进行编辑制造预期的变异。在完成对reads进行编辑的操作以后,同样会对reads进行重新的比对,所以环境中需要提前进行bwa和samtools的安装。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
python2 -m pip install pysam
python2 -m pip install numpy

# samtools (http://samtools.sourceforge.net/)
git clone https://github.com/samtools/htslib.git
make -C htslib

git clone https://github.com/samtools/samtools.git
make -C samtools
cp samtools/samtools $HOME/bin

git clone https://github.com/samtools/bcftools.git
make -C bcftools
cp bcftools/bcftools $HOME/bin


# bwa (http://bio-bwa.sourceforge.net/)
git clone https://github.com/lh3/bwa.git
make -C bwa
cp bwa/bwa $HOME/bin

使用方法

该软件支持多种不同类型的变异模拟。同时也支持多种不同的测序平台(Illumina、life 、MGI),根据要模拟的不同变异类型,需要准备不同的配置文件进行模拟。
因为仓库有相对系统的说明,在这里不进行展开,主要介绍一些容易被忽视的细节和局限性。

只处理Uniq Reads

NGS数据处理过程中,Markdup已经成为流程处理的一个标准化步骤,不管使用什么方法,这步其实都是难以避免的。
如果使用的是类似Picard的方法,进行Markdup处理(Bam数据中存在Duplication Reads)那么需要注意,数据模拟只会对Uniq Reads进行处理。实例如下图:
image
所以如果你的下游处理流程如果会对Dup簇信息进行整体的评估(例如评估变异的可信度),那么需要仔细的评估确认你检测流程是否适合使用这种方法进行模拟。
当然一个比较投机的方案是可以使用一个没有Markdup的Bam进行模拟(但是分析流程也无法进行模拟操作),比如在处理前,先把Bam中的Duplication Reads进行删除,使用rmDup而不是Markdup。所以对应的如果我们处理的是压缩的方法构建consensus Reads(数据中没有Duplication Reads)那就没什么影响,

Reads名称变化

在进行SV等类型模时,软件会对模拟Bam中的reads名称进行变更,可以在临时目录(tempDir/readname_convert.txt)查看进行reads名称变更(模拟过程中对reads进行编辑)的reads。示例如下:

1
2
3
4
5
6
7
8
9
10
(base)[Ben@ECS1 tempDir]$head tempDir/readname_convert.txt
dd10d900-02c6-48f5-8f71-10ea8844fea5: F350021481L2C003R05000360198_TGTCTGTGTCTG_TGTCTGTGTCTG, True, 43611912-43612005
dd10d900-02c6-48f5-8f71-10ea8844fea5: F350021481L2C003R05000360198_TGTCTGTGTCTG_TGTCTGTGTCTG, False, 43612017-43612110
8e77a7cc-baa4-4e2e-9f78-bb70ae3b53ff: F350021481L1C001R06500913354_TGTCTGTGTCTG_TGTCTGTGTCTG, True, 43611800-43611893
8e77a7cc-baa4-4e2e-9f78-bb70ae3b53ff: F350021481L1C001R06500913354_TGTCTGTGTCTG_TGTCTGTGTCTG, False, 43611976-43612069
3e37c602-85db-4018-8a6a-daaf252839af: F350021481L1C003R03701364230_TGTCTGTGTCTG_TGTCTGTGTCTG, True, 43611904-43611997
3e37c602-85db-4018-8a6a-daaf252839af: F350021481L1C003R03701364230_TGTCTGTGTCTG_TGTCTGTGTCTG, False, 43611953-43612046
09e73b73-c5c5-4c42-956a-61b42b0cab8c: F350021481L2C002R01000826519_TGTCTGTGTCTG_TGTCTGTGTCTG, False, 43611794-43611887
09e73b73-c5c5-4c42-956a-61b42b0cab8c: F350021481L2C002R01000826519_TGTCTGTGTCTG_TGTCTGTGTCTG, True, 43611952-43612045
4020f12a-3e66-42c6-bbd2-13ce87783f7d: F350021481L2C001R02700414361_TGTCTGACAGAC_TGTCTGACAGAC, False, 43611921-43612014

其实我们可以看到,前面第一列是Read更新后的ID,第二列是reads更新前的原始ID,第三列疑似正负链信息,第四列是该条Reads在染色体上的位置。
因为处理后支持融合的ReadsID名称及ID格式都发生了变化,所以在使用中需要确认这个ID变化对你下游检测流程的影响。
比如示例的数据如果想从Bam生成对应的Fastq测试完整流程显然是不可行的,因为ReadsID中的UMI信息由于ID变化已经完全丢失了,这时候,就需要基于上述ReadsID对应文件,将ReadsID逐一进行数据的还原。
相反,如果直接基于Bam进行检测,下游分析对ReadID格式不敏感,那就不需要考虑这部分影响。

CNV

CNV检测,目前常规检测方法都是基于深度的,而模拟过程其实就是基于Bam提取想要模拟CNV区域的Reads进行double,所以针对CNV,模拟后的数据不适合进行再次的Markdup,因为重新Markdup,会直接导致你的模拟数据都被剔除了。

-------------本文结束感谢您的阅读-------------