集団遺伝学

PSMCは、1個体のゲノムから集団の人口動態の変遷を推定する集団遺伝学のプログラムです。導入にあたり、日本語の説明が見当たらなかったので覚書をしておきます (2015年)。

PSMC

ダウンロード:https://github.com/lh3/psmc
今回は、version t 0.6.5-r67を使っています。

PSMCの作業工程は、以下のパートに分かれます。
  1. samtools mpileupを利用した変異検出
  2. 独自フォーマットの変換(psmcfa)
  3. PSMCの実行
  4. グラフの描写

1. PSMCではBAMファイルからsamtools mpileupを使って変異サイトをコールし、fastqに落としたものを使う必要がありますが、GATKでcallingしたSNP (vcf)を使うことは出来ませんでしたので、おとなしく、指示通りにsamtools mpileupを使いましょう。samtools mpileupにおけるdepthの上限(-D)と下限(-d)の設定は、公式ドキュメントによると、-Dがcoverage平均の2倍、-dが1/3が推奨されています。またlow coverage dataからの推定はお勧めしないとのことです。また0.1.19以下のsamtoolsでないとうまくいかないという話題が海外の掲示板に出ていましたのでご注意。

2. PSMCの入力ファイルは、psmcfaというファイル形式で、そのファイルコンバーターはPSMCのパッケージに含まれています( utils/fq2psmcfa )。中身はfasta-likeなファイル形式でIUPACで定められた塩基配列表記を使っているようです。出来上がる配列ファイルは小さいので不安になりますがw、下記の例で用いているHGDPのゲノムは31Mバイトくらいのファイルサイズでした。

ヒト全ゲノムを解析している印象としては、mpileup->PSMCの結果を出力するまでに合計で約1.5〜2日かかっています(ゲノム解析用スパコンを利用した場合)。mpileupに多大な時間を費やしますが、PSMC自体は数時間以内で終わります。

Low coverage dataでは何が起こるか?


low coverage data の代表である、1000Genomes Projectからランダムに日本人(JPT)ゲノムを1個体選んで、他の論文から引っ張ってきた25-30xのゲノムと比較したのが下図です。なお、Phase 3およびPCR-freeシーケンスデータを使ったPSMCの結果が、1000Genomes Projectの新しい論文1のExtended Figure 7に記載されています。

比較に用いたCHB, French, YorubaはいずれもHuman Genome Diversity Project由来の試料で、MPI-EVAの論文Meyer et al. 2012において、イルミナ社にてシーケンスされたものです。一方、1000Genomesのphase1の平均coverageは5〜6x程度ですので、推奨のdepthの下限は1/3、と言われても、無駄な抵抗とは思いつつ、2種類試してみました。

利用したsamtoolsは0.1.18です。
reference は変数 $refseqで指定、今回はHG19を利用しています。

(JPTパターン1) -d 3 -D 18
samtools mpileup -C50 -uf $refseq 1KG_JPT_low.bam |bcftools view -c -| vcfutils.pl vcf2fq -d 3 -D 18 | gzip > 1KG_JPT_low.d3D18.fq.gz

(JPTパターン2) -d 6 -D 50
samtools mpileup -C50 -uf $refseq 1KG_JPT_low.bam |bcftools view -c -| vcfutils.pl vcf2fq -d 6 -D 50 | gzip > 1KG_JPT_low.d6D50.fq.gz

HGDPのデータは下記に統一しました
samtools mpileup -C50 -uf $refseq HGDP_CHB.bam |bcftools view -c -| vcfutils.pl vcf2fq -d 8 -D 50 | gzip > HGDP_CHB.d8D50.fq.gz

コレ以降は、officialな例に従いました(具体的なコードは省略)。

utils/fq2psmcfa -q20 diploid.fq.gz > diploid.psmcfa psmc -N25 -t15 -r5 -p "4+25*2+4+6" -o diploid.psmc diploid.psmcfa utils/psmc2history.pl diploid.psmc | utils/history2ms.pl > ms-cmd.sh utils/psmc_plot.pl diploid diploid.psmc


low coverage genomeに対しては、plotの際にFalse Negative Ratio(FNR)を設定できるとのことで、1000Genomesでは確か20%と議論していた気がしたので、plotの際にはFNRを振って様子を見てみることにします。psmc_plot.plのオプション -M で各データのラベリングとFNRを設定できます。

日本人ゲノムだけで比較してみます。最初の4つのラベル1KG-JPTはパターン1で、FNRを0, 0.1, 0.2, 0.3 と振ってました。最後の1KG-JPT-deepがパターン2です。
psmc_plot.pl -G -P "right bottom" -M "1KG-JPT,1KG-JPT=0.1,1KG-JPT=0.2,1KG-JPT=0.3,1KG-JPT-deep" 1kg-low-JPT.2.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d6D50.psmc

結果がこちら
JPT low coverage

あってないようなdepthを使っているパターン1のFNR10%と、それよりかはdepthを少し厳しくしたパターン2がほぼ同じような推定値を示しているように見えます(水色と緑)。推定が106 年前までとなりました。105年前から最近のbottleneckのあたりが大きくぶれるようです。

今度は、HGDPゲノムと日本人ゲノムを一緒にしてみます。
psmc_plot.pl -G -P "right bottom" -M "CHB,French,Yoruba,1KG-JPT,1KG-JPT=0.1,1KG-JPT=0.2,1KG-JPT=0.3,1KG-JPT-deep" 1kg-low-JPT_middleHGDP.psmc HGDP00778.psmc HGDP00521.psmc HGDP00927.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d3D18.psmc 1kg_JPT.d6D50.psmc



やっぱり変ですね・・・。


参考文献
  1. The 1000 Geonmes Project Consortium, "A global reference for human genetic variation", Nature 526, 68–74 (01 October 2015) . Link



サブページ (1): 23andMe
Comments