跟着NC学数据分析:R语言用分子距离/环境距离/地理距离做mantel检验

论文

Genomic insights into local adaptation and future climate-induced vulnerability of a keystone forest tree in East Asia
https://www.nature.com/articles/s41467-022-34206-8
论文中基本提供了每个图的原始数据,也提供了大部分代码,原始数据在电脑上的存储名41467_2022_34206_MOESM8_ESM.zip
论文中提供的代码链接
https://github.com/jingwanglab/Populus_genomic_prediction_climate_vulnerability
但是在这部分代码里没有找到 做mantel的代码
mantel检验主要是用来做两个距离矩阵之间的相关性
论文中关于这部分分析的方法部分
To investigate and compare the role of geography and environment in shaping spatial genetic variation of adaptive (the 1779 adaptive variants identified by both LFMM and RDA) and neutral (the 535,191 LD-pruned SNPs as used for population structure analyses) variants, Mantel and partial Mantel tests were separately used to test for associations between _F_ST(_F_ST/1−_F_ST) and geographic (IBD) and environmental (IBE) distance (after accounting for the geographic distance) with significance determined using 999 permutations in the R package vegan87, where environmental distance was represented by Euclidean distance of all scaled environmental variables.
分子距离这里用的是FST相关,但是为什么用 _F_ST(_F_ST/1−_F_ST) 这个值暂时没有想明白。地理距离的算法在 Genetic subdivision and candidate genes under selection in North American grey wolves 这篇论文里有提到
https://datadryad.org/stash/dataset/doi:10.5061/dryad.c9b25
地理距离是用Genalex 这个软件算的
这篇论文里提到的分子距离是用plink算的
NC的论文里已经提供了算好的距离矩阵,今天的推文里就直接用距离矩阵做mantel检验,复现论文中的Fig2C
读取数据

library(tidyverse)

read.csv(“data/20221211/Sourcedata/Fig2c&d/Fst_adaption.csv”,
row.names = 1) %>%
as.matrix() %>%
Matrix::tril() %>%
as.dist() -> fst.adaption.dist

read.csv(“data/20221211/Sourcedata/Fig2c&d/Fst_neutral.csv”,
row.names = 1) %>%
as.matrix() %>%
Matrix::tril() %>%
as.dist() -> fst.neutral.dist

read.csv(“data/20221211/Sourcedata/Fig2c&d/geo_dist.csv”,
row.names = 1) %>%
as.matrix() %>%
Matrix::tril() %>%
as.dist() -> geo.dist
这里提供的数据是一个对称矩阵,读取进来以后是一个数据框,最后转换成了一个下三角矩阵用于 mantel函数的输入
做mantel检验的函数

vegan::mantel(fst.adaption.dist,geo.dist,permutations = 999)
vegan::mantel(fst.neutral.dist,geo.dist,permutations = 999)

image.png
输出结果和论文中的一致
做partial mantel检验

vegan::mantel.partial(fst.adaption.dist,env.dist,geo.dist,permutations = 999)
vegan::mantel.partial(fst.neutral.dist,env.dist,geo.dist,permutations = 999)
这个结果和论文中的也一致
作图代码

dat<-read_csv(“data/20221211/Sourcedata/Fig2c&d/plot_data.csv”)

tt1<-ttheme_minimal(core=list(fg_params=list(hjust=0,
parse=TRUE,
x=0,
fontsize=15,
col=c(“#db6786″,”#db6786”,
“#609ac6″,”#609ac6”)),
bg_params=list(fill=”#e1eff5″)))

table2c<-tibble(x=c(“italic(Mantel)minuteitalic(s)~italic(r)==’0.450′”,
“italic(Mantel)minuteitalic(s)~italic(p)==’0.001′”,
“italic(Mantel)minuteitalic(s)~italic(r)==’0.451′”,
“italic(Mantel)minuteitalic(s)~italic(p)==’0.001′”))

p2c<-dat %>%
select(geo_dit,fst1779,fst53w) %>%
pivot_longer(!geo_dit) %>%
ggplot(aes(x=geo_dit/100000,y=value))+
geom_point(aes(color=name,
fill=name),
shape=21,
size=5)+
scale_color_manual(values = c(“#609ac6″,”#db6786”))+
scale_fill_manual(values = c(“#d6dde2″,”#fbf3d5”))+
geom_smooth(aes(color=name),
method = “lm”,
formula = ‘y~x’)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank(),
legend.position = “none”)+
labs(x=”Geographic distance (100 km)”,
y=expression(italic(F)[ST]/(1-italic(F)[ST])))+
annotation_custom(tableGrob(table2c,
rows = NULL,
cols = NULL,
theme = tt1),
xmin=0.5, xmax=5, ymin=0.75, ymax=1)

p2c

table2d<-tibble(x=c(“italic(partial~Mantel)minuteitalic(s)~italic(r)==’0.676′”,
“italic(partial~Mantel)minuteitalic(s)~italic(p)==’0.002′”,
“italic(partial~Mantel)minuteitalic(s)~italic(r)==’0.180′”,
“italic(partial~Mantel)minuteitalic(s)~italic(p)==’0.117′”))

p2d<-dat %>%
select(env_dist,fst1779,fst53w) %>%
pivot_longer(!env_dist) %>%
ggplot(aes(x=env_dist,y=value))+
geom_point(aes(color=name,
fill=name),
shape=21,
size=5)+
scale_color_manual(values = c(“#609ac6″,”#db6786”))+
scale_fill_manual(values = c(“#d6dde2″,”#fbf3d5”))+
geom_smooth(aes(color=name),
method = “lm”,
formula = ‘y~x’)+
theme_bw(base_size = 15)+
theme(panel.grid = element_blank(),
legend.position = “none”)+
labs(x=”Environmental distance”,
y=expression(italic(F)[ST]/(1-italic(F)[ST])))+
annotation_custom(tableGrob(table2d,
rows = NULL,
cols = NULL,
theme = tt1),
xmin=0.5, xmax=5, ymin=0.75, ymax=1)

p2d

library(patchwork)

p2c+p2d

声明:文中观点不代表本站立场。本文传送门:https://eyangzhen.com/418330.html

(0)
联系我们
联系我们
分享本页
返回顶部