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

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

MENU

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

plinkの使い方

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

使用するRのパッケージ: tidyverse, dplyr, cowplot

使用するRのコマンド: read.delim, read.csv, colnames, merge, ggplot, plot_grid


この記事を読むと何ができるようになる?
plinkを使ったGWAS解析の演習を行います。この記事は前回の続きです。


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

前半では、国際的なヒトゲノムプロジェクト(1000 Genome Project)から遺伝子型データと出身地データを取得し、plinkを使って集団間の遺伝的距離と地理的な距離を主成分分析(PCA)を行いました。今回は、PCA結果をプロットした図をRを使って作成し、地域ごとにクラスターを形成していることを体験してみましょう。

eupatho-bioinfomatics.hatenablog.com

1. 前回取得したファイルを確認する。

今回の作業に必要なファイルは、以下の2つである。

qcvcf.eigenvec   #plinkでPCAを実行して生成されたファイル。
integrated_call_samples_v3.20200731.ALL.ped #サンプル名と出身地のデータ

ファイルの中身を最初の10行を表示するheadを使って表示すると、下のようにサンプル名とPC1からPC10の各成分値が表示される。

$ head qcvcf.eigenvec
#FID    IID PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
HG00096 HG00096 -0.0118373  -0.0242125  0.0151669   0.0174092   0.000316254 0.00041284  -0.000594013    0.0014243   0.00106578  -0.000591639
HG00097 HG00097 -0.0114178  -0.0267128  0.0147753   0.0187198   0.000772494 0.000983227 0.00026169  -0.00054953 3.9298e-05  -0.000449919
HG00099 HG00099 -0.0118923  -0.0259647  0.0156335   0.0183097   0.000429924 -0.000300371    -0.00174736 0.000186716 -0.000495259    -3.5582e-05
HG00100 HG00100 -0.0118524  -0.0259545  0.0157248   0.0154537   0.000780088 0.00082175  0.0012578   -0.000951344    0.00107604  0.000841985
HG00101 HG00101 -0.0120133  -0.0267428  0.0163999   0.0174753   0.00150347  0.000780043 -0.00135751 0.00105417  0.000284864 -0.000215038
HG00102 HG00102 -0.0111609  -0.0252704  0.0154327   0.0182801   0.000520204 0.000534406 -0.000548787    0.00041653  0.000135774 0.000475114
HG00103 HG00103 -0.0119031  -0.0256639  0.0147533   0.0193582   0.000486153 -0.00107121 -0.000550418    -0.000246204    0.000778476 -0.00145126
HG00105 HG00105 -0.0121201  -0.0257963  0.0155803   0.0213999   0.000925205 0.0011147   -0.00089272 0.000238308 0.000148134 0.000592621
HG00106 HG00106 -0.0125321  -0.0242146  0.0150248   0.0166768   0.00072697  -0.000118516    0.000203522 0.000848591 0.00106597  -0.000408958

続いて出身地データを確認する。
1行目にFamily ID, 2行目にIndividual ID, 7行目に国名が記されている。

$ head integrated_call_samples_v3.20200731.ALL.ped
Family ID,Individual ID,Paternal ID,Maternal ID,Gender,Phenotype,Population,Relationship,Siblings,Second,Order,Third,Order,Children,Other,Comments,phase,,,,,3,genotypes,related,genotypes,omni,genotypes,affy_genotypes,,
HG00096,HG00096,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00097,HG00097,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00098,HG00098,0,0,1,0,GBR,unrel,0,0,0,0,0,0,0,1,1,,,,,,,,,,,,,,,,,,,
HG00099,HG00099,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00100,HG00100,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00101,HG00101,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00102,HG00102,0,0,2,0,GBR,unrel,0,0,0,0,0,1,0,1,1,,,,,,,,,,,,,,,,,,,
HG00103,HG00103,0,0,1,0,GBR,unrel,0,0,0,0,0,1,0,1,0,,,,,,,,,,,,,,,,,,,
HG00104,HG00104,0,0,2,0,GBR,unrel,0,0,0,0,0,0,0,1,1,,,,,,,,,,,,,,,,,,,

2. PCA結果をプロットする。

RStudioを起動して、以下のコマンドを順に実行する。

1 掃除コマンド

rm(list=ls())
setwd("~/path/to/worlingdir")
#必要なパッケージをロードする。
library("tidyverse")
library("dplyr")

eupatho-bioinfomatics.hatenablog.com

2結果ファイルを読み込む。

pca <- read.delim("./qcvcf.eigenvec", header=TRUE)

3 FamilyIDが不要なので列を削除する。

pca <- pca[,-1]
head(pca)

4 出身地を記載した"integrated_call_samples_v3.20200731.ALL.csv"を読み込む。

data <- read.csv('integrated_call_samples_v3.20200731.ALL.csv', header=TRUE)
head(data)

5 必要な列を取り出す(https://qiita.com/hitsujisuke/items/d71ee40daa0786ae1680)。

columnList <- c("Family.ID","Population") #カラムリストの作成する。
data.selected <- data[, columnList] #データの抽出する。
head(data.selected)

6 結合するためにカラム名を揃える。

colnames(data.selected) <- c("IID","Population")
head(data.selected)

7 共通の列に基づいて、データテーブルを結合する。

D <- merge(pca, data.selected)
head(D)

8 ggplotで国別に色分けして描画する。

PC1,PC2
b1 <- ggplot(D, aes(PC1, PC2, col = Population)) + geom_point(alpha = 0.5, size = 2)
b1 <- b1 + coord_equal() + theme_light()
b1

PC3,PC4
b2 <- ggplot(D, aes(PC3, PC4, col = Population)) + geom_point(alpha = 0.5, size = 2)
b2 <- b2 + coord_equal() + theme_light()
b2

9 プロットを並べてpdfファイルに保存する。

library("cowplot")
pdf(file = "PCA.pdf")
plot_grid(b1, b2, align = "h")
dev.off()

同じフォルダ内に、PCA.pdfが生成される。

f:id:Harry-kun:20210627082330j:plain
ゲノムデータに対する主成分分析(PCA)で得られた成分値をプロットすると、人種によって集団構造を持つことがわかる。

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

お疲れ様でした、今回はこれで終わりです。

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