一,进行admixture分析的理解性问题
1.需不需需要对vcf进行去LD处理?
进化树、群体结构、PCA分析需要对LD进行质控吗? (qq.com)
注意看一下这篇文章的评论,我觉得也是,根据自己的数据量来确定,有文献去掉LD有文献说没有去掉:只要标记的覆盖度够强,哪个数据集应该都可以。感觉说白了只要你这个数据集里面包含所有的core snp就行。
2.admixture图怎么看怎么理解
群体结构分析structure原理 - 组学大讲堂问答社区 (omicsclass.com)
参看这篇文章足以。
3.admixture的操作流程
种群结构分析--ADMIXTURE和fastStructure - 知乎 (zhihu.com)
这篇文章写得很详细,这种做部分转载:
# step 1 vcf to plink format ,vcftools or plink
#vcftools --vcf sample_core_change_chr.vcf --plink --maf 0.05 --out plink_out
plink --vcf sample_core_change_chr.vcf -recode12 --out plink_out
# step 2
plink --noweb --file plink_out --hwe 0.0001 --make-bed --out QC
# step 3
for K in `seq 1 30`
do
admixture --cv QC.bed $K -j10 |tee log${K}.out
# -j 10 多线程 j和10之间不空格
done
# step 4 besk K 找cross-validation error值最低时所对应的K值为最合适的K值
grep -h CV log*.out |sort -nk4 -t ' ' > cross-validation_error.txt
#CV error (K=22): 0.06797
#CV error (K=20): 0.11223
#CV error (K=18): 0.13912
#CV error (K=21): 0.21628
#CV error (K=17): 0.25565
#CV error (K=11): 0.27448
#...
# step 5 画图
grep -h CV log*.out |sort -nk4 -t ' ' |tail -n +2 |perl -p -e 's/.*K=(.+?)\):.*/$1/g' |xargs -i rm QC.{}.Q QC.{}.P log{}.out
num=`grep -h CV log*.out |sort -nk4 -t ' ' |head -n 1 |perl -p -e 's/.*K=(.+?)\):.*/$1/g' `
paste <(cut -f 1 -d ' ' QC.fam ) <(cat QC.${num}.Q|tr ' ' '\t') > QC.${num}.Q.tab
二,pophelper包对 Admixture遗传结构分析可视化
Pophelper 2.3.1 • pophelper (royfrancis.com)
这里不多少说,参考官方的文件,这里直接给出我的代码流程
#先安装软件包,tidyverse,pophelper
library(conflicted)
library(dplyr)
library(tidyverse)
library(pophelper)
#切到使用的文件夹
setwd("G:..")
#颜色
clist <- list("shiny"=c("#8DD3C7","#1D72F5","#DF0101","#77CE61","#FF9326","#A945FF","#0089B2","#FDF060","#FFA6B2","#BFF217","#60D5FD","#CC1577","#F2B950","#7FB21D","#EC496F","#326397","#B26314","#027368","#A4A4A4","#610B5E"))
#给出样本名称,跟跑admixture的bed文件的fam文件一样的txt
label_list <- "G:.txt" #样本名
#举个例子类似
#fam文件
#1111 1111 0 0 0 0
#1122 1122 0 0 0 0
#txt文件(我把组别也带上了)
#1111 G1
#1122 G2
# 读取当前文件下所有的Q文件#将mul_slist中的每个元素的行名设置为label$V1的值。用了lapply函数将ID.txt应用到mul_slist的第一列。(这里大家只需要自己能够把当前文件下Q文件给弄好了就行)
mul_sfiles<- list.files(path=getwd(),pattern = '*.Q$')
mul_slists <- mul_sfiles[order(nchar(mul_sfiles),mul_sfiles)]
mul_slist <- readQ(files = mul_sfiles)
label <- read.delim(label_list,header=F,stringsAsFactors=F,sep = '\t')
if(length(unique(sapply(mul_slist,nrow)))==1)
mul_slist<-lapply(mul_slist,"rownames<-",label$V1)
#分群文件,其实刚刚的文件里面已经包含分群文件了
pop <- read.table('.txt', header = F, sep = '\t', stringsAsFactors = F,row.names = 1)
#1122 G2
#可视化,这里不能进行跨K,因为这个R包不具备跨K功能,但是也能出图
P2 <- plotQ((mul_slist),
imgoutput = "join",
imgtype = "pdf", #输出格式的变化
height = 1.2,width = 20,
showindlab = T,
useindlab = T,
sharedindlab = F,
showsp = T,
sppos = "left",
splabangle = 180,
showlegend = T,#图例
ordergrp = T,grplab = pop, #添加分组
sortind="all",
#sortind="Cluster1",
dpi = 800,
outputfilename = paste0("structure_barplot_pdf8"),units = "cm",
exportpath = getwd())
三,使用pong进行跨K画图
https://github.com/ramachandran-lab/pong
ADMIXTURE和STRUCTURE可视化:pong - 简书 (jianshu.com)
这里放上官网和参考文章,建议大家去官网里面下载参考数据
使用pong进行可视化的时候,需要注意一下几点
1.pong_filemap文件
第一列最好是用ABC字母表示,无任何意义
第二列是K值
第三列是Q文件的路径(我的是在当前文件下,所以直接把Q文件的名字放上去就行)
每列以Tab键分隔
A 4 c.4.Q
B 6 c.6.Q
C 5 c.5.Q
运行命令
pong -m file_map
运行命令的时候,你将会得到一个这个,复制这个在网页打开就可以得到图片,但是没有排序杂乱无章的
2.ind2pop.txt文件
这里只有一列,放上每一个样本的分组信息
注意的是,admixture的PQ文件是根据跑admixture的fam文件决定,这个地方只是给每一个样本带上一个分组信息
G6
G1
运行命令,这个时候你将会得到一个没有排序杂乱无章的但是有名称的图片
pong -m file_map -i indiv.txt
3.pop_order.txt文件
只有一列,刚刚补充分组信息,这个地方将会补充,把一个组的样本放在一起,呈现出图
G1
G2
G3
运行命令,将会出现下图,如果想要追求完美,可以去修图软件修改一下。
pong -m file_map -i indiv.txt -n pop_order.txt