https://academic.oup.com/nar/article/44/20/e154/2607979
该软件的源代码保存在sourceforge上,链接如下
https://sourceforge.net/projects/excavator2tool/
excavator2在计算测序深度时将reads分为了以下两个部分
in-target reads
off-target reads
in-target表示的是位于exon上的序列,off-target表示的是位于基因间区或者内含子区的序列,同样采用了滑动窗口的方式来统计每个区域的测序深度,只不过稍作变化,全称如下
mean windows read count
简称WMRC, 计算公式如下
单个外显子直接作为一个窗口,而非外显子区域则采用了一个固定长度的窗口,分开统计不同区域的测序深度,并进行校正,校正的时候考虑了GC含量,不同区域的mappability, 外显子的大小等因素。
利用归一化之后的测序深度,计算对照样本和实验样本的log2 ratio值,然后采用HSLM segmentation算法划分segment, 最后通过FastCall算法预测每个segment的拷贝数情况,细分为以下5种
two-copy deletion
one-copy deletion
normal
one-copy duplication
multiple-copy duplicaiton
该软件支持hg19和hg38两个版本,内置了对应的数据库,示意如下
该软件分成了3个模块,对应3个脚本,具体操作步骤如下
1. TargetPerla.pl提供一个捕获区域的bed文件,计算in-target和off-target区域的GC含量,mappability值,用于后续的归一化操作,用法如下
perl TargetPerla.pl SourceTarget.txt myTarget.bed MyTarget_w50000 50000 hg19
第一个参数为source target文件,该文件记录了基因组对应的bw文件和fasta文件的路径,内容示意如下
/data/ucsc.hg19.bw /data/hg19.fasta
空格分隔的两列,第一列为bw文件的路径,该文件是软件自带的,位于软件的安装目录,用于计算基因组不同区域的mappability,第二列为fasta文件的路径,用于计算不同区域的GC含量。
第二个参数为捕获区域的bed文件,第三个参数为输出结果的前缀,第四个参数为窗口的固定长度,第五个参数指定基因组的版本。
这一步类似于比对时建立参考基因组的索引,一个芯片建立一次即可,运行成功后,会生成一个文件夹,前缀为MyTarget_w50000。
2. EXCAVATORDataPrepare.pl计算测序深度,进行归一化处理,用法如下
perl EXCAVATORDataPrepare.pl ExperimentalFile.txt --processors 6 --target MyTarget_w50000 --assembly hg19
第一个参数是一个空格分隔的txt文件,指定了样本对应的bam文件,输出结果的路径,样本名称信息,内容示意如下
--processors
指定并行的线程数,--target
参数指定第一步生成的target的名称,--assembly
指定参考基因组的版本。
执行HSLM segmentation算法和FastCall算法,进行CNV分析,用法如下
perl EXCAVATORDataAnalysis.pl ExperimentalFileAnalysis.txt --processors 6 --target MyTarget_w50000 --assembly hg19 --output Results_MyProject_w50K --mode pooling
--mode
参数表示样本如何进行比较,支持pooling和paired两种模式,第一种模式将所有的实验样本混合与对照样本进行比较,第二种模式则是配对样本模式,比如癌和癌旁,两两之间进行比较,计算log2 ration值。
第一个参数ExperimentalFileAnalysis.txt是一个空格分隔的txt文件,指定了样本的比较操作,对于pooling模式,其内容示意如下
对于paired模式,其内容示意如下
T表示Treat, C表示Control, 后面的数字用于区分不同样本。
--output
指定了输出结果的目录,输出结果中提供了CNV区域对应的txt, VCF等文件,同时还提供了可视化的结果,示意如下
·end·