GWAS全流程#02 | 数据准备

写在前面

gwas整套流程按顺序如下：

• 1.测序数据质控（fastp）
• 2.比对&标记重复（bwa，samtools）
• 3.变异检测（gatk）
• 4.关联分析

前面一篇文章已经讲了以上的前3个内容，这次说最后关联分析的内容。

4.关联分析

在进行关联分析之前，需要对我们前面获得的变异位点信息文件.vcf进行数据质控过滤，并对测序个体进行群体结构和连锁不平衡分析(部分结果作为协变量输入)。

1.数据质控过滤

数据质控分为个体质控和位点质控，即有些测序个体是否存在大量位点缺失的情况，有些变异位点缺少部分个体的变异信息。这些个体和变异位点可能需要去除。

• ①个体质控（一般不做）
  事实上我们一般不进行个体过滤，只有在进行位点过滤后变异位点较少的情况下，才会考虑：是否因为某些个体因为有效测序深度过低导致大量位点缺失，进而在位点过滤时去掉这些位点而导致过少的变异位点。这些个体可以舍去。

#1.生成缺失率文件plink.imiss (按样本统计缺失率), plink.lmiss (按位点统计缺失率)
plink --vcf all.filtered.snp.vcf.gz \
     --allow-extra-chr \
     --missing \
     --out snp
#2.个体缺失率查看，这里可以看到变异位点缺失最多的样本也才5.7%的缺失率
#如果有20%的缺失率，那这个样本测得可能有问题，可以舍去，后面遇到再说
(base) [wangtao01@sonmi:~/gwas/tana/gatk_joint/merged_vcf]$ sort -k6,6gr snp.imiss | head
 SRR17015761   SRR17015761          Y   634107 11094943  0.05715
 SRR17015762   SRR17015762          Y   565257 11094943  0.05095
 SRR17015732   SRR17015732          Y   542475 11094943  0.04889
 SRR17015737   SRR17015737          Y   433534 11094943  0.03907
 SRR17015741   SRR17015741          Y   400748 11094943  0.03612
 SRR17015765   SRR17015765          Y   384860 11094943  0.03469
 SRR17015775   SRR17015775          Y   329450 11094943  0.02969
 SRR17015773   SRR17015773          Y   324598 11094943  0.02926
 SRR17015759   SRR17015759          Y   324127 11094943  0.02921
 SRR17015760   SRR17015760          Y   320872 11094943  0.02892

• ②位点质控
  一般用进行缺失率（missing genotypes）和最小等位基因频率（minor allele frequency，MAF）进行质控。后续进行gwas分析用plink质控，bcftools工具可用非gwas分析的数据质控比较方便。

plink --vcf merge.filtered.gatk.vcf.gz \
      --allow-extra-chr \    #允许不规范的染色体名称
      --set-missing-var-ids @:#   \   #生成变异位点ID，需加
      --biallelic-only strict \    #只保留二等位位点
      --geno 0.1 \  #去掉缺失率 >10% 的位点
      --mind 0.1 \   #去掉缺失率 >10% 的样本，一般不设置
      --maf 0.05 \
      --make-bed \   #输出bed格式文件，作为gwas主流输入文件
      --recode vcf-iid  \   #同时输出vcf文件备份
      --out gwas.qc

#过滤后的文件
-rw-r--r-- 1 wangtao01 wangtao01   22268443 Dec 18 16:26 gwas.qc.bed
-rw-r--r-- 1 wangtao01 wangtao01   81443920 Dec 18 16:26 gwas.qc.bim
-rw-r--r-- 1 wangtao01 wangtao01       1320 Dec 18 16:26 gwas.qc.fam
-rw-r--r-- 1 wangtao01 wangtao01       1328 Dec 18 16:26 gwas.qc.log
-rw-r--r-- 1 wangtao01 wangtao01        960 Dec 18 16:26 gwas.qc.nosex
-rw-r--r-- 1 wangtao01 wangtao01  455555358 Dec 18 16:26 gwas.qc.vcf

#输出结果日志：80%被过滤掉了
# 10841114 variants loaded
# 425867 variants removed due to missing genotype data (--geno)
# 8188403 variants removed due to minor allele threshold (--maf)
# 2226844 variants pass filters
$cat gwas.qc.log
PLINK v1.9.0-b.7.7 64-bit (22 Oct 2024)
Options in effect:
  --allow-extra-chr
  --biallelic-only strict
  --geno 0.1
  --maf 0.05
  --make-bed
  --out gwas.qc
  --recode vcf-iid
  --set-missing-var-ids @:#
  --vcf all.filtered.snp.vcf.gz

Hostname: sonmi
Working directory: /share/home/wangtao01/gwas/tana/gatk_joint/merged_vcf
Start time: Thu Dec 18 16:25:05 2025

Random number seed: 1766046305
515089 MB RAM detected; reserving 257544 MB for main workspace.
--vcf: gwas.qc-temporary.bed + gwas.qc-temporary.bim + gwas.qc-temporary.fam
written.
(253829 variants skipped.)
10841114 variants loaded from .bim file.
10841114 missing IDs set.
40 people (0 males, 0 females, 40 ambiguous) loaded from .fam.
Ambiguous sex IDs written to gwas.qc.nosex .
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 40 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.975406.
425867 variants removed due to missing genotype data (--geno).
8188403 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
2226844 variants and 40 people pass filters and QC.
Note: No phenotypes present.
--make-bed to gwas.qc.bed + gwas.qc.bim + gwas.qc.fam ... done.
--recode vcf-iid to gwas.qc.vcf ... done.

• ③准备表型数据
  我这里做的性别关联分析，样本的表型就是性别信息，是雌性还是雄性，这里把雄性数据编码为1，雌性数据编码为2。我这里准备一个包含3列（FID IID SEX）的表型数据文件sex.pheno.txt：

cat sex.pheno.txt
SRR17015732     SRR17015732     1
SRR17015733     SRR17015733     1
SRR17015734     SRR17015734     2
SRR17015735     SRR17015735     2
SRR17015736     SRR17015736     1
SRR17015737     SRR17015737     2
SRR17015738     SRR17015738     2
SRR17015739     SRR17015739     2
SRR17015740     SRR17015740     1
SRR17015741     SRR17015741     1
SRR17015742     SRR17015742     2
SRR17015743     SRR17015743     1
SRR17015744     SRR17015744     1
SRR17015745     SRR17015745     1
SRR17015746     SRR17015746     2
SRR17015747     SRR17015747     2
SRR17015748     SRR17015748     2
SRR17015749     SRR17015749     1
SRR17015750     SRR17015750     1
SRR17015751     SRR17015751     1
SRR17015752     SRR17015752     1
SRR17015753     SRR17015753     1
SRR17015754     SRR17015754     2
SRR17015755     SRR17015755     2
SRR17015756     SRR17015756     2
SRR17015757     SRR17015757     2
SRR17015758     SRR17015758     2
SRR17015759     SRR17015759     1
SRR17015760     SRR17015760     1
SRR17015761     SRR17015761     2
SRR17015762     SRR17015762     1
SRR17015763     SRR17015763     2
SRR17015765     SRR17015765     1
SRR17015767     SRR17015767     2
SRR17015769     SRR17015769     2
SRR17015771     SRR17015771     1
SRR17015773     SRR17015773     1
SRR17015775     SRR17015775     2
SRR17015777     SRR17015777     1
SRR17015779     SRR17015779     2

2.群体结构和连锁不平衡分析

全基因组关联分析gwas本质原理是，利用标记（变异位点）与目标位点之间存在的连琐不平衡，定位到与目标性状关联的连琐不平衡区域。当然，我们用于定位性状的各个位点显著性p值的是基于线性回归模型计算的，群体结构和个体间的亲缘关系会影响该值，因此使用关联分析的模型需要输入群体结构和亲缘关系数据对模型进行一定的矫正。

• ①输入数据调整
  尽管我么已经准备质控过的变异文件,vcf，1.但是由于我们的染色体名称不规范（规范染色体名称是1，2，3，4，5......，即纯数字），会导致下游一些分析关键报错，因此，需要将染色体名称进行调整；其次，2.下游structure分析需要排除连琐不平衡位点，因此也需要将这部分变异位点进行去除。

#1.调整染色体名称
#染色体名称在3个BED文件中的bim中修改
awk 'NR==FNR{m[$1]=$2; next}
     {if($1 in m){$1=m[$1]; $2=$1":"$4} print}' \
    chr_map.txt gwas.qc.bim > gwas.qc.chrnum.bim

ln -s  gwas.qc.bed  gwas.qc.chrnum.bed
ln -s  gwas.qc.fam  gwas.qc.chrnum.fam

#查看替换前后结果
(base) [wangtao01@sonmi:~/gwas/tana/gwas_workdir]$ head gwas.qc.bim gwas.qc.chrnum.bim
==> gwas.qc.bim <==
chr1A   chr1A:52447     0       52447   G       A
chr1A   chr1A:52891     0       52891   A       G
chr1A   chr1A:53186     0       53186   G       T
chr1A   chr1A:53391     0       53391   C       T
chr1A   chr1A:53434     0       53434   G       A
chr1A   chr1A:53517     0       53517   G       A
chr1A   chr1A:53930     0       53930   A       G
chr1A   chr1A:53950     0       53950   C       T
chr1A   chr1A:54060     0       54060   T       G
chr1A   chr1A:54198     0       54198   A       C

==> gwas.qc.chrnum.bim <==
1 1:52447 0 52447 G A
1 1:52891 0 52891 A G
1 1:53186 0 53186 G T
1 1:53391 0 53391 C T
1 1:53434 0 53434 G A
1 1:53517 0 53517 G A
1 1:53930 0 53930 A G
1 1:53950 0 53950 C T
1 1:54060 0 54060 T G
1 1:54198 0 54198 A C

#2.去除连琐不平衡LD位点
##第一步：计算 LD 剪枝位点列表（生成 prune.in / prune.out）
#注意！！--indep-pairwise 接受三个参数，
#第一个是窗口大小（可以是 位点数variant count 或者带 'kb' 指定为 kilobase），默认没有 'kb' 时，window size 是 SNP 个数 而不是物理距离
#第二个是 step size（位点数量）
#第三个是阈值（ r²）
plink --bfile gwas.qc.chrnum \
      --indep-pairwise 50 10 0.2 \
      --out gwas.qc.chrnum.ld
##第二步：按 prune.in 提取，生成 LD-pruned 的新数据集
plink --bfile gwas.qc.chrnum \
      --extract gwas.qc.chrnum.ld.prune.in \
      --make-bed \
      --out gwas.qc.chrnum.LD

#得到可作为群体结构structure分析的编译位点文件gwas.qc.chrnum.LD
(base) [wangtao01@sonmi:~/gwas/tana/gwas_workdir]$ ll -tr
total 643036
-rw-r--r-- 1 wangtao01 wangtao01      1041 Dec 18 17:07 sex.pheno.txt
-rw-r--r-- 1 wangtao01 wangtao01  81443920 Dec 18 17:09 gwas.qc.bim
-rw-r--r-- 1 wangtao01 wangtao01  22268443 Dec 18 17:09 gwas.qc.bed
-rw-r--r-- 1 wangtao01 wangtao01 455555358 Dec 18 17:09 gwas.qc.vcf
-rw-r--r-- 1 wangtao01 wangtao01       960 Dec 18 17:09 gwas.qc.nosex
-rw-r--r-- 1 wangtao01 wangtao01      1328 Dec 18 17:09 gwas.qc.log
-rw-r--r-- 1 wangtao01 wangtao01      1320 Dec 18 17:09 gwas.qc.fam
-rw-r--r-- 1 wangtao01 wangtao01       222 Dec 19 17:46 chr_map.txt
-rw-r--r-- 1 wangtao01 wangtao01  63629168 Dec 19 18:17 gwas.qc.chrnum.bim
lrwxrwxrwx 1 wangtao01 wangtao01        11 Dec 19 18:22 gwas.qc.chrnum.bed -> gwas.qc.bed
lrwxrwxrwx 1 wangtao01 wangtao01        11 Dec 19 18:22 gwas.qc.chrnum.fam -> gwas.qc.fam
-rw-r--r-- 1 wangtao01 wangtao01       960 Dec 19 18:41 gwas.qc.chrnum.ld.nosex
-rw-r--r-- 1 wangtao01 wangtao01   3033528 Dec 19 18:41 gwas.qc.chrnum.ld.prune.in
-rw-r--r-- 1 wangtao01 wangtao01  22100524 Dec 19 18:41 gwas.qc.chrnum.ld.prune.out
-rw-r--r-- 1 wangtao01 wangtao01      2426 Dec 19 18:41 gwas.qc.chrnum.ld.log
-rw-r--r-- 1 wangtao01 wangtao01       960 Dec 19 18:48 gwas.qc.chrnum.LD.nosex
-rw-r--r-- 1 wangtao01 wangtao01      1320 Dec 19 18:48 gwas.qc.chrnum.LD.fam
-rw-r--r-- 1 wangtao01 wangtao01   2698783 Dec 19 18:48 gwas.qc.chrnum.LD.bed
-rw-r--r-- 1 wangtao01 wangtao01   7686324 Dec 19 18:48 gwas.qc.chrnum.LD.bim
-rw-r--r-- 1 wangtao01 wangtao01       980 Dec 19 18:48 gwas.qc.chrnum.LD.log

• ②亲缘关系kinship计算
  使用GCTA软件基于IBS的方式计算亲缘关系，并行格式转换成常规的对角矩阵文件

##1.使用GCTA软件计算亲缘关系
#注意加上-autosome-num 参数，GCTA默认按照人的染色体名称，超过或者不匹配染色体名称则会报错
gcta --bfile gwas.qc.chrnum \
     --make-grm-gz \
     --make-grm-alg 0 \
     --autosome-num 24 \
     --out gcta_kinship
#输出如下两个文件，.grm.gz即grm阵列（genetic relationship matrix）文件，需要转成对角矩阵用于下游分析
-rw-r--r-- 1 wangtao01 wangtao01      9512 Dec 21 00:58 gcta_kinship.grm.gz
-rw-r--r-- 1 wangtao01 wangtao01       960 Dec 21 00:58 gcta_kinship.grm.id

##2.替换样品名称，转成对角矩阵
#用如下脚本grm_to_matrix.py进行替换，注意更改输入文件前缀：
——————————————————————————————————————————————————————————————-————————————————————
python3 - <<'PY'
import gzip
import numpy as np

prefix = "gcta_kinship"   # ⭐ 如果你的前缀不一样，就改这里

# 1. 读取样本ID（行列名）
ids = []
with open(prefix + ".grm.id") as f:
    for line in f:
        fid, iid = line.rstrip("\n").split()[:2]
        ids.append(iid)

n = len(ids)
K = np.eye(n, dtype=np.float32)

# 2. 读取 pairwise GRM
with gzip.open(prefix + ".grm.gz", "rt") as f:
    for line in f:
        i, j, _N, v = line.split()
        i = int(i) - 1
        j = int(j) - 1
        v = float(v)
        K[i, j] = v
        K[j, i] = v

# 3. 输出方阵
out = prefix + ".matrix.tsv"
with open(out, "w") as w:
    w.write("\t" + "\t".join(ids) + "\n")
    for i, name in enumerate(ids):
        w.write(name + "\t" + "\t".join(f"{x:.6g}" for x in K[i]) + "\n")

print("Done! Output:", out)
PY
——————————————————————————————————————————————————————————————————————————————————
#格式转换
python3  grm_to_matrix.py  
#最后拿到亲缘关系G矩阵
-rw-r--r-- 1 wangtao01 wangtao01     18181 Dec 21 01:26 gcta_kinship.matrix.tsv

gcta_kinship.matrix.tsv

• ③主成分分析PCA
  可使用plink直接分析，也可以用gcta生成亲缘关系再计算PCA

##1.plink
#--pca设置几个主成分，可设置样本数相同的数量
plink --bfile gwas.qc.chrnum \
     --pca 10 \
     --out plink_PCA  \
     --allow-extra-chr --set-missing-var-ids @:#
#输出主成分特征值和特征向量两个主要文件，eigenvec即包括10列主成分值
-rw-r--r-- 1 wangtao01 wangtao01        78 Dec 21 02:53 plink_PCA.eigenval
-rw-r--r-- 1 wangtao01 wangtao01      4999 Dec 21 02:53 plink_PCA.eigenvec

##2.gcta
#先进行亲缘关系计算。其实plink也做了，两步一起做了
#注意--make-grm参数而非--make-grm-gz
gcta --bfile gwas.qc.chrnum \
     --make-grm \
     --make-grm-alg 0 \
     --out gcta_PCA
#输出
-rw-r--r-- 1 wangtao01 wangtao01      3280 Dec 21 03:03 gcta_PCA.grm.bin
-rw-r--r-- 1 wangtao01 wangtao01       960 Dec 21 03:03 gcta_PCA.grm.id
-rw-r--r-- 1 wangtao01 wangtao01      3280 Dec 21 03:03 gcta_PCA.grm.N.bin

#再进行PCA分析
gcta --grm gcta_PCA \
     --pca 10 \
     --out gcta_PCA
#同样输出主成分特征值和特征向量两个主要文件
-rw-r--r-- 1 wangtao01 wangtao01       339 Dec 21 03:10 gcta_PCA.eigenval
-rw-r--r-- 1 wangtao01 wangtao01      5026 Dec 21 03:10 gcta_PCA.eigenvec

特征向量文件每一列代表每种pca值，特征值文件每种PCA变异量

在文献里找一个好看点的图，截给大模型，让其写一个用于可视化的R脚本pca_plot_plink.R，画出来效果不错，挺还原的。另外我这里画图的时候，我重新创建了一个R环境r_gwas,专门用于画图，防止污染我常用的base主环境。

############################################################
## PLINK PCA visualization (PC1 vs PC2)
## Environment: conda env "r_gwas"
############################################################

suppressPackageStartupMessages({
  library(data.table)
  library(dplyr)
  library(ggplot2)
})

## ---------- 1. read eigenvec ----------
eigvec <- fread("plink_PCA.eigenvec", header = FALSE)
# FID IID PC1 PC2 ...
setnames(eigvec, c("FID", "IID", paste0("PC", 1:(ncol(eigvec) - 2))))

## ---------- 2. read eigenval & explained variance ----------
eigval <- scan("plink_PCA.eigenval", quiet = TRUE)
var_exp <- eigval / sum(eigval)

pc1_lab <- sprintf("PC1 (%.1f%%)", 100 * var_exp[1])
pc2_lab <- sprintf("PC2 (%.1f%%)", 100 * var_exp[2])

## ---------- 3. read phenotype (sex) ----------
pheno <- fread("sex.pheno.txt", header = FALSE)
setnames(pheno, c("FID", "IID", "sex"))

## ⚠️ 按你的实际含义改这里
## 例如：1 = Female, 2 = Male（如果相反就调换）
pheno[, sex := factor(sex, levels = c(1, 2),
                      labels = c("Sex1", "Sex2"))]

## ---------- 4. merge ----------
df <- eigvec %>%
  inner_join(pheno, by = c("FID", "IID"))

## ---------- 5. PCA plot ----------
p <- ggplot(df, aes(x = PC1, y = PC2, color = sex)) +
  geom_hline(yintercept = 0, linetype = "dashed",
             linewidth = 0.4, color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed",
             linewidth = 0.4, color = "grey50") +
  geom_point(size = 2.6, alpha = 0.9) +
  stat_ellipse(
    aes(fill = sex),
    geom = "polygon",
    type = "norm",
    alpha = 0.18,
    color = NA
  ) +
  labs(
    title = "PLINK PCA (PC1 vs PC2)",
    x = pc1_lab,
    y = pc2_lab,
    color = "Sex",
    fill  = "Sex"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    panel.grid.minor = element_blank()
  )

## ---------- 6. save ----------
ggsave("PCA_PC1_PC2_sex.png", p, width = 7.2, height = 5.6, dpi = 300)
ggsave("PCA_PC1_PC2_sex.pdf", p, width = 7.2, height = 5.6)

message("Done: PCA_PC1_PC2_sex.png / PCA_PC1_PC2_sex.pdf")

• ④群体结构structure分析
  群体结构分析可以推断群体的分层关系，需要使用LD过滤后的数据gwas.qc.chrnum.LD进行分析。
  用如下for循环生成一个K2到K10的admixture脚本admixture.jobs.sh，后面用parallel并行运行。

#!/usr/bin/env bash
set -euo pipefail

PREFIX="gwas.qc.chrnum.LD"     # LD过滤后的bed前缀
K_MIN=2
K_MAX=10

# 输出命令列表
OUT_SH="admixture.jobs.sh"
: > ${OUT_SH}

for K in $(seq ${K_MIN} ${K_MAX}); do
  echo "admixture --cv=10 -j4 ${PREFIX}.bed ${K} 1> logs/admixture.K${K}.log 2>&1" >> ${OUT_SH}
done

echo "[OK] wrote ${OUT_SH}"
echo "Run: mkdir -p logs && parallel -j ${K_MAX} < ${OUT_SH}"

用parallel并行运行，得到Q矩阵，以及判断最合适（最小）K所用到的cv error值

mkdir -p logs
parallel -j 3 < admixture.jobs.sh

#输出Q矩阵，Q值代表来源某一亚群的可能性，所有亚群加起来等于1
-rw-r--r-- 1 wangtao01 wangtao01       720 Dec 21 14:15 gwas.qc.chrnum.LD.2.Q
-rw-r--r-- 1 wangtao01 wangtao01   4857804 Dec 21 14:15 gwas.qc.chrnum.LD.2.P
-rw-r--r-- 1 wangtao01 wangtao01      1440 Dec 21 14:16 gwas.qc.chrnum.LD.4.Q
-rw-r--r-- 1 wangtao01 wangtao01   9715608 Dec 21 14:16 gwas.qc.chrnum.LD.4.P
-rw-r--r-- 1 wangtao01 wangtao01      1080 Dec 21 14:17 gwas.qc.chrnum.LD.3.Q
-rw-r--r-- 1 wangtao01 wangtao01   7286706 Dec 21 14:17 gwas.qc.chrnum.LD.3.P
-rw-r--r-- 1 wangtao01 wangtao01      1800 Dec 21 14:18 gwas.qc.chrnum.LD.5.Q
-rw-r--r-- 1 wangtao01 wangtao01  12144510 Dec 21 14:18 gwas.qc.chrnum.LD.5.P
-rw-r--r-- 1 wangtao01 wangtao01      2160 Dec 21 14:19 gwas.qc.chrnum.LD.6.Q
-rw-r--r-- 1 wangtao01 wangtao01  14573412 Dec 21 14:19 gwas.qc.chrnum.LD.6.P
-rw-r--r-- 1 wangtao01 wangtao01      2880 Dec 21 14:21 gwas.qc.chrnum.LD.8.Q
-rw-r--r-- 1 wangtao01 wangtao01  19431216 Dec 21 14:21 gwas.qc.chrnum.LD.8.P
drwxr-xr-x 2 wangtao01 wangtao01      4096 Dec 21 14:21 logs
-rw-r--r-- 1 wangtao01 wangtao01      2520 Dec 21 14:25 gwas.qc.chrnum.LD.7.Q
-rw-r--r-- 1 wangtao01 wangtao01  17002314 Dec 21 14:25 gwas.qc.chrnum.LD.7.P
-rw-r--r-- 1 wangtao01 wangtao01      3240 Dec 21 14:26 gwas.qc.chrnum.LD.9.Q
-rw-r--r-- 1 wangtao01 wangtao01  21860118 Dec 21 14:26 gwas.qc.chrnum.LD.9.P
-rw-r--r-- 1 wangtao01 wangtao01      3600 Dec 21 14:37 gwas.qc.chrnum.LD.10.Q
-rw-r--r-- 1 wangtao01 wangtao01  24289020 Dec 21 14:37 gwas.qc.chrnum.LD.10.P

head gwas.qc.chrnum.LD.4.Q
0.515412 0.330016 0.000010 0.154562
0.000010 0.634053 0.006361 0.359577
0.000010 0.855060 0.000010 0.144920
0.000010 0.436062 0.000010 0.563918
0.000010 0.000010 0.000010 0.999970
0.000010 0.822138 0.002708 0.175144
0.000010 0.879779 0.000010 0.120201
0.000010 0.808053 0.191927 0.000010
0.000010 0.862587 0.137393 0.000010
0.000010 0.816857 0.029244 0.153889

#提取log文件里的cv error值，最小的那个对应着最合适的k值，这里结果K为2的时候最小
grep -h "CV error" logs/admixture.K*.log
CV error (K=10): 1.19629
CV error (K=2): 0.56840
CV error (K=3): 0.65797
CV error (K=4): 0.69925
CV error (K=5): 0.78295
CV error (K=6): 0.85496
CV error (K=7): 0.92420
CV error (K=8): 1.02147
CV error (K=9): 1.10976

同样切到r_gwas环境，用genk里的一个R脚本进行画图，这里调用的是pophelper包。

at draw_admixtureG.R
#!/usr/bin/env Rscript
# parse parameter ---------------------------------------------------------
library(argparser, quietly=TRUE)
# Create a parser
p <- arg_parser("draw sturcture figure for admixture result")

# Add command line arguments
p <- add_argument(p, "Qlist", help="input: Q matirx list file", type="character")
p <- add_argument(p, "sample", help="input: samplefile, ie  .fam", type="character")
p <- add_argument(p, "sample_order", help="input: sample order file with group, define order in output figure ", type="character")
p <- add_argument(p, "output", help="output prefix", type="character")

# Parse the command line arguments
argv <- parse_args(p)

sfiles <- argv$Qlist
samp <- argv$sample
ord  <- argv$sample_order
outpre <- argv$output


#sfiles <- "./demo.Q.list"
#samp <- "../00.prepare/demo.LD.fam"
#ord <- "../data/sample.pop"
#outpre <- "test"


library(pophelper)
library(ggplot2)


Qfiles <- read.table(sfiles, header = F)
qlist <- readQ(Qfiles$V1, filetype = "basic", indlabfromfile=F)

label <- read.table( samp, header = F,stringsAsFactors=F )

ordinfo <- read.table(ord, header = F, stringsAsFactors=F )
rownames(ordinfo) <- ordinfo$V1

for( i in 1:length(qlist) ){
    rownames(qlist[[i]]) <- label$V1
    qlist[[i]] <- qlist[[i]][as.character(ordinfo[,1]),]
}

# head(qlist[[3]])

xlist <- alignK(qlist)

## 每个K分别排序绘图
p1_sort  <- plotQ(
      xlist,
      sortind = "all",  ## ind每个K分别排序 
      imgoutput = "join",  ## 一张图还是多张
      returnplot=T,  
      exportplot=F,  
      useindlab = T ,    ## 显示ind名称
      sharedindlab= F ,  ## ind出现一次
      showindlab=T ,
      splab=paste0("K=",sapply(xlist,ncol))
      )

ggsave(filename = paste(outpre, "sorted.pdf",sep=".") , 
       p1_sort$plot[[1]],  
       width = 20,
       height = 10  
       )

## 按指定顺序绘图
p1_order  <- plotQ(
  xlist, 
  sortind = NA,  ## ind排序 
  imgoutput = "join",  ## 一张图还是多张
  returnplot=T,  
  exportplot=F,
  useindlab = T ,    ## 显示ind名称
  sharedindlab= T ,   ## ind出现一次
  showindlab=T,
  splab=paste0("K=",sapply(xlist,ncol)), ## 只显示K值
  grplab = ordinfo[,2,drop=FALSE], ## 添加grplab
  grplabsize=4,linesize=0.8,pointsize=4,
  ordergrp = F,  # 不安装grp排序
  )

ggsave(filename = paste(outpre, "ordered.pdf",sep=".") , 
       p1_order$plot[[1]],  
       width = 20,
       height = 10  
)



## 按指顺序
p1_group  <- plotQ(
  xlist, 
  sortind = NA,  ## ind排序 
  imgoutput = "join",  ## 一张图还是多张
  returnplot=T,  
  exportplot=F,
  useindlab = T ,    ## 显示ind名称
  sharedindlab= T ,   ## ind出现一次
  showindlab=T,
  splab=paste0("K=",sapply(xlist,ncol)), ## 只显示K值
  grplab = ordinfo[,2,drop=FALSE],
  grplabsize=4,linesize=0.8,pointsize=4,
  ordergrp = T)

ggsave(filename = paste(outpre, "orderedgroup.pdf",sep=".") , 
       p1_group$plot[[1]],  
       width = 20,
       height = 10  
)

这个脚本对你数据的输入要求（需要准备 3 个输入文件）
它的命令行参数是：
1.Qlist：一个文件，里面每行一个 .Q 文件路径
2.sample：fam 文件（gwas.qc.chrnum.LD.fam）
3.sample_order：一个两列文件：样本ID+分组名
4.output：输出前缀

#生成一个Q 文件路径列表
ls gwas.qc.chrnum.LD.{2,3,4,5}.Q > gwas.Q.list

#fam文件已存在

#手写或者用awk命令准备一个sample_order文件
awk 'BEGIN{OFS="\t"}
     {g=($3==1?"Male":($3==2?"Female":"Unknown")); print $1,g}' sex.pheno.txt \
| sort -k2,2 -k1,1 > sample.pop.order


##运行
Rscript 2.structure_plot_admixture.R \
  gwas.Q.list \
  gwas.qc.chrnum.LD.fam \
  sample.pop.order \
  structure
#输出3个pdf文件
-rw-r--r--  1 wangtao01 wangtao01       9780 Dec 22 00:15 structure.sorted.pdf
-rw-r--r--  1 wangtao01 wangtao01       9477 Dec 22 00:15 structure.ordered.pdf
-rw-r--r--  1 wangtao01 wangtao01       9477 Dec 22 00:15 structure.orderedgroup.pdf

• ⑤连琐不平衡衰减分析LD decay
  LD 衰减分析可以查看整个群体及亚群的变异位点间的连锁水平。
  软件：PopLDdecay ；conda安装的找不到画图脚本，所以去GitHub重新下载安装了。
  官网：https://github.com/BGI-shenzhen/PopLDdecay
  分析还是分两步：1.PopLDdecay程序计算群体内两两位点间的r^2值；2. Plot_OnePop.pl画图。如果有多个群体就分别进行r^2值计算，再用Plot_MultiPop.pl一起可视化

#1.计算位点间的r^2值
#分亚群的先生成一个亚群的样本列表文件，然后加参数 -SubPop
~/software/PopLDdecay/bin/PopLDdecay \
  -InVCF  gwas.qc.vcf \
  -SubPop female.list \
  -MaxDist 500 \
  -OutStat pop.Female.LDdecay   

~/software/PopLDdecay/bin/PopLDdecay \
  -InVCF  gwas.qc.vcf \
  -SubPop male.list \
  -MaxDist 500 \
  -OutStat pop.Male.LDdecay   

#输出.stat.gz结果
-rw-r--r--  1 wangtao01 wangtao01    5633612 Dec 23 03:15 pop.Female.LDdecay.stat.gz
-rw-r--r--  1 wangtao01 wangtao01    5639255 Dec 23 03:15 pop.Male.LDdecay.stat.gz

gzip -dc pop.Female.LDdecay.stat.gz |head -n 10
#Dist   Mean_r^2        Mean_D' Sum_r^2 Sum_D'  NumberPairs
1       0.7804  NA      47608.7121      NA      61006
2       0.7489  NA      26949.2520      NA      35983
3       0.7270  NA      23488.2912      NA      32310
4       0.7199  NA      22281.5461      NA      30949
5       0.7109  NA      21016.1581      NA      29561
6       0.7025  NA      21004.5115      NA      29901
7       0.6942  NA      19522.0154      NA      28123
8       0.6894  NA      19483.0819      NA      28259
9       0.6822  NA      19381.8086      NA      28411


#2. 单个群体或者亚群的画图
~/software/PopLDdecay/bin/Plot_OnePop.pl \
  -inFile pop.Female.LDdecay.stat.gz  \
  -output pop.Female.LDdecay \
  -keepR #保留画图的R脚本，后续图丑可以进行调整

#3. 多个群体或者亚群的画图
#先准备一个所以亚群LD分析结果文件列表Pop.ResultPath.list 
ls pop.*.stat.gz |awk -F"." '{ print $0"\t"$2 }' > Pop.ResultPath.list 

cat Pop.ResultPath.list
pop.Female.LDdecay.stat.gz      Female
pop.Male.LDdecay.stat.gz        Male

#画图
~/software/PopLDdecay/bin/Plot_MultiPop.pl \
  -inList  Pop.ResultPath.list  \
  -output pop.all.LDdecay \
  -keepR

• ⑥连琐不平衡区块LD block分析
  全基因组水平的LD block分析用plink计算，计算较慢，同时也可设置--chr --from-bp --to-bp三个参数来限定染色体特定区域的LD block分析；染色体局部区域的连锁情况使用LDBlockShow展示出图，-SeleVar参数选择d'或r^2值作为LD值。

##1 plink 计算全基因组LD block，另外--chr --from-bp   --to-bp 3个参数用来计算染色体局部连琐情况
plink --vcf gwas.qc.vcf \
        --out plink_ldblock \
        --chr chr1A --from-bp 1  --to-bp 500000 \ 
        --blocks no-pheno-req \
        --allow-extra-chr  --set-missing-var-ids @:#
#输出
-rw-r--r--  1 wangtao01 wangtao01        960 Dec 24 15:29 plink_ldblock.nosex
-rw-r--r--  1 wangtao01 wangtao01       4725 Dec 24 15:29 plink_ldblock.blocks.det
-rw-r--r--  1 wangtao01 wangtao01       2419 Dec 24 15:29 plink_ldblock.blocks
-rw-r--r--  1 wangtao01 wangtao01       1180 Dec 24 15:29 plink_ldblock.log

(main) [wangtao01@sonmi:~/gwas/tana/gwas_workdir]$ head plink_ldblock.blocks.det
 CHR          BP1          BP2           KB  NSNPS SNPS
chr1A        53391        55361        1.971     11 chr1A:53391|chr1A:53434|chr1A:53517|chr1A:53930|chr1A:53950|chr1A:54060|chr1A:54198|chr1A:54559|chr1A:54913|chr1A:54941|chr1A:55361
chr1A        58187        59474        1.288      6 chr1A:58187|chr1A:58477|chr1A:58478|chr1A:59167|chr1A:59283|chr1A:59474
chr1A        59600        59720        0.121      2 chr1A:59600|chr1A:59720
chr1A        60573        60605        0.033      2 chr1A:60573|chr1A:60605
chr1A        71835        72613        0.779      4 chr1A:71835|chr1A:72181|chr1A:72233|chr1A:72613
chr1A        73471        77176        3.706      3 chr1A:73471|chr1A:77097|chr1A:77176
chr1A        83372        90234        6.863      9 chr1A:83372|chr1A:84389|chr1A:84390|chr1A:86049|chr1A:86260|chr1A:87214|chr1A:87230|chr1A:88741|chr1A:90234
chr1A       107186       108656        1.471      2 chr1A:107186|chr1A:108656
chr1A       111987       112416         0.43      3 chr1A:111987|chr1A:112165|chr1A:112416

##2 LDBlockShow 用来进行染色体局部LDblock的可视化，-SeleVar 用来指定LD值计算方法（1: D'   2: r^2）
#另外默认输出svg图片，-OutPdf  -OutPng参数可用来转化，但默认会调用ImageMagick软件，其依赖GLIBCXX较新版本，本服务器不支持，可用其他软件如rsvg-convert替代。
#谨慎解决ImageMagick软件报错问题，我吃了大亏，解决依赖冲突的时候把我base环境污染了，conda都用不了，得重新装miniconda/miniforge。
LDBlockShow  -InVCF gwas.qc.vcf \
   -OutPut  LDBlockShow \
   -Region chr1A:300000-500000 \
   -SeleVar 2 \
   -OutPdf   -OutPng
#rsvg-convert软件转换图片
rsvg-convert -o LDBlockShow.png LDBlockShow.svg
rsvg-convert -f pdf -o LDBlockShow.pdf LDBlockShow.svg