PopCluster (Wang 2022) のAdmixture analysisの結果を自動で図示する

2023-02-13

新しい集団構造解析ソフトウェアPopCluster

Wang, J. Fast and accurate population admixture inference from genotype data from a few microsatellites to millions of SNPs. Heredity 129, 79–92 (2022). https://doi.org/10.1038/s41437-022-00535-z

TさんのPopCluster紹介転載

（Wangさん）本人が言うには，kinshipも考慮してくれて，STRUCTUREやADMIXTUREではでないもんがでる．と言うてました．STRUCTUREは時間かかるしデータが多いとアカン．ADMIXTUREはサンプルがアンバランスやと全然アカン．

それらをクリアしたのが「PopCluster」である．と言うてました．

PopClusterの難点のひとつは，結果の図示が面倒くさいところです。PopClusterの出力ファイルは人間が見やすいといえば見やすいのですが，図示に適したフォーマットにはなっていません。

そこで，PopClusterの出力結果とpopmap.txtさえあれば，図示まで自動でやってくれるスクリプトを書きました。

準備

PopClusterのAdmixture analysisを回しておく。各Kについて何回か繰り返しをするとよい（たとえば100回）。VCFからインプットファイルを作りたいなら， #VCFからPopCluster_(Wang_2022)_用data_fileを生成するスクリプトも参照。

gsedのインストール（Macならbrewで入る；Linuxならシェルスクリプト内のgsedをsedに変えれば動くはず）

Rのインストール（パッケージ tidyverseが必要）

popmap.txtの準備（サンプル名<tab>集団名）

以下の1と2のスクリプト（popcluster_plot.sh，popcluster_plot.R）をPopClusterの出力結果があるディレクトリに置く

実行

ターミナルでPopClusterの出力結果があるディレクトリに移動し

sh popcluster_plot.sh <INPUT_PREFIX> <OUT_DIR> <NUM_INDS> <Kmin> <Kmax> <path/to/popmap.txt>

を実行。<hoge>の部分は自分の解析に合わせて書き換える。

例：sh popcluster_plot.sh all ./plot 131 1 17 ./popmap.txt

結果はこんな感じ。

https://scrapbox.io/files/63ea4d16394d71001b327431.png

※ 不具合があったらご報告ください！

1｜PopClusterの出力から図示に必要な部分を抜粋・整形して保存し，2のRスクリプトをランするシェルスクリプト

以下の引数をとる。全ての引数が必要。

第1引数：PopClusterのパラメータファイル（〜.PcPjt）で指定したOutputfilename（INPUT_PREFIX）。

第2引数：図示結果を出力するディレクトリ（OUT_DIR）。

第3引数：解析した個体数（NUM_INDS）。

第4引数：図示したいKの最小値（Kmin）。

第5引数：図示したいKの最大値（Kmax）。99以下を想定している。

第6引数：popmap.txtへのパス（path/to/popmap.txt）。

出力されるファイルのフォーマット（.pngとか）や縦横比を変えたい場合は，このスクリプトの最終行をいじるとよい。

うまく走れば数秒で結果が出る。

1分以上かかる場合は，なんらかの異常が起きている。タスクキルして，引数にミスがないか確認すること。

code:popcluster_plot.sh

# usr/bin/sh

# Keisuke Onuki 2023-02-08

# usage: sh popcluster_plot.sh <INPUT_PREFIX> <OUT_DIR> <NUM_INDS> <Kmin> <Kmax> <path/to/popmap.txt>

# ex: sh popcluster_plot.sh all.131.miss0.9.ld0.5 ./plot 131 1 17 ./popmap.txt

# Kmaxは99以下を想定している

#################

INPUT_PREFIX=$1

OUT_DIR=$2

NUM_INDS=$3

Kmin=$4

Kmax=$5 #99以下を想定

popmap=$6

#################

mkdir $OUT_DIR

for i in seq 2 $((Kmax-Kmin+2))

# 並べ替えのためのID作成（ex: 5->05, 12->12）

ID=$((i-2+Kmin))

if $ID -lt 10 ; then

ID=0$ID

else

ID=$ID

# Best runファイル名を取得してINPUTに格納

INPUT=awk 'NR=='$i' {print $2}' ${INPUT_PREFIX}.K

# K=iのPopClusterの結果のうち，Admixture analysisの結果部分が始まる行を取得してADMIXTURE_STARTに格納

ADMIXTURE_START=grep -n "Inferred ancestry of individuals" ${INPUT} | awk -F ':' '{print $1}'

# 取得した行の次の行からサンプル数分の行を取得・整形（複数あるスペースの削除，行頭のスペースの削除，「: 」の削除）

# OUT_DIRに.r_inputファイルとして保存

awk 'NR=='$ADMIXTURE_START'+2,NR=='$ADMIXTURE_START'+'$NUM_INDS'+2 {print}' ${INPUT}| \

gsed 's/ */ /g' | \

gsed 's/^ //g'| \

gsed 's/: //g' > ${OUT_DIR}/${INPUT_PREFIX}_K_${ID}.r_input

done

# プロット

# 上だと出力のサイズやフォーマットを指定できる．

# 下だと出力サイズは自動計算してくれる．出力はPDF固定．

# Rscript popcluster_plot.R $OUT_DIR $popmap $Kmin $Kmax $OUT_DIR 20 7.5 .pdf

Rscript popcluster_plot.R $OUT_DIR $popmap $Kmin $Kmax $OUT_DIR

うまくいけば，第2引数（OUT_DIR）で指定したディレクトリ内に以下のファイルが出力される。

plot_K<Kmin>–<Kmax>.pdf （指定したKminからKmaxまでのプロット）

<INPUT_PREFIX>_K_X.r_input （図示に用いたシングルスペース区切りテキストファイル；図示するKの数だけ出力）

2｜整形された結果をプロットするRスクリプト

1のスクリプト内でランされるRスクリプト。このスクリプトではKが20以上になると色分けできないので注意。

original_20に任意の色セットを代入すれば，99色まではOK。

code:popcluster_plot.R

# PopClusterのRプロット

# Keisuke Onuki 2023-02-08

# dependency: tidyverse

# Usage: call script from shell with 5 (+ optional 3) arguments.

# Rscript popcluster_plot.R <path to edited PopCluster output files (.r_input)> <popmap> <low K to plot> <high K to plot> <output dir> <fig width> <fig height> <fig format>

# e.g.

# Rscript popcluster_plot.R ../output/tools/PopCluster ../data/popmap/popmap.txt 2 5 ../output/tools/PopCluster

# Rscript popcluster_plot.R ../output/tools/PopCluster ../data/popmap/popmap.txt 2 5 ../output/tools/PopCluster 50 4 .pdf

# rm(list = ls()); gc(); gc()

suppressPackageStartupMessages(library(tidyverse))

# get arguments from shell

args <- commandArgs(TRUE)

# default arguments

popcluster_out_dir <- "./" # path to edited PopCluster output directory which contain all .r_input files

popmap <- "./popmap.txt" # popmap file

k_range <- c(4:5) # c(2, 4, 6) is acceptable. # K values which want to plot

out_dir <- "./"

# substitute by arguments

popcluster_out_dir <- args1

popmap <- args2

k_range <- c(args3:args4)

out_dir <- args5

pop <- read.table(popmap, header = FALSE, col.names = c("ID", "pop"))

lf <- list.files(path = popcluster_out_dir, pattern = ".r_input\\>", ignore.case = TRUE, full.names = TRUE) # ディレクトリ内の.r_inputファイル名をリストで取得

data <- map(lf, read.table, header = FALSE) # .r_input ファイルをリストにまとめる。lapplyでも可。

kmax <- length(lf)

for (i in 1:kmax) {

colnames(datai) <- c("Index","Order","Label","%Miss","Cluster",paste0("V",sprintf("%02d", 1:i)))

datai <- datai %>%

mutate(K = i) %>%

cbind(pop)

}

data_long <- list()

data_long_all <- NULL

for (i in 1:kmax) {

data_longi <- datai %>%

pivot_longer(cols = starts_with("V"),

names_to = "clus",

names_prefix = "v",

values_to = "coef")

data_long_all <- rbind(data_long_all, data_longi)

}

original_20 <- c("#3366cc", "#dc3912", "#ff9900","#109618","#990099","#0099c6","#dd4477","#66aa00","#b82e2e","#316395","#994499","#22AA99","#AAAA11","#6633CC","#E67300","#8B0707","#651067","#329262","#5574A6","#3B3EAC")

p <- data_long_all %>%

filter(K %in% k_range) %>%

mutate(kname = str_c("K = ", K) %>% fct_inorder()) %>%

ggplot(aes(x = factor(ID), y = coef, fill = factor(clus))) +

geom_col(color = "gray", size = 0.1) + # ggplot2 v3.4.0以降では，sizeをlinewidthに変えた方がよい

facet_grid(kname ~ fct_inorder(pop), switch = "x", scales = "free", space = "free") +

theme_minimal() +

labs(x = "Individuals", y = "Ancestry") +

scale_y_continuous(expand = c(0, 0)) +

scale_x_discrete(expand = expansion(add = 1)) +

theme(panel.spacing.x = unit(0.1, "lines"),

# axis.text.x = element_blank(),

axis.text = element_blank(),

panel.grid = element_blank()) +

scale_fill_manual(values=original_20, guide = "none")

# default output format is pdf

# default fig size is derived from number of individuals and K

# users can adjust these values

width = ifelse(is.na(args6), nrow(pop) / 4, as.numeric(args6))

height = ifelse(is.na(args7), length(k_range) * 6, as.numeric(args7))

format = ifelse(is.na(args8), ".pdf", args8)

str_c(out_dir, "/plot_K", min(k_range), "-", max(k_range), format) %>%

ggsave(p,

width = width,

height = height,

units = "cm",

dpi = 300)

writeLines(str_c("Figure: written in ", out_dir, " plot_K", min(k_range), "_", max(k_range), format))