admixture全流程+跨K颜色处理

BioXplorer

已于 2024-05-28 20:23:20 修改

阅读量474

点赞数 4

文章标签： r语言学习

于 2023-10-10 13:43:23 首次发布

本文链接：https://blog.csdn.net/2302_79242191/article/details/133744911

版权

一，进行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