admixture 群體結構分析

tructure是與PCA、進化樹相似的方法，就是利用分子標記的基因型信息對一組樣本進行分類，分子標記可以是SNP、indel、SSR。相比於PCA，進化樹，群體結構分析可明確各個群之間是否存在交流及交流程度

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 簡單使用

第一步：將VCF變為plink格式

## vcftools 
vcftools --gzvcf SNP.vcf.gz --plink --out test
## 沒有壓縮，則為--vcf

## 或者plink 直接轉
plink --vcf SNP.vcf.gz--recode --out test--const-fid --allow-extra-chr
  
 # --vcf vcf 或者vcf.gz
 # --recode 輸出格式ped（默認bed）
 # --out 輸入前綴
 # --const-fid  添加群體信息familyID sampleID）(將family設置為0)；
 # --allow-extra-chr 允許非標准染色體編號

$\color{red}{**}$plink轉化，map文件，SNP那列為點,vcftools 則不是，但是ref為0

# 如果用vcftools轉換, 重新添加染色體
paste <(cut -f2 Test.map |awk -F ":" '{print $1}') <(cut -f2-4 Apple.map)  >map
## 如果用plink 轉換, 重新添加SNP編號
awk '{x+=1}{print $1"\tSNP"x"\t"$3"\t"$4}' Test.map >map
#重命名
mv map Test.map

因為plink本身是針對人類進行開發的，所以在運行時，加上--allow-extr-chr。此外對於vcf中的sampleID（familyID_sampleID）, plink 默認為下划線分隔（也可以通過參數--id-delim進行修改），分別作為family ID和sampleID。但是一般我們的樣本並不是那樣命名的，所以可以添加--double-id參數，將familyID和sampleID命名為相同，或者--const-fid,將familyID命名為0，表明-9

第二步 plink進行過濾，得到bed文件

plink --allow-extra-chr --noweb -file test --geno 0.05 --maf 0.05 --hwe 0.0001 --make-bed --out test1
# --noweb 不顯示網頁

第三步：尋找合適的K值

for i in {1..7};do admixture --cv test1.bed $i |tee log${i}.out;done
## 根據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

當k=4時，cv error 最小，選擇4
此外，
如果SNP數據集非常大，則可以隨機選擇SNP進行K值選擇分析，比如隨機選取20000個SNP進行分析，每個K值跑20次，確定最終的k值，然后分析

當利用plink轉的格式中，在運行上述命令是出現以下報錯

Invalid chromosome code!  Use integers

將*.bim中的第一列改為數值就可以了

第四步：多線程

admixture test1.bed 4 -j 5
## j， 線程數

第五步：作圖

ta1 = read.table("test1.4.Q")
head(ta1)
barplot(t(as.matrix(ta1)),col = rainbow(3),
        xlab = "Individual",
        ylab = "Ancestry",
        border = NA)

根據LD進行篩選SNP

首先 篩選合格的LD位點

## 對vcf進行
plink --vcf SNP.vcf.gz --indep-pairwise 100 50 0.2 -out Test-id-maf0.05-LD --allow-extra-chr

## 或者對bed進行操作也可以
plink --bfile  test1 --indep-pairwise 100 50 0.2

## --indep-pairwise 100 50 0.2； 100Kb窗口，50步長，0.2LD強度

然后進行提取

plink  --bfile test1 --extract plink.prune.in --make-bed --out test2

然后在進行和上述一樣的分析即可

歡迎掃碼交流

參考

--Admixture使用說明文檔cookbook