首先,对 vcf 文件的格式进行转换,转换为 plink 格式:
vcftools --gzvcf test.vcf.gz --plink --out test_out
而后根据连锁不平衡对数据进行过滤:
plink --noweb --file test_out --indep-pairwise 50 10 0.1 --geno 0 --out test_out
plink --file test_out --make-bed --extract test_out.prune.in --chr-set 12 no-xy --out test_out2
plink --bfile test_out2 --recode 12 --out test_out3
基于上述文件,首先计算 admixture:
for i in {2..8}; do admixture --cv test_out3.ped ${i} | tee log${i}.out; done
然后使用 gcta 进行 PCA 分析
gcta64 --make-grm --out test_PCA --bfile test_out2 --autosome-num 12
gcta64 --grm test_PCA --pca 3 --out test_PCA_out
最后构建 NJ 树:
plink --file Nip_Low_Depth3 --distance-matrix --out Nip_distance2
linux
vcftools --vcf 241218_235_samples_filter.recode.vcf --chr chr03 --from-bp 13231683 --to-bp 13252880 --recode --out ino80
vcftools --vcf ino80.recode.vcf --keep Niv2_Ruf2.list --recode --out ino80_Niv2
plink --vcf ino80_Niv2.recode.vcf --make-bed --out ino80 --double-id
plink --bfile ino80 --recode 12 --out ino80
plink --file ino80 --distance-matrix --out ino80
R
dist_matrix <- read.table("ino80.mdist", header = T)
headers <- read.table("ino80.mdist.id", header = FALSE, stringsAsFactors = FALSE)
colnames(dist_matrix) <- headers$V1
rownames(dist_matrix) <- headers$V1
dist_matrix <- dist_matrix[complete.cases(dist_matrix), complete.cases(t(dist_matrix))]
dist_obj <- as.dist(dist_matrix)
nj_tree <- nj(dist_obj)
write.tree(nj_tree, file = "ino80_nj_tree.nwk")