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

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

---

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

1. 所属目录：紫菜
3. 创建时间：2024/7/27
5. 作者：星云<XingYun>
7. 更新时间：2024/7/28
9. 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文件

1. #vcf转成ped和map
2. plink --vcf 149toTZC.2allele.filtered.missing0.1.maf0.01.recode.vcf --allow-extra-chr --recode --out 149toTZC_output
4. # ped、map转bed、fam和bim
5. plink --file 149toTZC_output --allow-extra-chr --make-bed --out 149toTZC_output

复制代码

1.2. 运行脚本adm.sh

添加运行权限：
chmod +x adm.sh

1. #!/bin/bash
2. 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值为最佳分群数

1. grep -h CV log*.out

复制代码

【注：这里作者只考虑了1-9，这里就临时定为k=6时有最小的cv值,】
1.4. 查看.Q文件(以149toTZC_output.3.Q为例)

列数与K值相对应，每行代表一个个体，数值代表在各个群体中的遗传组成

二、绘制堆积柱状图

2.1. 用*.Q 文件画图

1. setwd("")
2. data= read.table("149toTZC_output.6.Q");
3. head(data)
4. pdf("149toTZC_output_Q6.pdf",width = 9,height = 3)
5. barplot(t(as.matrix(data)),col = rainbow(6),
6.         xlab = "Individual",
7.         ylab = "Ancestry",
8.         border = NA)
9. dev.off()

复制代码

2.2. 使用pophelper包

pophelper的github链接：https://github.com/royfrancis/pophelper/
保举阅读：
R包：Pophelper 作堆积柱状图
用 R 从 GitHub 上安装包（很大程度上办理网络题目）
（本人是采用从github上直接clone下来举行本地安装的方法）
注意：大概会涉及版本不兼容的题目，请根据报错提示更新相关包，可以考虑先删除包去重新安装。（当然如果没有报错自然更好）

1. library(devtools)
2. setwd("C:/Users/Desktop")
3. devtools::install_local("pophelper.zip")
4. library(pophelper)
5. sfiles <- list.files(path="Q_matrix", full.names=T)
6. slist <- readQ(files=sfiles,indlabfromfile=T)
7. head(slist[[1]])
8. custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
9. slist <- lapply(slist, function(x) {
10.   rownames(x) <- custom_labels  # 将自定义标签应用于每个Q矩阵
11.   x
12. })
13. head(slist[[1]])
14. #调用plotQ()函数
15. plotQ(slist[2:5], imgoutput="join", showindlab=T,useindlab=T,
16.       selgrp="loc", ordergrp=T, showlegend=F,
17.       showtitle=F, showsubtitle=F, titlelab="The Great Structure",
18.       subtitlelab="The GWAS population structure of my organism.",
19.       height=0.7, indlabsize=0.6,  # 减小字体大小
20.       indlabheight=0.1,           # 增加标签高度
21.       indlabspacer=0,             # 调整标签间距
22.       barbordercolour="white", barbordersize=0,
23.       outputfilename="mypop2", imgtype="pdf",
24.       exportpath=getwd())
29. ##（以下代码是有关版本不兼容问题所进行的，如果上面的安装pophelper包未报错，就不必要看了）
30. # 检查rlang版本
31. rlang_version <- packageVersion("rlang")
32. cat("Current rlang version:", format(rlang_version), "\n")
33. # 删除原来的rlang包
34. remove.packages("rlang")
35. #下载rlang包
36. install.packages("rlang")

复制代码

主要代码解释
【整个脚本的主要目的是从一组Q矩阵文件中读取数据，然后生成一个群体结构图，展示个体在不同群体中的归属环境，并将结果保存为PDF文件】

1. #使用list.files()函数列出Q_matrix目录下的所有文件，并将结果存储在sfiles变量中。full.names=T表示返回完整的文件路径。
2. sfiles <- list.files(path="Q_matrix", full.names=T)
4. #读取Q矩阵文件：
5. #readQ()函数用于读取群体结构分析（如STRUCTURE分析）产生的Q矩阵文件，这些文件通常包含个体在不同群体中的归属概率。indlabfromfile=T表示从文件名中提取个体标签
6. slist <- readQ(files=sfiles,indlabfromfile=T)
8. #读取自定义标签
9. custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
11. #应用自定义标签到Q矩阵：
12. slist <- lapply(slist, function(x) {
13.   rownames(x) <- custom_labels  # 将自定义标签应用于每个Q矩阵
14.   x
15. })
18. #调用plotQ()函数
19. plotQ(slist[2:5], imgoutput="join", showindlab=T,useindlab=T,
20.       selgrp="loc", ordergrp=T, showlegend=F,
21.       showtitle=F, showsubtitle=F, titlelab="The Great Structure",
22.       subtitlelab="The GWAS population structure of my organism.",
23.       height=0.7, indlabsize=0.6,  # 减小字体大小
24.       indlabheight=0.1,           # 增加标签高度
25.       indlabspacer=0,             # 调整标签间距
26.       barbordercolour="white", barbordersize=0,
27.       outputfilename="mypop2", imgtype="pdf",
28.       exportpath=getwd())
29. #plotQ()函数用于可视化群体结构。这里，我们选择了slist列表中的第2到第5个元素进行绘制。参数解释如下：
30. # imgoutput="join"：表示合并所有选定的群体结构图到一个图像中。
31. # showindlab=T：显示个体标签。
32. # useindlab=T：使用自定义标签。
33. # selgrp="loc"：选择按位置（loc）分组。
34. # ordergrp=T：按组排序。
35. # showlegend=T：显示图例。
36. # showtitle=T：显示标题。
37. # showsubtitle=T：显示副标题。
38. # titlelab和subtitlelab：设置标题和副标题的文本。
39. # height：设置图像的高度。
40. # indlabsize：设置个体标签的大小。
41. # indlabheight：设置个体标签的高度。
42. # indlabspacer：设置个体标签之间的间距。
43. # barbordercolour和barbordersize：设置条形图边框的颜色和大小。
44. # outputfilename：设置输出文件名。
45. # imgtype：设置图像类型，这里是PDF格式。
46. # exportpath：设置图像的导出路径，这里使用getwd()获取当前工作目录。

复制代码

2.3. 使用TBtools(保举，比力方便)

保举文章：群体结构分析 | Pophelper 的“平替版”
最下面的 Sort Mode ，可以对个体举行排序，默认按照血统比例排序…

三、绘制地图分布

   把K=6的时间画到地图上
  保举文章：
[R包分享] mapmixture包优雅绘制个性化地图
工具 | R包mapmixture绘制群体结构与地图分布
批量转换地点为经纬度网站
3.1. 计算原理（可跳过）

   计算一个群体的遗传结构，是每一类加和取平均，比如CYS01到CYS11是一个群体，我就对这六列每列计算平均值，形成一行
  

   使用R代码来举行主动计算：

1. Q6 <- slist[[5]]
2. # 提取行名中的字母前缀
3. prefixes <- gsub("[0-9].*", "", rownames(Q6))
4. # 将提取的前缀作为新的行名的一部分（这里实际上不需要修改行名，因为分组只需要前缀）
5. prefixes <- trimws(prefixes) # 移除前缀字符串两端的空白字符
6. # 分组数据
7. grouped_data <- split(Q6, prefixes)
8. # 示例：打印每个群体的前几行
9. for (group in names(grouped_data)) {
10.   print(paste("Group:", group))
11.   print(head(grouped_data[[group]], 3))
12. }
13. # 计算每个群体的平均值
14. group_means_list <- lapply(grouped_data, function(x) colMeans(x))
15. # 将平均值列表转换为数据框
16. group_means_df <- do.call(rbind, group_means_list)
17. # 设置行名为群体名称
18. rownames(group_means_df) <- names(grouped_data)
19. # 改变输出格式，避免使用科学计数法
20. options(scipen = 999)
21. # 打印结果数据框
22. print(group_means_df)
23. rownames(group_means_df)

复制代码

  得到如许的（形成了一个新的数据框来存储求和取平均后的值）：
  

【注意：这里只是展示一下计算处理的过程，方便明确；现实上我们可以用后面mapmixture包来举行群体结构的地图分布】
3.2. 使用mapmixture包绘制分布地图

   使用mapmixture包绘制，要特定格式的输入数据，如下图所示：
需要一个存储了遗传结构信息的经admixture处理后的数据文件admixture1.csv，一个存储地理位置信息的数据文件coordinates.csv。
并且从中也可以看出它是会主动加和取平均
  

3.2.1. 先查看一下示例

1. install.packages("mapmixture")
2. library(tidyverse)
3. library(mapmixture)
4. library(terra)
5. library(gridExtra)
7. ##(使用包中自带的文件去查看)
8. file <- system.file("extdata", "admixture1.csv", package = "mapmixture")
9. admixture1<- read.csv(file)
10. file <- system.file("extdata", "coordinates.csv", package = "mapmixture")
11. coordinates<- read.csv(file)
12. map1 <- mapmixture(admixture1, coordinates, crs = 3035)
13. map1

复制代码

  得到的示例图：

    admixture1.csv的数据格式

    coordinates.csv的数据格式：

  3.2.2. 将之前的文件转为需要的数据格式

1. library(pophelper)
2. sfiles <- list.files(path="Q_matrix", full.names=T)
3. slist <- readQ(files=sfiles,indlabfromfile=T)
4. head(slist[[1]])
5. custom_labels <- read.table("custom_labels.txt", header = FALSE, stringsAsFactors = FALSE)$V1
6. slist <- lapply(slist, function(x) {
7.   rownames(x) <- custom_labels  # 将自定义标签应用于每个Q矩阵
8.   x
9. })
10. slist[[5]]
11. #调用plotQ()函数
12. Q6 <- slist[[5]]
16. #得到admixture1.csv
17. # 将原始行名复制为新列
18. Q6 <- Q6 %>% mutate(Ind = rownames(.))
19. # 提取行名中的字母前缀
20. Q6 <- Q6 %>% mutate(Site = gsub("[0-9].*", "", Ind))
21. # 移除前缀字符串两端的空白字符
22. Q6$Site <- trimws(Q6$Site)
23. # 移除原始的行名
24. rownames(Q6) <- NULL
25. # 重新排序列，使Site和Ind成为第一和第二列
26. Q6 <- Q6 %>% select(Site, Ind, everything())
27. Q6
28. write.csv(Q6, file = "admixture1.csv", row.names = FALSE)
31. #得到coordinates.csv
32. gps <- read.csv("gps.csv",header = F)  #读取原来的地理位置信息文件
33. gps
34. colnames(gps) <- c("Site", "Lat", "Lon")  # 更改列名
35. gps
36. write.csv(gps, file = "coordinates.csv", row.names = FALSE)

复制代码

  增补：这是之前存放地理位置信息的gps.csv文件部分截图：

  3.2.3. 正式绘图

1. earth <- terra::rast("D:/HYP_50M_SR_W/HYP_50M_SR_W.tif")  
2. #如果需要，可以使用自定义图层（同时后面的函数中的参数basemap要指定自定义的图层）
4. admixture11 <- read.csv("admixture11.csv")
5. coordinates1 <- read.csv("coordinates1.csv")
6. cluster_cols <- c("#e69f00", "red", "#009e73", "#f0e442", "purple", "#0072b2")
7. map3 <- mapmixture(admixture11, coordinates1,
8.                    cluster_cols = cluster_cols,
9.                    crs = 4326,
10.                    boundary = c(xmin = 115, xmax = 125, ymin = 22, ymax = 32),
11.                    pie_size = 0.3,basemap = earth)+
12.   theme(
13.     axis.title = element_text(size = 10),
14.     axis.text = element_text(size = 8))+
15.   guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))
16. map3

复制代码

  默认地图

    自定义图层地图

    mapmixture 函数是用于绘制群体结构图的 R 包。以下是函数的具体用法和参数解释：
  

1. mapmixture(
2.   admixture_df,
3.   coords_df,
4.   cluster_cols = NULL,
5.   cluster_names = NULL,
6.   boundary = NULL,
7.   crs = 4326,
8.   basemap = NULL,
9.   pie_size = 1,
10.   pie_border = 0.2,
11.   pie_border_col = "black",
12.   pie_opacity = 1,
13.   land_colour = "#d9d9d9",
14.   sea_colour = "#deebf7",
15.   expand = FALSE,
16.   arrow = TRUE,
17.   arrow_size = 1,
18.   arrow_position = "tl",
19.   scalebar = TRUE,
20.   scalebar_size = 1,
21.   scalebar_position = "tl",
22.   plot_title = "",
23.   plot_title_size = 12,
24.   axis_title_size = 10,
25.   axis_text_size = 8
26. )
29. 参数解释：
30. admixture_df: 包含 admixture 数据的数据框或 tibble。
31. coords_df: 包含坐标数据的数据框或 tibble。
32. cluster_cols: 长度与群数相同的颜色向量。如果为 NULL，则使用蓝色-绿色调色板。
33. cluster_names: 长度与群数相同的名称向量。如果为 NULL，则使用群列名称。
34. boundary: 定义地图边界的命名数值向量，例如 c(xmin = -15, xmax = 15, ymin = 30, ymax = 50)。如果为 NULL，则使用默认边界框。
35. crs: 坐标参考系统。默认为 WGS 84 - 世界地理坐标系 1984 (EPSG:4326)。
36. basemap: 用作底图的 SpatRaster 或 sf 对象。可以使用 terra 包的 rast() 函数从文件创建 SpatRaster 对象，或使用 sf 包的 st_read() 函数从文件创建 sf 对象。如果为 NULL，则使用世界国家边界。
37. pie_size: 饼图的大小（0 或更大）。
38. pie_border: 饼图边框的大小（0 或更大）。
39. pie_border_col: 饼图边框的颜色。
40. pie_opacity: 饼图的透明度（0 到 1）。
41. land_colour: 陆地颜色的字符串。
42. sea_colour: 海水颜色的字符串。
43. expand: 是否扩展坐标轴（TRUE 或 FALSE）。
44. arrow: 是否显示箭头（TRUE 或 FALSE）。使用 ggspatial 包的 annotation_north_arrow() 函数添加。
45. arrow_size: 箭头的大小（0 或更大）。
46. arrow_position: 箭头的position（"tl"、"tr"、"bl" 或 "br"）。
47. scalebar: 是否显示比例尺（TRUE 或 FALSE）。使用 ggspatial 包的 annotation_scale() 函数添加。
48. scalebar_size: 比例尺的大小（0 或更大）。
49. scalebar_position: 比例尺的position（"tl"、"tr"、"bl" 或 "br"）。
50. plot_title: 图像的主标题。
51. plot_title_size: 图像主标题的大小（0 或更大）。
52. axis_title_size: 坐标轴标题的大小（0 或更大）。
53. axis_text_size: 坐标轴文本的大小（0 或更大）。

复制代码

3.3. 亦可绘制堆积柱状图（增补）

   在使用这个mapmixture包时发现它也可以绘制之前的堆积柱状图

1. earth <- terra::rast("D:/HYP_50M_SR_W/HYP_50M_SR_W.tif")
2. admixture11 <- read.csv("admixture11.csv")
3. coordinates1 <- read.csv("coordinates1.csv")
4. cluster_cols <- c("#e69f00", "red", "#009e73", "#f0e442", "purple", "#0072b2")
5. map3 <- mapmixture(admixture11, coordinates1,
6.                    cluster_cols = cluster_cols,
7.                    crs = 4326,
8.                    boundary = c(xmin = 115, xmax = 125, ymin = 22, ymax = 32),
9.                    pie_size = 0.3,basemap = earth)+
10.   theme(
11.     axis.title = element_text(size = 10),
12.     axis.text = element_text(size = 8),legend.position = "top",plot.margin = margin(l = 10, r = 10))+
13.   guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1)))
14. map3
17. #绘制图1
18. structure_barplot <- structure_plot(admixture11,
19.                                     type = "structure",
20.                                     cluster_cols = cluster_cols,
21.                                     site_dividers = TRUE,
22.                                     divider_width = 0.4,
23.                                     site_order = c(
24.                                       "GYS","LJS","ND","NMW","YJC","DDC","DH","DLS","DMY","DXD","JJS","KNW","NJDHJ","PY","WYS",
25.                                       "WZ","YSD","CX","HND","NJD","SDC","XDC"
26.                                     ),
27.                                     labels = "site",
28.                                     flip_axis = FALSE,
29.                                     site_ticks_size = -0.05,
30.                                     site_labels_y = -0.35,
31.                                     site_labels_size = 2.2)+
32.   theme(
33.     axis.title.y = element_text(size = 8, hjust = 1),
34.     axis.text.y = element_text(size = 5),
35.   )
36. grid.arrange(map3, structure_barplot, nrow = 2, heights = c(4,1))
40. ## 绘制图2
41. facet_barplot <- structure_plot(admixture11,
42.                                 type = "facet",
43.                                 cluster_cols = cluster_cols,
44.                                 facet_col = 2,
45.                                 ylabel = "Admixture proportions",
46. )+
47.   theme(
48.     axis.title.y = element_text(size = 10),
49.     axis.text.y = element_text(size = 5),
50.     strip.text = element_text(size = 6, vjust = 1, margin = margin(t=1.5, r=0, b=1.5, l=0)),
51.   )
53. grid.arrange(map3, facet_barplot, ncol = 2, widths = c(3,2))

复制代码

  图1

    图2

  

---

总结

本文记录了一下群体遗传结构的分析并绘图过程中所做的笔记。
主要使用了admixture软件举行了群体遗传结构分析；
使用了pophelper包和TBtools软件举行了遗传结构堆积柱状图的绘制；
使用了mapmixture包举行了遗传结构的地图分布，同时其也可用于前面的堆积柱状图的绘制。
差不多就这些（遇到的主要困难还是前面pophelper包的安装，因为版本兼容等题目，花费了些许时间）。。。
–2024/7/28

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。