XHMM检测CNV

XHMM是一个用PCA降噪+HMM方法来检测全外显子CNV的软件。软件文档可在这里查看。

安装XHMM

XHMM安装依赖GCC 4.4以上，以及pthread和lapack。

wget https://bitbucket.org/statgen/xhmm/get/master.zip
unzip master.zip
cd statgen-xhmm-*
make

安装比较麻烦，可以到docker hub找一个现成的镜像。

下载测试数据

XHMM文档中从1000g弄了30个样本，每个样本取300个外显子作为测试数据

wget https://statgen.bitbucket.io/xhmm/EXAMPLE_BAMS.zip
unzip EXAMPLE_BAMS.zip

这些数据是使用hs37d5参考基因组进行比对，经过了排序、去重以及GATK BQSR的标准bam文件准备流程。

XHMM标准流程

获得覆盖深度

使用GATK3获得覆盖深度。由于GATK4的DepthOfCoverage与GATK3有较大差异，如--countType不能选择COUNT_FRAGMENTS等，这里还是使用GATK3进行。注意XHMM提供了一个interval_list，但我们日常使用还是习惯用bed格式多一点，而GATK的-L参数其实是能兼容bed的。

java -jar GenomeAnalysisTK.jar \
	-T DepthOfCoverage -I group1.bam.list -L bed/EXOME.bed \
	-R /b37/human_g1k_v37_decoy.fasta \
	-dt BY_SAMPLE -dcov 5000 -l INFO --omitDepthOutputAtEachBase --omitLocusTable \
	--minBaseQuality 0 --minMappingQuality 20 --start 1 --stop 5000 --nBins 200 \
	--includeRefNSites \
	--countType COUNT_FRAGMENTS \
	-o depth/group1

XHMM的教程中把bam分为了3组，来表明如果需要合并不同批次的数据应该怎么做。重复上面的流程获得3个批次的统计结果。也可以将list合并，一次获取结果。

合并结果

将GATK统计结果进行合并

xhmm --mergeGATKdepths \
	-o xhmm/data.rd.txt \
	--GATKdepths depth/group1.sample_interval_summary \
	--GATKdepths depth/group2.sample_interval_summary \
	--GATKdepths depth/group3.sample_interval_summary

获得target区域GC比例

使用GATK获得target区的GC比例，方便后面进行GC校正

java -jar GenomeAnalysisTK.jar \
	-T GCContentByInterval -L bed/EXOME.bed \
	-R /b37/human_g1k_v37_decoy.fasta \
	-o xhmm/data.locus.gc.txt
cat data.locus.gc.txt | awk '{if ($2 < 0.1 || $2 > 0.9) print $1}' > xhmm/extreme.gc.target.txt

计算重复碱基并过滤，可选

使用script文件夹中的interval_list_to_pseq_reg进行。这一步，在部分shell中会报错，可以直接看这个脚本的内容查看其修改的方式自己新建脚本进行修改。

interval_list_to_pseq_reg bed/EXOME.interval_list > bed/EXOME.target.reg

这里使用pseq进行，获取重复碱基多的区域并排除

pseq . loc-load \
	--locdb bed/EXOME.targets.LOCDB \
	--file bed/EXOME.target.reg \
	--group targets \
	--out bed/EXOME.targets.LOCDB.loc-load \
	--noweb

pseq . loc-stats \
	--locdb bed/EXOME.targets.LOCDB \
	--group targets --seqdb seqdb.hg19 | \
	awk '{if (NR > 1) print $_}' | \
	sort -k1 -g | awk '{print $10}' | \
	paste bed/EXOME.interval_list - | \
	awk '{print $1"\t"$2}' \
	> bed/DATA.locus_complexity.txt

cat bed/DATA.locus_complexity.txt | awk '{if ($2 > 0.25) print $1}' > bed/low_complexity_targets.txt

其中的seqdb.hg19可以在plinkseq官网下载，但是我发现官网下载都崩了，也可以自己构建

pseq . seq-load --seqdb seqdb.hg19 \
	--file hg19.fa.gz \
	--format build=hg19 repeat-mode=lower \
	--name hg19 --description UCSC-hg19

过滤

使用xhmm进行过滤

xhmm --matrix \
	-r xhmm/data.rd.txt --centerData --centerType target \
	-o xhmm/data.filtered_centered.RD.txt \
	--outputExcludedTargets xhmm/data.filtered_centered.RD.txt.filtered_targets.txt \
	--outputExcludedSamples xhmm/data.filtered_centered.RD.txt.filtered_samples.txt \
	--excludeTargets xhmm/extreme.gc.target.txt --excludeTargets bed/low_complexity_targets.txt \
	--minTargetSize 10 --maxTargetSize 10000 \
	--minMeanTargetRD 10 --maxMeanTargetRD 500 \
	--minMeanSampleRD 25 --maxMeanSampleRD 200 \
	--maxSdSampleRD 150

PCA聚类分析

使用xhmm进行PCA聚类

xhmm --PCA -r xhmm/data.filtered_centered.RD.txt --PCAfiles xhmm/data.RD.PCA

降噪

使用聚类结果进行降噪

xhmm --normalize -r xhmm/data.filtered_centered.RD.txt --PCAfiles xhmm/data.RD.PCA \
	--normalizeOutput xhmm/data.PCA_normalized.txt \
	--PCnormalizeMethod PVE_mean --PVE_mean_factor 0.7

Z检验

使用xhmm进行Z检验

xhmm --matrix \
	-r xhmm/data.PCA_normalized.txt --centerData --centerType sample --zScoreData \
	-o xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt \
	--outputExcludedTargets xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt.filtered_targets.txt \
	--outputExcludedSamples xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt.filtered_samples.txt \
	--maxSdTargetRD 30

原始深度过滤

使用xhmm进行过滤

xhmm --matrix \
	-r xhmm/data.rd.txt \
	--excludeTargets xhmm/data.filtered_centered.RD.txt.filtered_targets.txt \
	--excludeTargets xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt.filtered_targets.txt \
	--excludeSamples xhmm/data.filtered_centered.RD.txt.filtered_samples.txt \
	--excludeSamples xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt.filtered_samples.txt \
	-o xhmm/data.same_filtered.RD.txt

CNV检测

使用xhmm检测CNV检测，做到这一步其实在xhmm/data.xcnv文件里已经能看到结果了。

xhmm --discover \
	-p params.txt \
	-r xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt -R xhmm/data.same_filtered.RD.txt \
	-c xhmm/data.xcnv -a xhmm/data.aux_xcnv -s xhmm/data

我获得的结果某些值与XHMM标准答案有点差异，不过总体结果是相同的。

SAMPLE	CNV	INTERVAL	KB
HG00121	DEL	22:18898403-18913235	14.83
HG00113	DUP	22:17071769-17073440	1.67

结果Vcf导出

将结果转为VCF导出，通过看每个Sample的最后一个INFO是Y或N可以看出是否检出。

xhmm --genotype \
	-p params.txt \
	-r xhmm/data.PCA_normalized.filtered.sample_zscores.RD.txt -R xhmm/data.same_filtered.RD.txt \
	-g xhmm/data.xcnv -F /b37/human_g1k_v37_decoy.fasta \
	-v xhmm/data.xhmm.vcf

docker

使用docker

docker pull pzweuj/xhmm:1.0

写了一个脚本来运行

xhmmPipe.sh -i <bam dir> -o <output dir> -s <seqdb> -b <bed file> -r <ref.fa>