plinkを使ってゲノムの集団構造を体験する-パート2
plinkの使い方
導入難易度★☆☆☆☆
使用難易度★★★☆☆
使用するRのパッケージ: tidyverse, dplyr, cowplot
使用するRのコマンド: read.delim, read.csv, colnames, merge, ggplot, plot_grid
この記事を読むと何ができるようになる?
plinkを使ったGWAS解析の演習を行います。この記事は前回の続きです。
Twitterで記事の更新をお知らせしているので、興味を持たれた方は是非フォローをお願いします。
フォローする @harrykun_blog
前半では、国際的なヒトゲノムプロジェクト(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が生成される。
なお、本記事で扱う実習内容は「ゼロから実践する遺伝統計学セミナー」(羊土社, p128-135)の演習内容を参考にさせていただきました。この本は初めてGWAS解析を実施する人が知りたい内容が丁寧かつコンパクトに説明されており、他にも色々なplinkの使用方法、QTL解析やeQTL解析からマンハッタンプロットの描画まで、多くの演習内容が詰まっています。GWAS解析をこれから始めてみたい方には、非常におすすめです。本記事の一番下にAmazonのリンクを貼っておきます。
お疲れ様でした、今回はこれで終わりです。
よければ他の記事も見ていってください。