tructure是與PCA、進(jìn)化樹相似的方法饰及，就是利用分子標(biāo)記的基因型信息對(duì)一組樣本進(jìn)行分類蔗坯，分子標(biāo)記可以是SNP、indel燎含、SSR宾濒。相比于PCA，進(jìn)化樹屏箍，群體結(jié)構(gòu)分析可明確各個(gè)群之間是否存在交流及交流程度
1 軟件安裝
conda install -c bioconda admixture
admixture
**** ADMIXTURE Version 1.3.0 ****
**** Copyright 2008-2015 ****
**** David Alexander, Suyash Shringarpure, ****
**** John Novembre, Ken Lange ****
**** ****
**** Please cite our paper! ****
**** Information at www.genetics.ucla.edu/software/admixture ****
Usage: admixture <input file> <K>
See --help or manual for more advanced usage.
2 簡(jiǎn)單使用
第一步：將VCF變?yōu)閜link格式
## vcftools
vcftools --gzvcf SNP.vcf.gz --plink --out test
## 沒(méi)有壓縮绘梦，則為--vcf
## 或者plink 直接轉(zhuǎn)
plink --vcf SNP.vcf.gz--recode --out test--const-fid --allow-extra-chr
# --vcf vcf 或者vcf.gz
# --recode 輸出格式ped（默認(rèn)bed）
# --out 輸入前綴
# --const-fid 添加群體信息familyID sampleID）(將family設(shè)置為0)橘忱；
# --allow-extra-chr 允許非標(biāo)準(zhǔn)染色體編號(hào)
plink轉(zhuǎn)化，map文件卸奉，SNP那列為點(diǎn),vcftools 則不是钝诚，但是ref為0
# 如果用vcftools轉(zhuǎn)換, 重新添加染色體
paste <(cut -f2 test.map |awk -F ":" '{print $1}') <(cut -f2-4 test.map) >map
## 如果用plink 轉(zhuǎn)換, 重新添加SNP編號(hào)
awk '{x+=1}{print $1"\tSNP"x"\t"$3"\t"$4}' Test.map >map
#重命名
mv map Test.map
因?yàn)閜link本身是針對(duì)人類進(jìn)行開發(fā)的，所以在運(yùn)行時(shí)榄棵，加上--allow-extr-chr凝颇。此外對(duì)于vcf中的sampleID（familyID_sampleID）, plink 默認(rèn)為下劃線分隔（也可以通過(guò)參數(shù)--id-delim進(jìn)行修改），分別作為family ID和sampleID疹鳄。但是一般我們的樣本并不是那樣命名的拧略，所以可以添加--double-id參數(shù)，將familyID和sampleID命名為相同瘪弓，或者--const-fid,將familyID命名為0鳍徽，表明-9
第二步 plink進(jìn)行過(guò)濾淆两，得到bed文件
plink --allow-extra-chr --noweb -file test --geno 0.05 --maf 0.05 --hwe 0.0001 --make-bed --out test1
# --noweb 不顯示網(wǎng)頁(yè)
第三步：尋找合適的K值
for i in {1..7};do admixture --cv test1.bed $i |tee log${i}.out;done
## 根據(jù)CV error 確定K值
grep -h 'CV' log*.out
CV error (K=1): 0.22317
CV error (K=2): 0.15018
CV error (K=3): 0.12804
CV error (K=4): 0.12109
CV error (K=5): 0.12656
當(dāng)k=4時(shí)尤辱，cv error 最小迹淌，選擇4
此外齿诉，
如果SNP數(shù)據(jù)集非常大流部，則可以隨機(jī)選擇SNP進(jìn)行K值選擇分析摹蘑，比如隨機(jī)選取20000個(gè)SNP進(jìn)行分析惊畏，每個(gè)K值跑20次舀透，確定最終的k值栓票，然后分析
當(dāng)利用plink轉(zhuǎn)的格式中，在運(yùn)行上述命令是出現(xiàn)以下報(bào)錯(cuò)
Invalid chromosome code! Use integers
將*.bim中的第一列改為數(shù)值就可以了
第四步：多線程
admixture test1.bed 4 -j 5
## j愕够， 線程數(shù)
第五步：作圖
ta1 = read.table("test1.4.Q")
head(ta1)
barplot(t(as.matrix(ta1)),col = rainbow(3),
xlab = "Individual",
ylab = "Ancestry",
border = NA)
根據(jù)LD進(jìn)行篩選SNP
首先 篩選合格的LD位點(diǎn)
## 對(duì)vcf進(jìn)行
plink --vcf SNP.vcf.gz --indep-pairwise 100 50 0.2 -out Test-id-maf0.05-LD --allow-extra-chr
## 或者對(duì)bed進(jìn)行操作也可以
plink --bfile test1 --indep-pairwise 100 50 0.2
## --indep-pairwise 100 50 0.2走贪； 100Kb窗口，50步長(zhǎng)惑芭，0.2LD強(qiáng)度
然后進(jìn)行提取
plink --bfile test1 --extract plink.prune.in --make-bed --out test2
然后在進(jìn)行和上述一樣的分析即可
