Admixture做群体结构分析是好用,但也有一些不顺手的地方。最大的问题是不支持非整数的染色体号!
相信我们手里绝大部分vcf或plink格式文件,染色体ID基本是文本类型的吧。注意plink处理时加上-allow-extra-chr
,若染色体数超过人类,可使用--chr-set
设置。
scaffold或chr等字符不谈,我们可以手动处理为数字。但异常多倍体呢?比如小麦、油菜等,必须有chrA、ChrB。能否把亚组基因组也直接改为其他数字来代替?
官方好像没找到回答,一些网友认为是可以的。
比如异源四倍体的油菜,将A亚基因组改为数字,C亚基因组改为另一组数字。
sed 's/^chrA0//g;s/^chrA//g' snp.bim |sed 's/^chrC01/11/g;s/^chrC02/12/g;s/^chrC03/13/g;s/^chrC04/14/g;s/^chrC05/15/g;s/^chrC06/16/g;s/^chrC07/17/g;s/^chrC08/18/g;s/^chrC09/19/g' >new.snp.bim
mv new.snp.bim snp.bim
for K in {2..10};do
admixture --cv snp.bed $K | tee log${K}.out;
done
grep -h 'CV' log*.out >CV.log
网友评论