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)行和上述一樣的分析即可

參考

--Admixture使用說(shuō)明文檔cookbook