admixture 群体结构分析

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/

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,242评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,769评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,484评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,133评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,007评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,080评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,496评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,190评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,464评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,549评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,330评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,205评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,567评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,889评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,160评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,475评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,650评论 2 335

推荐阅读更多精彩内容

  • tructure是与PCA、进化树相似的方法,就是利用分子标记的基因型信息对一组样本进行分类,分子标记可以是SNP...
    斩毛毛阅读 14,274评论 7 23
  • 群体遗传学中测的很多个个体,得到了最终的SNP vcf文件,需要将其分成群体,看那几个物种聚在一起,一般使用的软件...
    zhusitao阅读 13,277评论 10 14
  • 1. 软件下载及安装admixture: 使用conda进行软件安装 conda installadmixture...
    雨林课堂阅读 4,185评论 0 5
  • 目的:对群体结构和亲缘关系进行评估以确定使用的统计模型和获得相应的矩阵 评估内容(遗传上差异过大应剔除,相似性高的...
    滴滴_54f1阅读 7,402评论 0 8
  • 群体结构是指材料的亚群分化情况,会导致标记间的非连锁关联,进而导致关联分析结果出现假阳性。 地理隔离、人工选择、移...
    Wei_Sun阅读 15,118评论 3 27