admixture 群体结构分析

2023-04-28 本文已影响0人 iBioinformatics

http://dalexander.github.io/admixture/admixture-manual.pdf

群体遗传学中测的很多个个体，得到了最终的SNP vcf文件，需要将其分成群体，看哪几个物种聚在一起，一般使用的软件就是STRUCTURE，但是STREUTURE运行速度极慢，后面frappe软件提升了速度，但是也不是很快；admixture凭借其运算速度，成为了主流的分析软件。

一、软件安装

conda install -c bioconda admixture

wget http://dalexander.github.io/admixture/binaries/admixture_linux-1.3.0.tar.gz
tar -zvxf admixture_linux-1.3.0.tar.gz
cd admixture_linux-1.3.0
./admixture --version


Usage: admixture <input file> <K>
See --help or manual for more advanced usage.

二、简单使用

对于每一个K值，可以跑20次
admixture 软件一共分为5步：

1. vcf文件转化成plink格式

ADMIXTURE输入文件为二进制的PLINK文件（.bed），同时将对应的bim和fam文件放在同一个文件夹下，对于我们得到到非bed格式文件，使用Plink1.9软件进行转化，即可得到*.bed, *.bim, *.fam文件

将vcf文件转化成ped文件

/soft/vcftools-master/bin/vcftools \
   -vcf total.final.snp.vcf  \
   --plink  \
   --out xj

利用 plink 对 ped 文件进行过滤生成 bed 文件

/soft/plink/plink  \
   --allow-extra-chr   \
   --noweb  \
   --file xj  \
   --geno 0.05  \
   --maf 0.05  \
   --hwe 0.0001  \
   --make-bed  \
   --out QC

--const-fid 添加群体信息familyID sampleID）(将family设置为0)；
--allow-extra-chr 允许非标准染色体编号
--noweb 不显示网页

因为plink本身是针对人类进行开发的，所以在运行时，加上--allow-extr-chr。此外对于vcf中的sampleID（familyID_sampleID）, plink 默认为下划线分隔（也可以通过参数--id-delim进行修改），分别作为family ID和sampleID。但是一般我们的样本并不是那样命名的，所以可以添加--double-id参数，将familyID和sampleID命名为相同，或者--const-fid,将familyID命名为0，表明-9

2. 寻找合适的K值

利用 admixture 计算

for K in  {1..7};do /soft/admixture --cv QC.bed $K | tee log${K}.out;done

besk 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中的第一列改为数值就可以了

3. 多线程

admixture test1.bed 4 -j 5

j：线程数

4. R画图

data= read.table("./result/QC.4.Q");
head(data)
pdf("./result/Q4.pdf")
barplot(t(as.matrix(data)),col = rainbow(3),
        xlab = "Individual",
        ylab = "Ancestry",
        border = NA,space = 0)
dev.off()

使用pophelper R包进行群体结构绘制

mkdir Q_mattrix
cd Q_matrix
ln -s ../*Q .

library(pophelper)
sfiles <- list.files(path="Q_matrix", full.names=T)
slist <- readQ(files=sfiles,indlabfromfile=T)
plotQ(slist[2:5],imgoutput="join",showindlab=T,
      selgrp="loc",ordergrp=T,showlegend=T,
      showtitle=T,showsubtitle=T,titlelab="The Great Structure",
      subtitlelab="The GWAS population structure of my organism.",
      height=1.6,indlabsize=2.3,indlabheight=0.08,indlabspacer=-1,
      barbordercolour="white",barbordersize=0,outputfilename="mypop",imgtype="pdf",
      exportpath=getwd())

并且有讨论是否需要进行根据连锁情况进行标记的过滤以及如何使用R进行结果图形的绘制
https://vcru.wisc.edu/simonlab/bioinformatics/programs/admixture/admixture-manual.pdf

# removal each SNP that has an R2 value of greater than 0.1 with any 
# other SNP within a 50-SNP sliding window (advanced by 10 SNPs each time). 
plink \
   --bfile rawData  \
   --indep-pairwise 50 10 0.1 
plink  \
   --bfile rawData  \
   --extract plink.prune.in  \
   --make-bed  \
   --out prunedDatatbl=read.table("hapmap3.3.Q") 

# 画图 
barplot(t(as.matrix(tbl)), col=rainbow(3),xlab="Individual #", ylab="Ancestry", border=NA)

参考：
https://www.jianshu.com/p/5105f0edc152
https://mp.weixin.qq.com/s/nvgKkjzUyQ0Oy4mHoeZTOg
https://mp.weixin.qq.com/s/EV4KjzFqe6AYqC9lQ7YSmw
https://www.bilibili.com/read/cv6093335/