R语言用指数函数拟合泛基因组曲线

泛基因组分析里泛基因家族曲线通常用幂函数（y=ax^b+c）来拟合，核心基因家族曲线通常是用指数函数（y=ae^bx+c）来拟合
论文 The pangenome of an agronomically important crop plant Brassica oleracea 中提到了这两个拟合方程
也有的论文里泛基因家族和核心基因家族曲线都用指数函数来拟合，比如论文 Chromosome-level assemblies of multiple Arabidopsis genomes reveal hotspots of rearrangements with altered evolutionary dynamics
拟南芥这篇论文里提供了拟合函数曲线的代码
https://github.com/schneebergerlab/AMPRIL-genomes/blob/master/pangenome/pan-genome.plot.r
这里用到的是R语言的nls() 函数，这个函数需要指定一个start函数，指定函数曲线ABC三个常数。这个最开始应该怎么选择这个初始值还没有想明白
最近看到一个R包 ggtrendline 用来画拟合曲线的图，非常方便，这个R包里也提供了拟合指数函数的方法，看了看这个R包的代码，也是用到的nls()函数，但是不用指定start函数，他这里应该是自定义了函数去算这个参数的初始值。
以下是R包中给出的示例代码

library(ggtrendline)
x<-1:5
y<-c(2,4,8,16,28)
xy<-data.frame(x,y)
getInitial(y ~ SSexp3P(x,a,b,c), data = xy)
拟合方程的代码
fitexp3P <- nls(y~SSexp3P(x,a,b,c), data=xy)
summary(fitexp3P)
这个R包还直接提供了作图函数，可以把拟合方程直接添加到图上
ggtrendline(x,y,model = “exp3P”,eSize = 5)+
theme_bw()

image.png
这里如果要改拟合方程的字体的细节，可以用如下代码
update_geom_defaults(“text”, list(family = “Times New Roman”,hjust=0,
vjust=0))
添加拟合方程文本的代码
param <- substitute(expression(italic(y) == 5~”e”^{8.25~italic(x)}~+10))[2]
ggplot()+
geom_text(aes(x=1,y=1,label=as.character(param)),parse = TRUE,size=10)

image.png
下面用拟南芥中的数据试试
泛基因家族曲线

library(tidyverse)
library(ggtrendline)

dat03<-read.table(“data/20230318/Source_Data.Figure2/Fig2d/pan-genome.gene.clustering.pan-genome.txt”,
header = FALSE)

dat03 %>% head()
x.pan<-dat03 %>% pull(V1)
y.pan<-dat03 %>% pull(V2)

ggtrendline(x.pan,y.pan,model = “exp3P”,
eSize = 10,
eq.y = 29200,
rrp.y = 29000,
CI.fill=”darkgreen”,
CI.color=”darkgreen”)+
theme_bw(base_size = 15)+
labs(x=NULL,y=NULL)+
theme(panel.grid = element_blank())+
ggtitle(“Pan”)

image.png
核心基因家族曲线

dat01<-read.table(“data/20230318/Source_Data.Figure2/Fig2d/pan-genome.gene.clustering.core-genome.txt”,
header = FALSE)

dat01 %>% head()

x.core<-dat01 %>% pull(V1)
y.core<-dat01 %>% pull(V2)

ggtrendline(x.core,y.core,model = “exp3P”,
eSize = 10,
eq.y = 27000,
rrp.y = 26800,
CI.fill=”darkgreen”,
CI.color=”darkgreen”)+
theme_bw(base_size = 15)+
labs(x=NULL,y=NULL)+
theme(panel.grid = element_blank())+
ggtitle(“Core”)

image.png
组合图

library(patchwork)
p1+p2

image.png