バイオインフォマティクスでゲノムワイド関連解析(GWAS)

バイオインフォマティクスを頑張っている方が、本ブログの内容を真似することで、自分のデータで解析ができる情報を提供することが目標です! 今はGATKの解説をメインテーマにしています。

MENU

plinkを使ってゲノムの集団構造を体験する-パート1

plinkの使い方

導入難易度★☆☆☆☆
使用難易度★★★☆☆

この記事を読むと何ができるようになる?
今回から2回に分けて、plinkを使ったGWAS解析の演習を行います。まず、国際的なヒトゲノムプロジェクト(1000 Genome Project)により得られたシークエンスデータを公表しているサイトから、遺伝子型データと出身地データを取得します。次に、plinkを使って、集団間の遺伝的距離と地理的な距離を主成分分析(PCA)を行って解析します。最後に、結果をプロットした図をRを使って作成し、地域ごとにクラスターを形成していることを体験してみましょう。

今回は、解析に必要なデータの入手と主成分分析までを行います。

なお、本記事で扱う実習内容は「ゼロから実践する遺伝統計学セミナー」(羊土社, p128-135)の演習内容を参考にさせていただきました。この本は初めてGWAS解析を実施する人が知りたい内容が丁寧かつコンパクトに説明されており、他にも色々なplinkの使用方法、QTL解析やeQTL解析からマンハッタンプロットの描画まで、多くの演習内容が詰まっています。GWAS解析をこれから始めてみたい方には、非常におすすめです。


Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。

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

f:id:Harry-kun:20210626105148p:plain
1000 Genome Projectのトップページ

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を使ってプロットを作成します。

よければ他の記事も見ていってください。

eupatho-bioinfomatics.hatenablog.com