1. 軟件下載及安裝admixture：

使用conda進(jìn)行軟件安裝

conda install?admixture

2. VCF文件格式轉(zhuǎn)換為bed格式文件（似乎admixture 可以直接識(shí)別ped/map文件格式的輸入文件）

vcf文件轉(zhuǎn)為ped文件：

方法1：

使用vcftools支持將vcf文件轉(zhuǎn)換成plink對(duì)應(yīng)的ped/map格式，如下

vcftools? --vcf input.vcf --plink --out output

方法2：

plink支持直接讀取vcf文件格式，基本用法如下：

plink --vcf input.vcf --recode --out output?

map文件? ? ? 染色體編號(hào)為數(shù)字, 未知為0SNP名稱為字符或數(shù)字, 如果不重要, 可以從1編號(hào), 注意要和bed文件SNP列一一對(duì)應(yīng)染色體的摩爾未知(可選項(xiàng), 可以用0)SNP物理坐標(biāo)

重要！?因?yàn)檗D(zhuǎn)換成的ped和map文件無(wú)法匹配，需要手動(dòng)更改上一步轉(zhuǎn)換好的map文件

map數(shù)據(jù)格式為四列：

bed文件? ? 第一列: Family ID # 如果沒(méi)有, 可以用個(gè)體ID代替第二列: Individual ID # 個(gè)體ID編號(hào)第三列: Paternal ID # 父本編號(hào)第四列: Maternal ID # 母本編號(hào)第五列: Sex (1=male; 2=female; other=unknown) # 性別, 如果未知, 用0表示第六列: Phenotype # 表型數(shù)據(jù), 如果未知, 用0表示第七列以后: 為SNP分型數(shù)據(jù), 可以是AT CG或11 12, 或者A T C G或1 1 2 2————————————————版權(quán)聲明：本文為CSDN博主「育種數(shù)據(jù)分析之放飛自我」的原創(chuàng)文章，遵循CC 4.0 BY-SA版權(quán)協(xié)議，轉(zhuǎn)載請(qǐng)附上原文出處鏈接及本聲明。原文鏈接：https://blog.csdn.net/yijiaobani/article/details/83017730

使用plink將ped/map轉(zhuǎn)換為二進(jìn)制的bed文件，命令行如下：

plink --file inputfile --make-bed --out filename

第一個(gè)FILENAME的后綴為.ped和.map，生成的第二個(gè)FILENAME的后綴為.bed、.bim、.fam

3.1.vcftools去除或保留vcf文件中的樣品

例1：只保留1和10號(hào)兩個(gè)樣品，執(zhí)行以下代碼：

vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--indv ?1--indv ?10 ?> out.vcf

例2：刪除1號(hào)樣品，執(zhí)行以下代碼：

vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--remove-indv ?1?> out.vcf

例3：如果樣品較多，也可將樣品保存到文件 id.txt 中，每行為一個(gè)樣品ID，格式如下：

sample1

2

..

然后使用下面兩個(gè)選項(xiàng)對(duì)vcf文件保留或者刪除樣品。

--keep<filename>保留樣品

--remove

<filename> ??刪除樣品

代碼如下：

vcftools --vcf in.vcf --recode --recode-INFO-all --stdout ?--keep id.txt ? > out.vcf

作者：花事Le

鏈接：http://www.itdecent.cn/p/542d9b63dcd1

來(lái)源：簡(jiǎn)書(shū)

著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處。

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位點(diǎn)：

plink --bfile file --extract snp.txt --make-bed --out snp

其中，snp.txt的文件格式如下,一個(gè)SNP位點(diǎn)一行：

rs1

rs2

rs3

4. 如何選擇合適的K值

可以同時(shí)運(yùn)行多個(gè)程序, 每個(gè)程序不同的k值, 比如, 想要k值選擇1,2,3,4,5, 可以寫(xiě)為:

?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(交叉驗(yàn)證的誤差)值:

grep?-h?CV??*.out

結(jié)果如下：（這個(gè)K值顯示是否有誤？應(yīng)該從第一開(kāi)始分別是K=1，2，3依次往下）

對(duì)這個(gè)K值出現(xiàn)這樣的情況？為何K10開(kāi)始，個(gè)人覺(jué)得這個(gè)K值顯示有誤，應(yīng)該從第一開(kāi)始分別是K=1，2，3依次往下

5. 繪制Q值的百分比柱狀圖

使用R語(yǔ)言

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)載請(qǐng)附上原文出處鏈接及本聲明。?

原文鏈接：https://blog.csdn.net/yijiaobani/article/details/83017730