By Tingting, 16 September, 2024
Forums

由于又对原文件进行处理,之前操作的步骤笔记记录的比较分散,所以整理汇总一部分。

 

文件:combine_all.pass.recode.clear.vcf

 一、删除突变位点(--snps-only),只保留染色体上的SNP(--chr 1-13)

plink --vcf combine_all.pass.recode.clear.vcf --snps-only --chr 1-13 --make-bed --out combine_all.pass.recode.clear.snp.only.13chr.vcf

生成后缀为bed fam nosex log文件

二、将vcf文件转换为HapMap格式,后缀为hmp.txt,使用tassel

#在conda中安装
conda install -c biobuilds tassel
#将二进制文件转换为vcf文件
plink --bfile combine_all.pass.recode.clear.snp.only.13chr.vcf --recode vcf --out combine_all.pass.recode.clear.snp.only.13chr.vcf
#将vcf文件转化hapmap
run_pipeline.pl -Xms10g -Xmx100g  -vcf combine_all.pass.recode.clear.snp.only.13chr.vcf.vcf -sortPositions -export combine_all.pass.recode.snp.only.13chr.hmp.txt -exportType HapmapDiploid

三、群体结构分析

#conda下载admixture
conda install admixture
#运行admixture
for K in 1 2 3 4 5 7 8 9 10; do admixture --cv combine_all.pass.recode.clear.snp.only.13chr.vcf.bed $K | tee log${K}.out; done
#查看结果
grep -h CV log*.out
#查看最佳K值文件,输出最佳K值文件:3.Q
#在文件中第一列添加ID名字,以及表头Q1、Q2、Q3
awk 'NR==FNR{a[NR]=$0; next} {print a[FNR], $0}' ID.txt combine_all.pass.recode.clear.snp.only.13chr.vcf.3.Q > CV_3.txt
#在R中需检查文件的格式,转换为数字格式才可运行