PLINKで主成分分析＆Rでプロット

主成分分析（PCA: principal components analysis）はSNPデータに基づく集団構造の推定や系統推定のための予備解析によく使われる手法です。MIG-seqやRAD-seq、GRAS-Diなどのゲノム縮約法だけでなく、全ゲノム解析でもPCAは頻繁に使われます。そのため、解析から結果の可視化まで簡単に実装できると嬉しいはず。ここでは、機能が多すぎるゲノム解析プログラムであるPLINKを用いてPCAを実行し、Rを使って結果をプロットするシェルスクリプトを紹介します。

環境

macOS 11.6.1
PLINK v1.90b6.24 64-bit (6 Jun 2021)  
R version 4.1.1 (2021-08-10)

インストール

公式HPからダウンロードしてください。パスを通しておくと便利です。

インプットファイル

ここではMIG-seqやRAD-seqのデータを想定して、Stacksのpopulationsから書き出したSNPデータを使います。
populationsの書き出しで--plinkフラッグを付けるとpopulations.plink.mapとpopulations.plink.pedが出力されます。この2つのファイルがインプットとなります。

Stacksの使い方は以下を参照してください。

Stacks以外で出力したinputを使う場合は、スクリプトのinput=の部分でinputファイル名を適宜修正してください。
集団名の指定方法が違う可能性があり、スクリプトがうまく動作しないかもしれません。気合で頑張ってください。

解析

以下のシェルスクリプトをinputファイルと同じディレクトリに置き、sh plink-pca-plot.shで実行します。
PCAの条件や図の出力サイズを変更したい場合はシェルスクリプトの# 変数の定義の所を書き換えます。プロットの色や形を変えたい場合はRスクリプトの中身を直接編集してください。

このスクリプトはPLINKを実行した後に、プロットのためのRscriptを生成し、ターミナル上で実行します。

plink-pca-plot.sh

#!/bin/zsh

# Author: YF (2022-01-07)
# Stacksのpopulationsで書き出したSNPデータをplinkを使ってPCA解析する
# Rscriptを生成・実行して、PCAの結果をプロットする

# 変数の定義。PCAやプロットの条件を変える時はここを変更する
input="populations.plink"  #入力ファイル名を指定
out="pca"      #PCAの出力ファイル名を指定
npc=20          #PCの数を指定
format="pdf"    #プロットの出力フォーマット e.g., pdf, jpg, png
prefix="out"    #プロットのファイル名
width=20        #プロットの幅（cm）
height=15       #プロットの高さ（cm）

# PCA
plink --file ${input} --out ${out} --pca header ${npc} --allow-extra-chr
# --allow-extra-chr フラッグはStacksで出力した時以外は不要かもしれない

# generate Rscript
cat << EOS > pca-plot.R
rm(list=ls())
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(viridis))

pca <- read.table("${out}.eigenvec", header=T)
eigenval <- read.csv("${out}.eigenval", header=F)
df <- data.frame(pc = 1:nrow(eigenval), eigenval/sum(eigenval)*100)

p1 <- ggplot(df, aes(pc, V1)) +
  geom_bar(stat = "identity")+
  xlab("Percentage components")+
  ylab("Variance explained (%)")+
  theme_bw(base_size = 15)

ggsave("${prefix}_barplot.${format}", p1,
         width = ${width}, height = ${height},
         units = "cm", dpi = 300)

p2 <- ggplot(pca, aes(PC1, PC2, color=FID, shape=FID))+
  geom_point(size=5, alpha=0.9)+
  scale_shape_manual(values=c(rep(15:18, 20)))+
  scale_color_viridis(option = "D", discrete = T,end = 0.99)+
  xlab(paste0("PC1 (", round(df[1,2], digits = 2), "%)"))+
  ylab(paste0("PC2 (", round(df[2,2], digits = 2), "%)"))+
  theme_bw(base_size = 15)

ggsave("${prefix}_PC1PC2plot.${format}", p2, 
         width = ${width}, height = ${height},
         units = "cm", dpi = 300)

p3 <- ggplot(pca, aes(PC3, PC4, color=FID, shape=FID))+
  geom_point(size=5, alpha=0.9)+
  scale_shape_manual(values=c(rep(15:18, 20)))+
  scale_color_viridis(option = "D", discrete = T,end = 0.99)+
  xlab(paste0("PC3 (", round(df[3,2], digits = 2), "%)"))+
  ylab(paste0("PC4 (", round(df[4,2], digits = 2), "%)"))+
  theme_bw(base_size = 15)

ggsave("${prefix}_PC3PC4plot.${format}", p3, 
         width = ${width}, height = ${height},
         units = "cm", dpi = 300)

p4 <- ggplot(pca, aes(PC2, PC3, color=FID, shape=FID))+
  geom_point(size=5, alpha=0.9)+
  scale_shape_manual(values=c(rep(15:18, 20)))+
  scale_color_viridis(option = "D", discrete = T,end = 0.99)+
  xlab(paste0("PC2 (", round(df[2,2], digits = 2), "%)"))+
  ylab(paste0("PC3 (", round(df[3,2], digits = 2), "%)"))+
  theme_bw(base_size = 15)

ggsave("${prefix}_PC2PC3plot.${format}", p4, 
         width = ${width}, height = ${height},
         units = "cm", dpi = 300)

EOS

# plot
Rscript pca-plot.R

PLINKで3つの出力、Rで4つの出力があります。

1. pca.log　PLINKのログ
2. pca.eigenvec　PCAの結果
3. pca.eigenval　PCスコアの一覧
4. out_barplot.pdf　主成分の寄与率の棒グラフ
5. out_PC1PC2plot.pdf　PC1とPC2の散布図
6. out_PC2PC3plot.pdf　PC2とPC3の散布図
7. out_PC3PC4plot.pdf　PC3とPC4の散布図

出力される主成分の寄与率の棒グラフの例

出力される散布図の例　ラベルに寄与率が入ります。

---

不具合があればお知らせください。
ただし、Stacks以外での出力ファイルを使った場合の不具合とWindowsやLinux環境での不具合には対応できない可能性が高いです。

参考サイト