一、先介绍一下选择性清除和图片的解读
选择消除
在自然选择下，一种新的正向选择或有利性状会在群体中增加频率(普遍率)，这是适应性进化的主要驱动力。一个性状要经历正向选择，必须有两个显著特征。首先，这种特性必须是有益的；换句话说，它必须增加有机体生存和繁殖的可能性。第二，这种特征必须是可遗传的，这样才能遗传给生物体的后代。如果一个特征导致更多的后代拥有相同的特征，那么这个特征比随机产生的特征更可能在群体中变得普遍。在分子水平上，当一种特定的DNA变异由于对携带它的生物体的影响而变得更加普遍时，选择就发生了。
简单的说就是基因组某区域由于受到了选择而消除多态性，即遗传多样性降低，在群体中出现高频的等位基因和低频的等位基因。图1显示了在选择前后沿一条染色体的多态性，包括受选择基因。祖先的等位基因（Ancestral alleles）显示为灰色，衍生的等位基因（Non-ancestral alleles）显示为蓝色。当一个新的正向选择基因(红色)上升到高频率时，染色体上邻近的连锁等位基因也随之上升到高频率，形成了选择消除（Selective sweep）。
图片
图1 物种驯化过程中多态性的变化趋势
选择消除分析
选择消除分析有助于我们对物种自然历史的理解。在驯化物种中，通常可以利用对品种或品种间表型多样性的理解，结合选择消除分析，来更好地理解控制表型变异的基本生物学意义。例如，在狗(所有陆生脊椎动物中体型多样性最大的物种)中，研究人员试图在最小和最大的品种之间进行选择消除分析。发现与大型犬种相比，小型犬种几乎只拥有一种单一基因。该基因编码胰岛素样生长因子1(IGF1)，已有研究表明它同样影响着小鼠和人类的体型。类似的研究设计已应用于几种不同的驯化物种，如牛、猪、羊、兔子和水稻。
如何进行选择消除分析研究呢？目前选择消除分析常用的分析方法主要有以下几种：
1. 基于群体分化分析（F-statistics）
群体间遗传分化系数（Fst），由F统计量演变而来，反映群体等位基因杂合性水平，用于衡量种群分化程度，取值从0到1。为0则认为两个种群间是随机交配的，基因型完全相似；为1则表示是完全隔离的，完全不相似。它往往从基因的多样性来估计，比如SNP或者微卫星(串联重复序列一种，长度小于等于10bp)。估算群体间的遗传分化系数可以区分群体间和群体内相对遗传变异大小，是解释群体遗传变异程度的主要来源。
图片
H, 杂合度，杂合基因型占的比例；
HT种群期望杂合度，把所有亚群当作一个种群，利用Hardy-Weinberg平衡预计种群的等位基因频率，计算杂合度：
图片
Hs亚群预计杂合度，通过每个亚群预计的杂合度计算（假设3个亚群体）：
图片
Hexp1表示1号群体的期望杂合度，其余类推。
实际研究中，若群体 Fst 值在 0~0.05 之间，表明其各亚群间遗传分化很小，可以不考虑；若 Fst 值在 0.05~0.15 之间，为中度分化；若 Fst 值在 0.15~0.25 之间，则亚群间遗传分化较大；Fst值在0.25以上，则说明亚群间有很大的遗传分化。
图片
图2 群体Fst分布图图片说明：横坐标代表不同的染色体名称，纵坐标代表相应染色体窗口内Fst值，两条虚线代表两种选择阈值（top 5%或1%）。
2. 基于群体多态性分析（Pi）
从核苷酸多态性数据推断适应性突变是进化遗传学的一个主要研究目标，核苷酸多样性（nucleotide diversity）是指样本中所有可能匹配成对的序列间核苷酸位点差异百分比的平均值，用Pi或π表示。核苷酸多样性常用于衡量种群内或种群间的多样性，是遗传变异的一个量化值，通常与其它衡量种群多样性的数值，如期望杂合度等有关联。
-
人工选择的群体，遗传多样性相对单一，Pi值较小
-
野生群体遗传多样性大，Pi值比较大
-
单个群体内部基因型多样性（0-1），多样性越大，Pi越大。
图片
图3 群体Pi分布图图片说明：横坐标代表不同的染色体名称，纵坐标代表相应染色体窗口内Pi值，两条虚线代表两种选择阈值（top 5%或1%）。
3. 基于群体中性进化分析（中性检验）
Tajima Test（Tajima’s D）的目的是鉴定目标DNA序列在进化过程中是否遵循中性进化模型。在标准中性进化模型下，Tajima's D的理论值为零。如果Tajima's D值为正，表明存在大量的中等频率的等位基因，这可能是由于群体瓶颈效应、群体结构，或者平衡选择引起的。如果Tajima's D值为负，表明存在大量的低频等位基因位点。Tajima's D分析表示中性进化，越偏离0，受选择程度越高。
图片
图4 群体Tajima's D分布图图片说明：横坐标代表不同的染色体名称，纵坐标代表相应染色体窗口内Tajima's D值。
上述提到的各种方法均有利弊，例如，Tajima's D在检验正向选择的同时容易受到群体历史和背景选择的干扰，但是Tajima's D所检验的低频突变丰度的信号能够在选择发生位点被固定后持续较长的一段时间。因此同时利用两种或多种检验方法使它们的优缺点得以互补，从而能够确保筛选到受选择区域的准确性。
Fst和Pi结合是已被证实是一种很有效力地检测选择消除区域的方法，特别是在挖掘与生存环境密切相关的功能区时，往往可以得到较强的选择信号。
图片
图5 群体间Fst和Pi的选择消除分析图片说明：横轴为Pi ratio（Log10）值，纵坐标表示Fst值，分别对应上面的频率分布图和右侧的频率分布图，中部的点图则代表不同窗口内的相应的Fst和Pi比值。其中最上方绿色和红色区域为Pi选择出来的top 5%区域，右侧橘色区域为Fst所选择top 5%区域，中间红色和紫色区域为Fst和Pi的交集，即为候选的位点。
参考文献
1、Kliman R M. Encyclopedia of evolutionary biology[M]. Academic Press, 2016.2、Tishkoff S A, Reed F A, Ranciaro A, et al. Convergent adaptation of human lactase persistence in Africa and Europe[J]. Nature genetics, 2007, 39(1): 31-40.3、Sutter N B, Bustamante C D, Chase K, et al. A single IGF1 allele is a major determinant of small size in dogs[J]. Science, 2007, 316(5821): 112-115.4、Hufford M B, Xu X, Van Heerwaarden J, et al. Comparative population genomics of maize domestication and improvement[J]. Nature genetics, 2012, 44(7): 808-811.5、Qiu Q, Wang L, Wang K, et al. Yak whole-genome resequencing reveals domestication signatures and prehistoric population expansions[J]. Nature communications, 2015, 6(1): 1-7.6、Li M, Tian S, Jin L, et al. Genomic analyses identify distinct patterns of selection in domesticated pigs and Tibetan wild boars[J]. Nature genetics, 2013, 45(12): 1431-1438.7、Posbergh C J, Huson H J. All sheeps and sizes: a genetic investigation of mature body size across sheep breeds reveals a polygenic nature[J]. Animal Genetics, 2020.8、Qanbari S, Pausch H, Jansen S, et al. Classic selective sweeps revealed by massive sequencing in cattle[J]. PLoS Genet, 2014, 10(2): e1004148.
接下来介绍组合图片的绘制方法！
二、再介绍计算方法和绘图代码
背景
小编：一见钟情这种图，搜了几篇文章，但没找到现成的R包可以画，各个社群里这一年时间遇到了4次同样的问题，不少人想学。图片
应通神邀请，复现一下选择清除里常用的一种Fst-θπ图。示例数据使用pinfsc50包的vcf文件：https://github.com/knausb/pinfsc50，下载地址：https://github.com/knausb/pinfsc50/blob/master/inst/extdata/pinf_sc50.vcf.gz。
图片
上游计算
这里我们随机把样本分为两群，编写popA和popB两个群体文件
#popA
BL2009P4_us23
DDR7602
IN2009T1_us22
LBUS5
NL07434
P10127
P10650
P11633
#popB
P12204
P13527
P1362
P13626
P17777us22
P6096
P7722
RS2009P1_us8
blue13
t30-4
然后使用vcftools进行计算。
FST
vcftools --gzvcf pinf_sc50.vcf.gz \
--weir-fst-pop popA \
--weir-fst-pop popB \
--out popA_vs_popB \
--fst-window-size 10000 \
--fst-window-step 1000
结果文件
# popA_vs_popB.windowed.weir.fst
CHROM BIN_START BIN_END N_VARIANTS WEIGHTED_FST MEAN_FST
Supercontig_1.50 1 10000 34 0.011472 0.0207269
Supercontig_1.50 1001 11000 19 -0.0161601 0.00668657
Supercontig_1.50 2001 12000 4 -0.0320845 -0.0166505
Supercontig_1.50 29001 39000 6 -0.153846 -0.153846
Supercontig_1.50 30001 40000 63 -0.0305726 -0.0340977
Supercontig_1.50 31001 41000 72 0.0460198 0.000469097
Supercontig_1.50 32001 42000 90 0.0348269 -0.001657
Supercontig_1.50 33001 43000 99 0.0280698 -0.00348489
Supercontig_1.50 34001 44000 103 0.0366487 0.00265645
pi
注意窗口和步长需要和上面一样长
vcftools \
--gzvcf pinf_sc50.vcf.gz \
--window-pi 10000 \
--window-pi-step 1000 \
--out popA \
--keep popA
vcftools \
--gzvcf pinf_sc50.vcf.gz \
--window-pi 10000 \
--window-pi-step 1000 \
--out popB
--keep popB
结果文件
# popA.windowed.pi
CHROM BIN_START BIN_END N_VARIANTS PI
Supercontig_1.50 1 10000 28 0.000728333
Supercontig_1.50 1001 11000 16 0.000416667
Supercontig_1.50 2001 12000 1 3.25e-05
Supercontig_1.50 30001 40000 12 0.000337657
Supercontig_1.50 31001 41000 19 0.000420406
Supercontig_1.50 32001 42000 33 0.000848475
Supercontig_1.50 33001 43000 39 0.00100639
Supercontig_1.50 34001 44000 43 0.00108092
Supercontig_1.50 35001 45000 43 0.00108092
# popB.windowed.pi
CHROM BIN_START BIN_END N_VARIANTS PI
Supercontig_1.50 1 10000 19 0.000270648
Supercontig_1.50 1001 11000 12 0.000209009
Supercontig_1.50 2001 12000 5 5.68522e-05
Supercontig_1.50 29001 39000 6 6.31826e-05
Supercontig_1.50 30001 40000 62 0.000819953
Supercontig_1.50 31001 41000 67 0.000880988
Supercontig_1.50 32001 42000 83 0.0011814
Supercontig_1.50 33001 43000 91 0.00128539
Supercontig_1.50 34001 44000 91 0.00128539
下面直接使用上述三个文件绘图。对于没有群体遗传分析经验的小伙伴我们也直接提供一份计算好的示例文件。
绘图
准备工作
主要是准备一个合适的主题还有预处理数据。预处理数据比较简单，就是将两个群体的π值合并到一个文件中，并计算一下Population A VS Population B和Population B VS Population A之间的比值。主题参考了ggtree::theme_dendrogram, 因此保留这个名字
library(tidyverse)
theme_dendrogram <- function(){theme_bw() +
theme(
legend.background = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.background = element_rect(fill = "white", colour = "white"),
panel.border = element_blank(),
axis.line = element_line(color = "black"),
axis.line.y = element_line(color = "black"),
axis.ticks.y = element_line(color = "black"),
axis.text.y = element_text(color = "black")
)
}
# 读入并清洗数据
pi_A <- read_tsv("popA.windowed.pi")
pi_B <- read_tsv("popB.windowed.pi")
fst <- read_tsv('popA_vs_popB.windowed.weir.fst')
pi <- pi_A %>%
left_join(pi_B, by=c("CHROM", "BIN_START", "BIN_END")) %>%
dplyr::rename(pi_A=PI.x, pi_B=PI.y) %>%
left_join(fst, by=c("CHROM", "BIN_START", "BIN_END")) %>%
select(-c("N_VARIANTS.x", "N_VARIANTS.y", "N_VARIANTS", "WEIGHTED_FST")) %>%
mutate(AVSB=pi_A/pi_B,
BVSA=pi_B/pi_A)
# 计算一下几个阈值，可以根据自己的数据调整
fst90 <- quantile(pi$MEAN_FST, 0.90)
pi90 <- quantile(pi$AVSB, 0.90, na.rm = TRUE)
pi10 <- quantile(pi$AVSB, 0.10, na.rm = TRUE)
# 计算下xlimits和ylimits，保证每张图的limit是一致的，方便后期拼图
limit_x <- c(min(log(pi$AVSB), na.rm=TRUE)*0.95, max(log(pi$AVSB), na.rm=TRUE)*1.05)
limit_y <- c(min(pi$MEAN_FST, na.rm=TRUE)*0.95, max(pi$MEAN_FST, na.rm=TRUE)*1.05)
主图
主图是最简单的一个，直接根据计算好的值画就行
p_main <- ggplot(pi, aes(x=log(AVSB), y=MEAN_FST)) +
geom_point(
aes(
color=case_when(
log(AVSB) > log(pi90) & MEAN_FST > fst90 ~ 'Selected region (B region)',
log(AVSB) < log(pi10) & MEAN_FST > fst90 ~ 'Selected region (A region)',
TRUE ~ 'Whole genome'
)
)
) +
geom_vline(xintercept=log(pi90), color='grey66', linetype="dashed", size=1) +
geom_vline(xintercept=log(pi10), color='grey66', linetype="dashed", size=1) +
geom_hline(yintercept=fst90, color='grey66', linetype="dashed", size=1) +
lims(x=limit_x, y=limit_y) +
labs(x=expr(theta[pi]*"ratio("*pi[italic('A')]/pi[italic('B')]*")"),
y=expr(italic("F"["ST"]))) +
scale_color_manual('', values = c('Selected region (A region)'='#39539d',
'Selected region (B region)'='#55ad57',
'Whole genome'='grey')) +
theme_bw() +
theme(panel.grid = element_blank(),
legend.position = c(0.75, 0.25),
legend.background = element_blank())
右边的图
图的难点主要有两个：
- 需要预先计算频率直方图。因为图中的频率直方图是着色的，而直接使用
geom_histogram()计算的不好着色，所以需要预先计算频率直方图 - 为了和累积分布函数的曲线图相匹配，从而使图更美观，频率直方图的Frequency需要缩放到
[0,1]的区间
其他细节见代码注释
# cut() function返回的值是区间，而为了在图上画柱状图，需要取其区间中间点作为坐标，因此编写此函数
get_midpoint <- function(cut_label) {
tmp_lst <- strsplit(gsub("\\(|\\)|\\[|\\]", "", as.character(cut_label)), ",")
sapply(tmp_lst, function(x)mean(as.numeric(x)))
}
#手工计算频率直方图，解决难点1
data_right <-pi %>%
# cut_width是根据区间长度划分并统计，还有对应的两个函数可以根据指定的bins数或指定区间
# 进行计算，可以查看cut_width的帮助文档
mutate(cut_group = cut_width(MEAN_FST, width = 0.01)) %>%
group_by(cut_group) %>%
summarize(counts=n()) %>%
mutate(prop=counts/sum(counts),
mid_pos=get_midpoint(as.character(cut_group))) %>%
mutate(plot_prop = 1/max(prop) * prop) # 解决难点2
# 便于绘制图例与坐标轴
max1<-max(data_right$prop)
fst90_str <- as.character(round(fst90, 3))
p_right <- ggplot(data_right) +
stat_ecdf(data=pi,
aes(y=MEAN_FST),
geom="smooth",
se=F,
size=0.5,
color="black") +
geom_col(aes(y = mid_pos,
x=plot_prop,
fill=ifelse(mid_pos>fst90,
'selected',
'other')
),
orientation = 'y') +
geom_hline(yintercept=fst90, color='grey66', linetype="dashed", size=1) +
lims(y=limit_y) +
scale_x_continuous("Cumulative (%)",
labels = scales::percent_format(suffix = ''),
position = "top",
sec.axis = sec_axis(~.*max1,
name="Frequency (%)",
labels = scales::percent_format(suffix = ''))
) +
scale_fill_manual('',
values = c('selected'= '#e15e37','other'='grey66'),
labels = c('selected'= expr(italic("F"["ST"])>!!fst90_str),
'other'=expr(italic("F"["ST"])<=!!fst90_str)
)
) +
theme_dendrogram() +
theme(axis.title.y=element_blank(),
axis.ticks.y=element_blank(),
axis.text.y=element_blank(),
axis.line.y = element_blank(),
legend.position = c(0.75, 0.25))
上边的图
如法炮制上面的图
data_top <-pi %>%
mutate(cut_group = cut_width(log(AVSB), width = 0.01)) %>%
group_by(cut_group) %>%
summarize(counts=n()) %>%
mutate(prop=counts/sum(counts),
mid_pos=get_midpoint(as.character(cut_group))) %>%
mutate(plot_prop = 1/max(prop) * prop)
max2<-max(data_top$prop)
pi90_str <- as.character(round(log(pi90), 3))
pi10_str <- as.character(round(log(pi10), 3))
p_top <- ggplot(data_top) +
stat_ecdf(data=pi,
aes(x=log(AVSB)),
geom="smooth",
se=F,
size=0.5,
color="black") +
geom_col(aes(x= mid_pos,
y=plot_prop,
fill=case_when(
mid_pos > log(pi90) ~ 'selectedB',
mid_pos < log(pi10) ~ 'selectedA',
TRUE ~ 'other'
)
)
) +
geom_vline(xintercept=log(pi90), color='grey66', linetype="dashed", size=1) +
geom_vline(xintercept=log(pi10), color='grey66', linetype="dashed", size=1) +
lims(x=limit_x) +
scale_y_continuous("Cumulative (%)",
labels = scales::percent_format(suffix = ''),
position = "right",
sec.axis = sec_axis(~.*max2,
name="Frequency (%)",
labels = scales::percent_format(suffix = ''))
) +
scale_fill_manual('',
values = c('selectedA'='#39539d',
'selectedB'='#55ad57',
'other'='grey'),
labels = c('selectedA'= expr(theta[pi]*"ratio"<=!!pi10_str),
'selectedB'= expr(theta[pi]*"ratio">=!!pi90_str),
'other'=expr(!!pi10_str<theta[pi]*"ratio<"*!!pi90_str))
) +
theme_dendrogram() +
theme(axis.title.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.x=element_blank(),
axis.line.x=element_blank(),
legend.position = c(0.75, 0.25)
)
拼图
使用patchwork拼图
library(patchwork)
p_top + plot_spacer() + p_main + p_right + plot_layout(ncol=2, heights = c(1,4), widths=c(4,1))
ggsave("FST.pdf", width=18, height=16, units="cm")
手工调整下缝隙（缝隙是坐标轴导致的，调整了半天没用代码搞成功，哪位大神会的话可以后台留言）
代码图：
图片
成品图：
图片
小编测试
测试成功！与正文数据的物种不同，细菌数据进行了双向选择，但是对小编的样本来说野生vs栽培有选择意义，体现了驯化过程，栽培vs野生没有意义，所以我的横坐标和正文的不一样，只有一条线，详细的复现的复现代码明天单独发一版推文。
图片
网友评论