对 SNP 进行 LD 过滤

vcf=all.vcf
plink \
  --vcf ${vcf} \
  --geno 0.1 \
  --maf 0.01 \
  --biallelic-only strict \
  --out all.missing_maf \
  --recode vcf-iid \
  --allow-extra-chr \
  --set-missing-var-ids @:# \
  --keep-allele-order
vcftools \
  --vcf ${vcf} \
  --maf 0.01 \
  --max-missing 0.9 \
  --min-alleles 2 \
  --max-alleles 2 \
  --recode \
  --recode-INFO-all \
  --out all.missing_maf
plink \
  --vcf all.missing_maf.vcf \
  --extract tmp.ld.prune.in \
  --out all.LDfilter \
  --recode vcf-iid \
  --keep-allele-order \
  --allow-extra-chr \
  --set-missing-var-ids @:#

进化树构建

将过滤后的vcf文件，转换为 phylip 文件，再用多种软件进行建树。

fasttree
iqtree
raxml

NJ 法构树

fasttree

iqtree

raxml-ng

PCA 分析

plink

对过滤后的vcf文件进行分析。

gcta

eigensoft

structure 分析

admixture
structure

CJchen | 群体结构分析 | Pophelper 的“平替版”

plink \
  --vcf all.LDfilter.vcf \
  --make-bed \
  --out all \
  --allow-extra-chr \
  --keep-allele-order \
  --set-missing-var-ids @:#
for i in `seq 2 4`
do
admixture \
  --cv \
  -j2 \
  all.bed \
  $i \
  1> admix.$i.log 2>&1
done
# seq 2 4 \
#   | awk '{print "admixture --cv -j2 all.bed "$1" 1>admix."$1".log 2>&1"}' \
#   > admixture.sh

LDdecay 分析

PopLDdecay \
  -InVCF ${vcf} \
  -SubPop pop.SC.table \
  -MaxDist 500 \
  -OutStat pop.SC.stat
PopLDdecay \
  -InVCF ${vcf} \
  -SubPop pop.YZR.table \
  -MaxDist 500 \
  -OutStat pop.YZR.stat

LDblock 分析

# 提取需要展示区域的 vcf 文件
awk '$1~/ˆ#/ || ($1=="1" && $2 > 300000 && $2 < 500000)' all.missing_maf.vcf > chr1_300K-500K.vcf
# 使用 plink 计算两两 SNP R^2 值
plink \
  --vcf chr1_300K-500K.vcf \
  --r2 square \
  --out chr1_300K-500K.square \
  --allow-extra-chr
# 生成 SNP 坐标表格
# 绘制 LD 热图
LDBlockShow \

参考