1. 軟件下載及安裝admixture:
使用conda進行軟件安裝
conda install?admixture
2. VCF文件格式轉(zhuǎn)換為bed格式文件(似乎admixture 可以直接識別ped/map文件格式的輸入文件)
vcf文件轉(zhuǎn)為ped文件:
方法1:
使用vcftools支持將vcf文件轉(zhuǎn)換成plink對應的ped/map格式,如下
vcftools? --vcf input.vcf --plink --out output
方法2:
plink支持直接讀取vcf文件格式常柄,基本用法如下:
plink --vcf input.vcf --recode --out output?
重要端幼!?因為轉(zhuǎn)換成的ped和map文件無法匹配邻眷,需要手動更改上一步轉(zhuǎn)換好的map文件
map數(shù)據(jù)格式為四列:
使用plink將ped/map轉(zhuǎn)換為二進制的bed文件代咸,命令行如下:
plink --file inputfile --make-bed --out filename
第一個FILENAME的后綴為.ped和.map,生成的第二個FILENAME的后綴為.bed陵吸、.bim玻墅、.fam
3.1.vcftools去除或保留vcf文件中的樣品
例1:只保留1和10號兩個樣品,執(zhí)行以下代碼:
vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--indv ?1--indv ?10 ?> out.vcf
例2:刪除1號樣品壮虫,執(zhí)行以下代碼:
vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--remove-indv ?1?> out.vcf
例3:如果樣品較多澳厢,也可將樣品保存到文件 id.txt 中环础,每行為一個樣品ID,格式如下:
sample1
2
..
然后使用下面兩個選項對vcf文件保留或者刪除樣品剩拢。
--keep<filename>保留樣品
--remove
<filename> ??刪除樣品
代碼如下:
vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--keep id.txt ? > out.vcf
作者:花事Le
鏈接:http://www.reibang.com/p/542d9b63dcd1
來源:簡書
著作權(quán)歸作者所有线得。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權(quán),非商業(yè)轉(zhuǎn)載請注明出處徐伐。
3.2 plink提取指定樣本和指定SNP的數(shù)據(jù)(keep贯钩,extract函數(shù)
plink --bfile inputfile --noweb --keep sampleID.txt --recode --make-bed --out fileout
inputfile為不加.bed后綴的bed文件
其中,sampleID.txt第一列為提取的樣本Family ID办素,第二列為Within-family ID(IID)
plink提取SNP位點:
plink --bfile file --extract snp.txt --make-bed --out snp
其中角雷,snp.txt的文件格式如下,一個SNP位點一行:
rs1
rs2
rs3
4. 如何選擇合適的K值
可以同時運行多個程序, 每個程序不同的k值, 比如, 想要k值選擇1,2,3,4,5, 可以寫為:
?for?K?in?1?2?3?4?5;?do?admixture?--cv?hapmap3.bed?$K?|?tee?log${K}.out;?done
例子:
for K in 1 2 3 4 5 6 7 8 9 10 11 12; do admixture --cv 10729bed2.bed $K | tee log${K}.out; done
多線程: admixture??hapmap3.bed?3?-j?4
使用grep命令去查看*out文件的cv error(交叉驗證的誤差)值:
grep?-h?CV??*.out
結(jié)果如下:(這個K值顯示是否有誤?應該從第一開始分別是K=1性穿,2勺三,3依次往下)
5. 繪制Q值的百分比柱狀圖
使用R語言
ta1?=?read.table("D:/files.3.Q")
head(ta1)
barplot(t(as.matrix(ta1)),col?=?rainbow(3),
????????xlab?=?"Individual",
????????ylab?=?"Ancestry",
????????border?=?NA)
————————————————————————————————————————————
本文部分分析步驟參考了CSDN博主「育種數(shù)據(jù)分析之放飛自我」的原創(chuàng)文章,遵循CC 4.0 BY-SA版權(quán)協(xié)議绊含,轉(zhuǎn)載請附上原文出處鏈接及本聲明桑嘶。?
原文鏈接:https://blog.csdn.net/yijiaobani/article/details/83017730