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()

这里如果要改拟合方程的字体的细节，可以用如下代码

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)

下面用拟南芥中的数据试试

泛基因家族曲线

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")

核心基因家族曲线

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")

组合图

library(patchwork)
p1+p2

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

微信公众号好像又有改动，如果没有将这个公众号设为星标的话，会经常错过公众号的推文，个人建议将 小明的数据分析笔记本 公众号添加星标，添加方法是

点开公众号的页面，右上角有三个点

点击三个点，会跳出界面

直接点击 设为星标 就可以了

微信扫一扫分享文章