2023年7月10日发(作者:)
群体结构分析软件admixture安装及使⽤经验1. 软件下载及安装admixture:使⽤conda进⾏软件安装conda install admixture2. VCF⽂件格式转换为bed格式⽂件(似乎admixture 可以直接识别ped/map⽂件格式的输⼊⽂件)vcf⽂件转为ped⽂件:⽅法1:使⽤vcftools⽀持将vcf⽂件转换成plink对应的ped/map格式,如下vcftools --vcf --plink --out output⽅法2:plink⽀持直接读取vcf⽂件格式,基本⽤法如下:plink --vcf --recode --out output
map⽂件 染⾊体编号为数字, 未知为0SNP名称为字符或数字, 如果不重要, 可以从1编号, 注意要和bed⽂件SNP列⼀⼀对应染⾊体的摩尔未知(可选项, 可以⽤0)SNP物理坐标重要! 因为转换成的ped和map⽂件⽆法匹配,需要⼿动更改上⼀步转换好的map⽂件map数据格式为四列: 使⽤plink将ped/map转换为⼆进制的bed⽂件,命令⾏如下:plink --file inputfile --make-bed --out filename第⼀个FILENAME的后缀为.ped和.map,⽣成的第⼆个FILENAME的后缀为.bed、.bim、.fam
3. plink提取指定样本和指定SNP的数据(keep,extract函数plink --bfile inputfile --noweb --keep --recode --make-bed --out fileoutinputfile为不加.bed后缀的bed⽂件其中,第⼀列为提取的样本Family ID,第⼆列为Within-family ID(IID)plink提取SNP位点:plink --bfile file --extract --make-bed --out snp其中,的⽂件格式如下,⼀个SNP位点⼀⾏:rs1rs2rs34. 如何选择合适的K值可以同时运⾏多个程序, 每个程序不同的k值, ⽐如, 想要k值选择1,2,3,4,5, 可以写为: for K in 1 2 3 4 5; do admixture --cv $K | tee log${K}.out; done例⼦:for K in 1 2 3 4 5 6 7 8 9 10 11 12; do admixture --cv $K | tee log${K}.out; done多线程: admixture 3 -j 4使⽤grep命令去查看*out⽂件的cv error(交叉验证的误差)值:grep -h CV *.out结果如下:(这个K值显⽰是否有误?应该从第⼀开始分别是K=1,2,3依次往下)对这个K值出现这样的情况?为何K10开始,个⼈觉得这个K值显⽰有误,应该从第⼀开始分别是K=1,2,3依次往下5. 绘制Q值的百分⽐柱状图使⽤R语⾔ta1 = ("D:/files.3.Q")head(ta1)barplot(t((ta1)),col = rainbow(3), xlab = "Individual", ylab = "Ancestry", border = NA)
————————————————————————————————————————————本⽂部分分析步骤参考了CSDN博主「育种数据分析之放飞⾃我」的原创⽂章,遵循CC 4.0 BY-SA版权协议,转载请附上原⽂出处链接及本声明。
发布者:admin,转转请注明出处:http://www.yc00.com/web/1688984642a191586.html
评论列表(0条)