群体遗传结构的分析并绘图

星石传说

已于 2024-07-28 14:56:34 修改

阅读量572

点赞数 27

分类专栏：实验室——生信分析文章标签：生信 R linux

于 2024-07-28 14:55:49 首次发布

本文链接：https://blog.csdn.net/2301_78630677/article/details/140739916

版权

实验室——生信分析专栏收录该内容

21 篇文章 19 订阅

订阅专栏

群体遗传结构的分析并绘图

文章目录

群体遗传结构的分析并绘图
总结

群体遗传结构的分析并绘图

所属目录：紫菜

创建时间：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/

（本人是采用从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的时候画到地图上

批量转换地址为经纬度网站

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

星石传说

关注

27
点赞
踩
8

收藏

觉得还不错? 一键收藏
打赏
0
评论
群体遗传结构的分析并绘图

之前得到了经过SNP过滤后的vcf文件，本文就使用该文件进行后续的群体遗传结构分析之前的博文：SNP过滤本文记录了一下群体遗传结构的分析并绘图过程中所做的笔记。主要使用了admixture软件进行了群体遗传结构分析；使用了pophelper包和TBtools软件进行了遗传结构堆积柱状图的绘制；使用了mapmixture包进行了遗传结构的地图分布，同时其也可用于前面的堆积柱状图的绘制。差不多就这些（遇到的主要难题还是前面pophelper包的安装，因为版本兼容等问题，花费了些许时间）。。。
复制链接

扫一扫