plinkを使ってゲノムの集団構造を体験する-パート1
plinkの使い方
導入難易度★☆☆☆☆
使用難易度★★★☆☆
この記事を読むと何ができるようになる?
今回から2回に分けて、plinkを使ったGWAS解析の演習を行います。まず、国際的なヒトゲノムプロジェクト(1000 Genome Project)により得られたシークエンスデータを公表しているサイトから、遺伝子型データと出身地データを取得します。次に、plinkを使って、集団間の遺伝的距離と地理的な距離を主成分分析(PCA)を行って解析します。最後に、結果をプロットした図をRを使って作成し、地域ごとにクラスターを形成していることを体験してみましょう。
今回は、解析に必要なデータの入手と主成分分析までを行います。
なお、本記事で扱う実習内容は「ゼロから実践する遺伝統計学セミナー」(羊土社, p128-135)の演習内容を参考にさせていただきました。この本は初めてGWAS解析を実施する人が知りたい内容が丁寧かつコンパクトに説明されており、他にも色々なplinkの使用方法、QTL解析やeQTL解析からマンハッタンプロットの描画まで、多くの演習内容が詰まっています。GWAS解析をこれから始めてみたい方には、非常におすすめです。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
1. ツールのインストール
conda install -c bioconda plink2
インストールが正しく導入されたか確認する。
plink2
正しく導入されていれば、下記が表示される。
PLINK v2.00a2.3 64-bit (24 Jan 2020) www.cog-genomics.org/plink/2.0/ (C) 2005-2020 Shaun Purcell, Christopher Chang GNU General Public License v3 plink2 <input flag(s)...> [command flag(s)...] [other flag(s)...] plink2 --help [flag name(s)...] Commands include --rm-dup list, --make-bpgen, --export, --freq, --geno-counts, --sample-counts, --missing, --hardy, --indep-pairwise, --ld, --sample-diff, --make-king, --king-cutoff, --write-samples, --write-snplist, --make-grm-list, --pca, --glm, --adjust-file, --score, --variant-score, --genotyping-rate, --pgen-info, --validate, and --zst-decompress. "plink2 --help | more" describes all functions.
2. 遺伝子型データと出身地データを入手する。
まず、Genome1000プロジェクトのサイトにアクセスして、必要なファイルをダウンロードする。
1000 Genomes | A Deep Catalog of Human Genetic Variation
Data > Download data from the IGSR FTP site > release > 20130502の順に移動する。 ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b・・・と続くファイルのURLをコピーする。
#そのままペーストすると http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz #となるので、httpの部分をftpに変更する。さらに、wgetコマンドを頭につける。(1.17GB) wget \ ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz
続いて、同じページの下部にある、integrated_call_samples_v3.20200731.ALL.pedをダウンロードする。
このファイルは小さいので、クリックしてダウンロードすればよい。
3. plinkの実行
まず、vcf形式をbed形式に変換する。
cd $workingdir plink --make-bed -vcf ALL.chr1.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz --recode --out ALL.chr1
すると、同じディレクトリー内に、ALL.chr1.bed, ALL.chr1.bim, ALL.chr1.famが生成される。
次に、ジェノタイプのサンプル数やSNP数の確認してみる。
これらのファイルに対して、wcコマンドで行数を確認すると、元のファイルに含まれるSNP数やサンプル数などの情報を簡単に知ることができる。*はワイルドカード。
wc ALL.chr1.*
1462893 3001925 4049026847 ALL.chr1.bed 6468094 38808564 126931853 ALL.chr1.bim 2504 15024 62600 ALL.chr1.fam 28 134 1074 ALL.chr1.log 6468094 25872376 100282585 ALL.chr1.map 2504 5008 40064 ALL.chr1.nosex
bimファイルの6468094がSNP数を示す。 famファイルの2504がサンプル数を示す。
4. PCA解析を実行する。
plinkの機能を使って、主成分分析(PCA)を行う。
plink2 --bfile ALL.chr1 --pca approx 10 --out qcvcf
解析が完了すると、qcvcf.eigenval, qcvcf.eigenvec, qcvcf.eigenvec.txtが生成される。
お疲れ様でした、今回はこれで終わりです。
次回はRを使ってプロットを作成します。
よければ他の記事も見ていってください。