集団構造解析を行うための準備

ここでは、集団構造解析を行うにあたり必要な準備について解説します。主に解析用の遺伝子型ファイル（VCFファイル）のフィルタリングの手順について述べます。

Contents

1. フィルタリング用ツールのインストール
2. VCFファイルのフィルタリング
3. References

1. フィルタリング用ツールのインストール

1.1. VCFtoolsのインストール

GitHub - vcftools/vcftools: A set of tools written in Perl and C++ for working with VCF files, such as those generated by the 1000 Genomes Project.

A set of tools written in Perl and C++ for working with VCF files, such as those generated by the 1000 Genomes Project. - GitHub - vcftools/vcftools: A set of ...

git clone https://github.com/vcftools/vcftools.git
cd vcftools
./autogen.sh
./configure
make
make install

1.2. BCFtoolsのインストール

SAMtools/BCFtools/HTSlib - Downloads

Samtools

wget https://github.com/samtools/bcftools/releases/download/1.16/bcftools-1.16.tar.bz2
tar xvf bcftools-1.16.tar.bz2
cd bcftools-1.x 
./configure --prefix=/where/to/install
make
make install

1.3. PLINKのダウンロード

PLINK 1.9

mkdir plink
cd plink
wget https://s3.amazonaws.com/plink1-assets/plink_linux_x86_64_20220402.zip      #Linuxの場合
unzip plink_linux_x86_64_20220402.zip

それぞれ、インストール先のディレクトリにパスを通しておいてください。

1.4. Beagleのダウンロード

Beagle 5.4

VCFファイルのインピュテーションツールです。これ一つでgenotype imputation（欠測した遺伝子型を統計学的に推定）からhaplotype phasing（ハプロタイプの推定、selective sweep解析で必要になります）まで実行してくれます。ただし、穴ボコだらけのVCFでもとりあえず穴埋めしてしまうため、使用前にあまりにも信頼性の低い座位は除外しておく必要があります。

wget https://faculty.washington.edu/browning/beagle/beagle.22Jul22.46e.jar

インストールは必要ありません。ダウンロードしたファイルをパスが通ったディレクトリに移しておきましょう。

2. VCFファイルのフィルタリング

どの集団構造解析（PCA、系統解析、ADMIXTURE解析）にも言えることですが、「信頼性が高く」、「情報のダブり（重複）がない」SNPsをゲノムワイドに満遍なく配置したVCFファイルを準備することが重要です。また、PCAはその性質上、missing data（ある座位において、genotypingできている個体とできていない個体がある状態、つまりデータに穴がある状態）を許さないので、そのような座位は事前に解析から除外するか、統計的に推定して穴埋め（imputation）しておく必要があります。筆者がウメの集団構造解析で踏んだ手順は以下の通りです。

デプスが基準を下回るSNPsを除外する
Genotyping率が基準を下回るSNPsを除外する
VCFファイルのimputationを行う
Biallelic SNPs以外を除外する
マイナーアリル頻度が基準を下回るSNPsを除外
ハーディ・ワインベルグ平衡から著しく逸脱するSNPsを除外
強い連鎖不平衡状態にあるSNPsは間引く

それぞれについて具体的に解説していきます。

2.1. デプスが基準を下回るSNPsを除外する

「デプス」とは、簡単に言うと、「そのSNPの確からしさが何本のリードによって保証されているか」を表す指標です。例えば、ある座位が「A」と「G」のアリルで成り立っていたとして、その塩基を保証するリードの数が1本だったらどうでしょう？ちょっと心許ないですよね。ここでは、最小デプスを8を下回るSNPsは除外しておきます。なお、SNP callingにVarScan2を用いた場合は、デフォルトで最小デプス8のSNPsが選抜されていますので、この処理は不要です。

処理前のVCFファイルの名称を「before_filtering.vcf」とします。

vcftools \
--vcf before_filtering.vcf \
--minDP 8 \          #最小デプスを8に設定
--recode \           # vcfファイルを出力
--out minDP8_applied          # minDP_applyed.recode.vcfというファイルが出力されます

2.2. Genotyping率が基準を下回るSNPsを除外する

あまりにもgenotyping率の低いSNPsはそもそも信用できません。ただし、100%を目指そうとするとほとんど除外されて無くなってしまいますので、ある程度の妥協は必要です。ここではgenotyping率8割以上のSNPsを残すことにします。

vcftools \
--vcf minDP8_applied.recode.vcf \
--max-missing 0.8                    # 1-(missing rate)の値を入力
--recode
--out minDP8_mm0.2_applied

2.3. VCFファイルのimputationを行う

次に、Beagleを用いてVCFファイルの欠測値を埋めていきます。あくまでも「確率的に推定される遺伝子型」を当てはめていく作業であるため、少しくらいの間違いはあるということを頭の片隅に置いておきましょう。

java -jar beagle.22Jul22.46e.jar \
gt=minDP8_mm0.2_applied.recode.vcf \
out=minDP8_mm0.2_imputed \
nthreads=2 　　　　　　　　　　# 使用CPUコア数　あまり多くすると解析が走らないので注意

imputation後のVCFファイルは「minDP_mm0.2_imputed.vcf.gz」という名前で出力されます。

2.4. Biallelic SNPs以外を除外する

集団構造解析を行うためのソフトウェアには、Biallelic（アリルが2種類）の座位しか読み込めないものが多くあります。そこで、予めアリルが3つ以上ある座位を除外しておきます。操作はとても簡単で、VCFファイルをPLINK形式（pedファイルとmapファイル）に変換するだけでOKです。PLINKそのものがbiallelic SNPにしか対応していないので、変換時に自動的にそれ以外の座位を除外してくれます。

vcftools \
--gzvcf minDP8_mm0.2_imputed.vcf.gz \        #gzip圧縮ファイルを扱う場合は--gzvcfを指定
--plink                 　　　　　　　　　　#PLINK形式に変換
--out minDP8_mm0.2_imputed

2.5. マイナーアリル頻度が基準を下回るSNPsを除外

マイナーアリル頻度（Minor Allele Frequency; MAF）とは、あるbiallelicの座位において、「頻度の少ない方のアリルの頻度」のことです。例えば下表のような遺伝子型の表があったとそれば、SNP1座では100サンプルのうち、95サンプルがAアリル、5サンプルがGアリルを持っているので、MAF = 5/100 =0.05となります。サンプルサイズがどのくらいあるかにもよりますが、あまりにもMAFが小さいSNPは、シークエンシングエラー由来の可能性が否定できないため、MAF > 0.03～0.05程度でフィルタリングすることが多いです。

	SNP1	SNP2	SNP3	SNP4	SNP5
Sample 1	A	A	A	A	A
Sample 2	A	A	A	A	A
・・・・・	・・・・・	・・・・・	・・・・・	・・・・・	・・・・・
・・・・・	・・・・・	・・・・・	・・・・・	・・・・・	・・・・・
・・・・・	・・・・・	・・・・・	・・・・・	・・・・・	・・・・・
Sample 96	G	A	A	A	A
Sample 97	G	G	A	A	A
Sample 98	G	G	G	A	A
Sample 99	G	G	G	G	A
Sample 100	G	G	G	G	G
MAF	0.05	0.04	0.03	0.02	0.01

ここでは、MAF > 0.03でフィルタリングを行います。

plink \
--file minDP8_mm0.2_imputed \
--maf 0.03 \
--recode \
--out minDP8_mm0.2_imputed_maf0.03

2.6. ハーディ・ワインベルグ平衡から著しく逸脱するSNPsを除外

ある集団における遺伝子頻度が、ハーディ・ワインベルグの法則を満たしているとき、ハーディ・ワインベルグ平衡（Hardy-Weinberg Equilibrium; HWE）にあるといいます。集団構造解析の手法は、基本的にはHWE状態にあるマーカー遺伝子座を用いることを前提としていますので、あまりにもそこから逸脱したマーカーは解析から除外します（ただし、解析の目的によっては、敢えて含める場合もあります）。ここでは、p < 0.0001となる遺伝子座（SNPs）を除外します。

plink \
--file minDP8_mm0.2_imputed_maf0.03 \
--hwe 0.0001 \
--recode \
--out minDP8_mm0.2_imputed_maf0.03hwe

2.7. 強い連鎖不平衡状態にあるSNPsは間引く

連鎖不平衡（Linkage Disequilibrium; LD）とは、集団において、複数の遺伝子座の対立遺伝子（多型）の間にランダムでない相関が見られる現象をいいます。例えば、下の図のように、ある遺伝子座の対立遺伝子が「T」であるとき、別の特定の遺伝子座の対立遺伝子は「C」になっていることが多い、というようなイメージです。このような遺伝領域は対立遺伝子の並び順が保存されており、「LDブロック」、「ハプロタイプブロック」などと呼ばれます。

少しわかりにくいかも知れませんが、「並び順が保存されている」ということは、そこにある多型はすべて「同じ情報を持っている」ということを意味します。同じLDブロック内にあるSNPは代表した1つに絞っておかないと、データに偏りが生じてしまいます。ですので、近傍（10 kbとか50 SNPs）にあるSNPs同士の相関関係（R²）を調べ、相関が基準値以上に高い場合は1つ間引くという処理を行います。

ここでは、50 SNPsの枠を3 SNPsずつずらしながらR²を計算し、R² > 0.5となった場合にSNPを1つ間引くという処理を行います。

plink \
--file minDP8_mm0.2_imputed_maf0.03hwe \
--indep 50 3 2 \           # 50 SNPsの枠を3 SNPsずつずらしながらR²を計算し、R² > 0.5となった場合にSNPを1つ間引く
--out indep_50_3_2

「indep_50_3_2.prune.in」というファイルに処理の結果残ったSNPsが、「indep_50_3_2.prune.out」というファイルに除外されたSNPsが出力されます。この結果を元の遺伝子型ファイルに適用し、LD pruning後の遺伝子型ファイルを出力します。

plink \
--file minDP8_mm0.2_imputed_maf0.03hwe \
--extract indep_50_3_2.prune.in \
--recode \
--out minDP8_mm0.2_imputed_maf0.03hwe_pruned

得られたファイルが全てのフィルタリングを終えたファイルです。必要に応じてVCFファイルに戻したり、bed形式に変換したりして解析に使用してください。

# VCFに戻す（SNPhyloに入力する場合）
plink \
--file minDP8_mm0.2_imputed_maf0.03hwe_pruned \
--recode vcf \
--out minDP8_mm0.2_imputed_maf0.03hwe_pruned

# bed形式に変換する（flushpca2、ADMIXTUREに入力する場合）
plink \
--file minDP8_mm0.2_imputed_maf0.03hwe_pruned \
--make-bed \
--out minDP8_mm0.2_imputed_maf0.03hwe_pruned

お疲れ様でした！

3. References

Petr Danecek, Adam Auton, Goncalo Abecasis, Cornelis A. Albers, Eric Banks, Mark A. DePristo, Robert E. Handsaker, Gerton Lunter, Gabor T. Marth, Stephen T. Sherry, Gilean McVean, Richard Durbin, 1000 Genomes Project Analysis Group, The variant call format and VCFtools, Bioinformatics, Volume 27, Issue 15, 1 August 2011, Pages 2156–2158, https://doi.org/10.1093/bioinformatics/btr330

Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011 Nov 1;27(21):2987-93. doi: 10.1093/bioinformatics/btr509. Epub 2011 Sep 8. PMID: 21903627; PMCID: PMC3198575.

Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, Sham PC. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007 Sep;81(3):559-75. doi: 10.1086/519795. Epub 2007 Jul 25. PMID: 17701901; PMCID: PMC1950838.

B L Browning, X Tian, Y Zhou, and S R Browning (2021) Fast two-stage phasing of large-scale sequence data. Am J Hum Genet 108(10):1880-1890. doi:10.1016/j.ajhg.2021.08.005

B L Browning, Y Zhou, and S R Browning (2018). A one-penny imputed genome from next generation reference panels. Am J Hum Genet 103(3):338-348. doi:10.1016/j.ajhg.2018.07.015

Follow me!