群体遗传结构的分析并绘图
文章目录
群体遗传结构的分析并绘图
所属目录：紫菜
创建时间：2024/7/27
作者：星云<XingYun>
更新时间：2024/7/28
URL：https://blog.csdn.net/2301_78630677/article/details/140739916
前言
之前得到了经过SNP过滤后的vcf文件，本文就使用该文件进行后续的群体遗传结构分析
之前的博文：SNP过滤
本文的分析记录流程：
一、群体遗传结构分析
群体遗传学中测的很多个个体，得到了最终的SNP过滤后的vcf文件，需要将其分成群体，看哪几个物种聚在一起，一般使用软件admixture进行分析
推荐文章：
admixture 群体结构分析;
5. GWAS：群体结构——Admixture;
群体结构分析软件admixture安装及使用经验
1.1. 获得输入文件：plink的bed文件
#vcf转成ped和map
plink --vcf 149toTZC.2allele.filtered.missing0.1.maf0.01.recode.vcf --allow-extra-chr --recode --out 149toTZC_output
# ped、map转bed、fam和bim
plink --file 149toTZC_output --allow-extra-chr --make-bed --out 149toTZC_output
1.2. 运行脚本adm.sh
添加运行权限：
chmod +x adm.sh
#!/bin/bash
for K in 1 2 3 4 5 6 7 8 9; do admixture --cv 149toTZC_output.bed $K | tee log${K}.out; done
1.3. 最佳K值确定
每个K值对应的CV值（Cross-validation error）输出在标准输出文件中;
查看cv值，cv error最小的K值为最佳分群数
grep -h CV log*.out
【注：这里作者只考虑了1-9，这里就暂且定为k=6时有最小的cv值,】
1.4. 查看.Q文件(以149toTZC_output.3.Q为例)
列数与K值相对应，每行代表一个个体，数值代表在各个群体中的遗传组成
二、绘制堆积柱状图
2.1. 用*.Q 文件画图
setwd("")
data= read.table("149toTZC_output.6.Q");
head(data)
pdf("149toTZC_output_Q6.pdf",width = 9,height = 3)
barplot(t(as.matrix(data)),col = rainbow(6),
xlab = "Individual",
ylab = "Ancestry",
border = NA)
dev.off()
2.2. 使用pophelper包
pophelper的github链接：https://github.com/royfrancis/pophelper/
推荐阅读：
R包：Pophelper 作堆积柱状图
用 R 从 GitHub 上安装包（很大程度上解决网络问题）
（本人是采用从github上直接clone下来进行本地安装的方法）
注意：可能会涉及版本不兼容的问题，请根据报错提示更新相关包，可以考虑先删除包去重新安装。（当然如果没有报错自然更好）
library(devtools)
setwd("C:/Users/Desktop")
devtools::install_local("pophelper.zip")
library(pophelper)
sfiles <- list.files(path="Q_matrix", full.names=T)
slist <- readQ(files=sfiles,indlabfromfile=T)
head(slist[[1]])
custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
slist <- lapply(slist, function(x) {
rownames(x) <- custom_labels # 将自定义标签应用于每个Q矩阵
x
})
head(slist[[1]])
#调用plotQ()函数
plotQ(slist[2:5], imgoutput="join", showindlab=T,useindlab=T,
selgrp="loc", ordergrp=T, showlegend=F,
showtitle=F, showsubtitle=F, titlelab="The Great Structure",
subtitlelab="The GWAS population structure of my organism.",
height=0.7, indlabsize=0.6, # 减小字体大小
indlabheight=0.1, # 增加标签高度
indlabspacer=0, # 调整标签间距
barbordercolour="white", barbordersize=0,
outputfilename="mypop2", imgtype="pdf",
exportpath=getwd())
##（以下代码是有关版本不兼容问题所进行的，如果上面的安装pophelper包未报错，就不必要看了）
# 检查rlang版本
rlang_version <- packageVersion("rlang")
cat("Current rlang version:", format(rlang_version), "\n")
# 删除原来的rlang包
remove.packages("rlang")
#下载rlang包
install.packages("rlang")
主要代码解释
【整个脚本的主要目的是从一组Q矩阵文件中读取数据，然后生成一个群体结构图，展示个体在不同群体中的归属情况，并将结果保存为PDF文件】
#使用list.files()函数列出Q_matrix目录下的所有文件，并将结果存储在sfiles变量中。full.names=T表示返回完整的文件路径。
sfiles <- list.files(path="Q_matrix", full.names=T)
#读取Q矩阵文件：
#readQ()函数用于读取群体结构分析（如STRUCTURE分析）产生的Q矩阵文件，这些文件通常包含个体在不同群体中的归属概率。indlabfromfile=T表示从文件名中提取个体标签
slist <- readQ(files=sfiles,indlabfromfile=T)
#读取自定义标签
custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
#应用自定义标签到Q矩阵：
slist <- lapply(slist, function(x) {
rownames(x) <- custom_labels # 将自定义标签应用于每个Q矩阵
x
})
#调用plotQ()函数
plotQ(slist[2:5], imgoutput="join", showindlab=T,useindlab=T,
selgrp="loc", ordergrp=T, showlegend=F,
showtitle=F, showsubtitle=F, titlelab="The Great Structure",
subtitlelab="The GWAS population structure of my organism.",
height=0.7, indlabsize=0.6, # 减小字体大小
indlabheight=0.1, # 增加标签高度
indlabspacer=0, # 调整标签间距
barbordercolour="white", barbordersize=0,
outputfilename="mypop2", imgtype="pdf",
exportpath=getwd())
#plotQ()函数用于可视化群体结构。这里，我们选择了slist列表中的第2到第5个元素进行绘制。参数解释如下：
# imgoutput="join"：表示合并所有选定的群体结构图到一个图像中。
# showindlab=T：显示个体标签。
# useindlab=T：使用自定义标签。
# selgrp="loc"：选择按位置（loc）分组。
# ordergrp=T：按组排序。
# showlegend=T：显示图例。
# showtitle=T：显示标题。
# showsubtitle=T：显示副标题。
# titlelab和subtitlelab：设置标题和副标题的文本。
# height：设置图像的高度。
# indlabsize：设置个体标签的大小。
# indlabheight：设置个体标签的高度。
# indlabspacer：设置个体标签之间的间距。
# barbordercolour和barbordersize：设置条形图边框的颜色和大小。
# outputfilename：设置输出文件名。
# imgtype：设置图像类型，这里是PDF格式。
# exportpath：设置图像的导出路径，这里使用getwd()获取当前工作目录。
2.3. 使用TBtools(推荐，比较方便)
推荐文章：群体结构分析 | Pophelper 的“平替版”
最下面的 Sort Mode ，可以对个体进行排序，默认按照血统比例排序…
三、绘制地图分布
把K=6的时候画到地图上
推荐文章：
[R包分享] mapmixture包优雅绘制个性化地图
工具 | R包mapmixture绘制群体结构与地图分布
3.1. 计算原理（可跳过）
计算一个群体的遗传结构，是每一类加和取平均，比如CYS01到CYS11是一个群体，我就对这六列每列计算平均值，形成一行
使用R代码来进行自动计算：
Q6 <- slist[[5]]
# 提取行名中的字母前缀
prefixes <- gsub("[0-9].*", "", rownames(Q6))
# 将提取的前缀作为新的行名的一部分（这里实际上不需要修改行名，因为分组只需要前缀）
prefixes <- trimws(prefixes) # 移除前缀字符串两端的空白字符
# 分组数据
grouped_data <- split(Q6, prefixes)
# 示例：打印每个群体的前几行
for (group in names(grouped_data)) {
print(paste("Group:", group))
print(head(grouped_data[[group]], 3))
}
# 计算每个群体的平均值
group_means_list <- lapply(grouped_data, function(x) colMeans(x))
# 将平均值列表转换为数据框
group_means_df <- do.call(rbind, group_means_list)
# 设置行名为群体名称
rownames(group_means_df) <- names(grouped_data)
# 改变输出格式，避免使用科学计数法
options(scipen = 999)
# 打印结果数据框
print(group_means_df)
rownames(group_means_df)
得到这样的（形成了一个新的数据框来存储求和取平均后的值）：
【注意：这里只是展示一下计算处理的过程，方便理解；实际上我们可以用后面mapmixture包来进行群体结构的地图分布】
3.2. 使用mapmixture包绘制分布地图
使用mapmixture包绘制，要特定格式的输入数据，如下图所示：
需要一个存储了遗传结构信息的经admixture处理后的数据文件admixture1.csv，一个存储地理位置信息的数据文件coordinates.csv。
并且从中也可以看出它是会自动加和取平均
3.2.1. 先查看一下示例
install.packages("mapmixture")
library(tidyverse)
library(mapmixture)
library(terra)
library(gridExtra)
##(使用包中自带的文件去查看)
file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
admixture1<- read.csv(file)
file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
coordinates<- read.csv(file)
map1 <- mapmixture(admixture1, coordinates, crs = 3035)
map1
得到的示例图：
admixture1.csv的数据格式
coordinates.csv的数据格式：
3.2.2. 将之前的文件转为需要的数据格式
library(pophelper)
sfiles <- list.files(path="Q_matrix", full.names=T)
slist <- readQ(files=sfiles,indlabfromfile=T)
head(slist[[1]])
custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
slist <- lapply(slist, function(x) {
rownames(x) <- custom_labels # 将自定义标签应用于每个Q矩阵
x
})
slist[[5]]
#调用plotQ()函数
Q6 <- slist[[5]]
#得到admixture1.csv
# 将原始行名复制为新列
Q6 <- Q6 %>% mutate(Ind = rownames(.))
# 提取行名中的字母前缀
Q6 <- Q6 %>% mutate(Site = gsub("[0-9].*", "", Ind))
# 移除前缀字符串两端的空白字符
Q6$Site <- trimws(Q6$Site)
# 移除原始的行名
rownames(Q6) <- NULL
# 重新排序列，使Site和Ind成为第一和第二列
Q6 <- Q6 %>% select(Site, Ind, everything())
Q6
write.csv(Q6, file = "admixture1.csv", row.names = FALSE)
#得到coordinates.csv
gps <- read.csv("gps.csv",header = F) #读取原来的地理位置信息文件
gps
colnames(gps) <- c("Site", "Lat", "Lon") # 更改列名
gps
write.csv(gps, file = "coordinates.csv", row.names = FALSE)
补充：这是之前存放地理位置信息的gps.csv文件部分截图：
3.2.3. 正式绘图
earth <- terra::rast("D:/HYP_50M_SR_W/HYP_50M_SR_W.tif")
#如果需要，可以使用自定义图层（同时后面的函数中的参数basemap要指定自定义的图层）
admixture11 <- read.csv("admixture11.csv")
coordinates1 <- read.csv("coordinates1.csv")
cluster_cols <- c("#e69f00", "red", "#009e73", "#f0e442", "purple", "#0072b2")
map3 <- mapmixture(admixture11, coordinates1,
cluster_cols = cluster_cols,
crs = 4326,
boundary = c(xmin = 115, xmax = 125, ymin = 22, ymax = 32),
pie_size = 0.3,basemap = earth)+
theme(
axis.title = element_text(size = 10),
axis.text = element_text(size = 8))+
guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))
map3
默认地图
自定义图层地图
mapmixture 函数是用于绘制群体结构图的 R 包。以下是函数的详细用法和参数解释：
mapmixture(
admixture_df,
coords_df,
cluster_cols = NULL,
cluster_names = NULL,
boundary = NULL,
crs = 4326,
basemap = NULL,
pie_size = 1,
pie_border = 0.2,
pie_border_col = "black",
pie_opacity = 1,
land_colour = "#d9d9d9",
sea_colour = "#deebf7",
expand = FALSE,
arrow = TRUE,
arrow_size = 1,
arrow_position = "tl",
scalebar = TRUE,
scalebar_size = 1,
scalebar_position = "tl",
plot_title = "",
plot_title_size = 12,
axis_title_size = 10,
axis_text_size = 8
)
参数解释：
admixture_df: 包含 admixture 数据的数据框或 tibble。
coords_df: 包含坐标数据的数据框或 tibble。
cluster_cols: 长度与群数相同的颜色向量。如果为 NULL，则使用蓝色-绿色调色板。
cluster_names: 长度与群数相同的名称向量。如果为 NULL，则使用群列名称。
boundary: 定义地图边界的命名数值向量，例如 c(xmin = -15, xmax = 15, ymin = 30, ymax = 50)。如果为 NULL，则使用默认边界框。
crs: 坐标参考系统。默认为 WGS 84 - 世界地理坐标系 1984 (EPSG:4326)。
basemap: 用作底图的 SpatRaster 或 sf 对象。可以使用 terra 包的 rast() 函数从文件创建 SpatRaster 对象，或使用 sf 包的 st_read() 函数从文件创建 sf 对象。如果为 NULL，则使用世界国家边界。
pie_size: 饼图的大小（0 或更大）。
pie_border: 饼图边框的大小（0 或更大）。
pie_border_col: 饼图边框的颜色。
pie_opacity: 饼图的透明度（0 到 1）。
land_colour: 陆地颜色的字符串。
sea_colour: 海水颜色的字符串。
expand: 是否扩展坐标轴（TRUE 或 FALSE）。
arrow: 是否显示箭头（TRUE 或 FALSE）。使用 ggspatial 包的 annotation_north_arrow() 函数添加。
arrow_size: 箭头的大小（0 或更大）。
arrow_position: 箭头的position（"tl"、"tr"、"bl" 或 "br"）。
scalebar: 是否显示比例尺（TRUE 或 FALSE）。使用 ggspatial 包的 annotation_scale() 函数添加。
scalebar_size: 比例尺的大小（0 或更大）。
scalebar_position: 比例尺的position（"tl"、"tr"、"bl" 或 "br"）。
plot_title: 图像的主标题。
plot_title_size: 图像主标题的大小（0 或更大）。
axis_title_size: 坐标轴标题的大小（0 或更大）。
axis_text_size: 坐标轴文本的大小（0 或更大）。
3.3. 亦可绘制堆积柱状图（补充）
在使用这个mapmixture包时发现它也可以绘制之前的堆积柱状图
earth <- terra::rast("D:/HYP_50M_SR_W/HYP_50M_SR_W.tif")
admixture11 <- read.csv("admixture11.csv")
coordinates1 <- read.csv("coordinates1.csv")
cluster_cols <- c("#e69f00", "red", "#009e73", "#f0e442", "purple", "#0072b2")
map3 <- mapmixture(admixture11, coordinates1,
cluster_cols = cluster_cols,
crs = 4326,
boundary = c(xmin = 115, xmax = 125, ymin = 22, ymax = 32),
pie_size = 0.3,basemap = earth)+
theme(
axis.title = element_text(size = 10),
axis.text = element_text(size = 8),legend.position = "top",plot.margin = margin(l = 10, r = 10))+
guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))
map3
#绘制图1
structure_barplot <- structure_plot(admixture11,
type = "structure",
cluster_cols = cluster_cols,
site_dividers = TRUE,
divider_width = 0.4,
site_order = c(
"GYS","LJS","ND","NMW","YJC","DDC","DH","DLS","DMY","DXD","JJS","KNW","NJDHJ","PY","WYS",
"WZ","YSD","CX","HND","NJD","SDC","XDC"
),
labels = "site",
flip_axis = FALSE,
site_ticks_size = -0.05,
site_labels_y = -0.35,
site_labels_size = 2.2)+
theme(
axis.title.y = element_text(size = 8, hjust = 1),
axis.text.y = element_text(size = 5),
)
grid.arrange(map3, structure_barplot, nrow = 2, heights = c(4,1))
## 绘制图2
facet_barplot <- structure_plot(admixture11,
type = "facet",
cluster_cols = cluster_cols,
facet_col = 2,
ylabel = "Admixture proportions",
)+
theme(
axis.title.y = element_text(size = 10),
axis.text.y = element_text(size = 5),
strip.text = element_text(size = 6, vjust = 1, margin = margin(t=1.5, r=0, b=1.5, l=0)),
)
grid.arrange(map3, facet_barplot, ncol = 2, widths = c(3,2))
图1
图2
总结
本文记录了一下群体遗传结构的分析并绘图过程中所做的笔记。
主要使用了admixture软件进行了群体遗传结构分析；
使用了pophelper包和TBtools软件进行了遗传结构堆积柱状图的绘制；
使用了mapmixture包进行了遗传结构的地图分布，同时其也可用于前面的堆积柱状图的绘制。
差不多就这些（遇到的主要难题还是前面pophelper包的安装，因为版本兼容等问题，花费了些许时间）。。。
–2024/7/28
标签：文件,遗传,--,mapmixture,使用,绘图,群体,size From： https://blog.csdn.net/2301_78630677/article/details/140739916