XHMM (Fromer et al., 2012) は CNV を検出するソフトウェアの一つである。XHMM は、WES のデータを GATK などで処理してから得られる depth of coverage file を入力データとして、CNV を検出する。詳細な XHMM の使い方は XHMM ウェブサイトあるいは Fromer et al., 2014 の論文で確認出来る 。

XHMM の使い方

データの準備

XHMM では、WES データをマッピングした結果である BAM ファイルを GATK で処理した後のデータを入力データとしている。サンプルデータ（EXAMPLE_BAMS.zip）は XHMM ウェブサイトでダウンロードできる。また、マッピング時に利用した FASTA ファイルは 1000 genome プロジェクトのサイトでダウンロードできる。

• EXAMPLE_BAMS.zip
  展開後に BAM ファイルと BAM ファイルのインデックス、EXOME.interval_list、group1.READS.bam.list などのファイルが見られる。
• human_g1k_v37.fasta.gz
• human_g1k_v37.fasta.fai

GATK DepthOfCoverage 処理

BAM ファイルに対して、GATK DepthOfCoverage 処理を行い、sequencing depth を計算する。BAM ファイルが多い場合は、分割して GATK を実行することができる。この例の場合、サンプルデータの中には BAM ファイルが 30 個あり、これを 3 つのグループに分けて GATK の処理を行う。各グループにどのような BAM が入っているのかは group*.READS.bam.list のファイルに記載されている。

java -Xmx3072m -jar GenomeAnalysisTK.jar \
     -T DepthOfCoverage -I group1.READS.bam.list -L EXOME.interval_list \
     -R ./human_g1k_v37.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 group1.DATA

java -Xmx3072m -jar GenomeAnalysisTK.jar \
     -T DepthOfCoverage -I group2.READS.bam.list -L EXOME.interval_list \
     -R ./human_g1k_v37.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 group2.DATA

java -Xmx3072m -jar GenomeAnalysisTK.jar \
     -T DepthOfCoverage -I group3.READS.bam.list -L EXOME.interval_list \
     -R ./human_g1k_v37.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 group3.DATA

GATK を実行し終えた後に、3 つの結果をマージする。

xhmm --mergeGATKdepths -o ./DATA.RD.txt \
     --GATKdepths group1.DATA.sample_interval_summary \
     --GATKdepths group2.DATA.sample_interval_summary \
     --GATKdepths group3.DATA.sample_interval_summary

XHMM 処理

XHMM を用いた処理では、

1. 外れ値を除去
2. PCA
3. PCA 結果を用いて CNV を検出

の順で行う。

まずは sequencing depth に外れ値が存在するサンプルや検出ターゲット領域を除去する。

フィルタリングの 1 つである GC 含量を計算する。

java -Xmx3072m -jar GenomeAnalysisTK.jar \
    -T GCContentByInterval -L EXOME.interval_list \
    -R ./human_g1k_v37.fasta \
    -o ./DATA.locus_GC.txt

cat ./DATA.locus_GC.txt | awk '{if ($2 < 0.1 || $2 > 0.9) print $1}' > ./extreme_gc_targets.txt

次に、フィルタリングの 1 つである repeat-masked bases を計算する。

./sources/scripts/interval_list_to_pseq_reg EXOME.interval_list > ./EXOME.targets.reg

pseq . loc-load --locdb ./EXOME.targets.LOCDB --file ./EXOME.targets.reg --group targets --out ./EXOME.targets.LOCDB.loc-load

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

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

フィルタリング＆外れ値除外を行う。

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

次に、外れ値を除去したのちに、PCA を行う。

xhmm --PCA \
     -r ./DATA.filtered_centered.RD.txt \
     --PCAfiles ./DATA.RD_PCA

PCA を行ったのちに、分散の大きい主成分を取り除く。これにより、CNV に由来しない主成分を除去することができる。具体的にある主成分の分散が、全主成分の分散よりも 0.7 倍よりも大きければ、その主成分を取り除く。そのような主成分を取り除いてから、残りのデータをもう一度正規化する。

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

正規化に失敗したいくつかのデータをさらに取り除くために、Z スコアによるフィルタリングを行う。

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

ここまでの手順において、入力データに対して、主成分の分散や Z スコアによるフィルタリングが行われた。この過程において、いくつかのサンプルとターゲット領域が取り除かれた。ここで得られた解析後のデータと最初の入力データ DATA.RD.txt を比較可能にするために、同様なサンプルとターゲット領域を DATA.RD.txt からも取り除く。

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

最後に HMM のアルゴリズムを利用して CNV を検出する。

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

出力ファイルである DATA.xcnv には CNV の検出結果が保存されている。1 つの CNV につき 1 行で記載されている。

上で検出できた CNV に対して、統計的な手法により検出結果の評価を行う。

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

xcnv フォーマット

XHMM は最終的に拡張子が xcnv のタブ区切りのテキストファイルを出力する。

SAMPLE	CNV	INTERVAL	KB	CHR	MID_BP	TARGETS	NUM_TARG	Q_EXACT	Q_SOME	Q_NON_DIPLOID	Q_START	Q_STOP	MEAN_RD	MEAN_ORIG_RD
HG00121	DEL	22:18898402-18913235	14.83	22	18905818	104..117	14	9	94	94	84	-2.54	36.73
HG00113	DUP	22:17071768-17073440	1.67	22	17072604	4..11	8	25	99	99	54	25	4.04	193.79