論文

Chromosome-level assemblies of multiple Arabidopsis genomes reveal hotspots of rearrangements with altered evolutionary dynamics

https://www.nature.com/articles/s41467-020-14779-y

擬南芥NC_panGenome.pdf

分析代碼的github主頁

https://github.com/schneebergerlab/AMPRIL-genomes

論文中組裝了7個擬南芥的基因組，做了一些泛基因組相關(guān)的分析凛俱，數(shù)據(jù)和大部分代碼都公開了炭玫，我們試著復(fù)現(xiàn)一下其中的圖和一些分析過程胚吁，今天的推文復(fù)現(xiàn)一下論文中的figure2d 下側(cè)的小圖鳖粟，泛基因組分析的論文里通常都有這個圖

image.png

這個展示的也是基因家族瞧壮，先用orthorfinder做聚類蝌麸，然后利用orthorfinder的結(jié)果進(jìn)行統(tǒng)計莉恼，作圖數(shù)據(jù)整理成了如下格式痊硕，數(shù)據(jù)有兩個赊级，一個是泛基因組的基因家族數(shù)量，一個是核心基因家族的數(shù)量岔绸，如何根據(jù)orthorfinder結(jié)果統(tǒng)計得到這個表格理逊，論文對應(yīng)的github主頁也提供了相應(yīng)的腳本，今天主要介紹畫圖

論文中提供的代碼是用R語言的基礎(chǔ)繪圖函數(shù)做的亭螟，這里我們用ggplot2來作圖

首先是核心基因家族

可以用連續(xù)的點或者我也看到有用箱線圖做的

這里我用箱線圖挡鞍，看起來可能會好看一點

library(tidyverse)
library(ggplot2)

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

ggplot()+
  geom_boxplot(data=dat01,aes(x=factor(V1),y=V2),fill="#3ba889")

image.png

擬合模型并添加擬合曲線

xvalue_core<-dat01 %>% pull(V1)
yvalue_core<-dat01 %>% pull(V2)

model_core<-nls(yvalue_core~A*exp(B*xvalue_core)+C,
                start = list(A=800,B=-0.3,C=800))
summary(model_core)

dat_core<-data.frame(x=seq(1,8,by=0.1),
                     y=predict(model_core,newdata = data.frame(xvalue_core=seq(1,8,by=0.1))))

ggplot()+
  geom_boxplot(data=dat01,aes(x=factor(V1),y=V2),fill="#3ba889")+
  geom_line(data=dat_core,aes(x=x,y=y),
            lty="dashed")

image.png

泛基因組

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

xvalue_pan<-dat02 %>% pull(V1)
yvalue_pan<-dat02 %>% pull(V2)

model_pan<-nls(yvalue_pan~A*exp(B*xvalue_pan)+C,
             start = list(A=800,B=-0.3,C=800))
model_pan
summary(model_pan)
dat_pan<-data.frame(x=seq(1,8,by=0.1),
                    y=predict(model_pan,newdata = data.frame(xvalue_pan=seq(1,8,by=0.1))))

合起來作圖和美化

ggplot()+
  geom_boxplot(data=dat01,aes(x=factor(V1),y=V2),fill="#3ba889")+
  geom_line(data=dat_core,aes(x=x,y=y),
            lty="dashed")+
  geom_boxplot(data=dat02,aes(x=factor(V1),y=V2),fill="#4593c3")+
  geom_line(data=dat_pan,aes(x=x,y=y),
            lty="dashed")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x="Number of genomes",y="Pan-enome size\n(gene number)") -> p1

ggplot()+
  geom_boxplot(data=dat01,aes(x=factor(V1),y=V2),fill="#f18e0c")+
  geom_line(data=dat_core,aes(x=x,y=y),
            lty="dashed")+
  geom_boxplot(data=dat02,aes(x=factor(V1),y=V2),fill="#af2934")+
  geom_line(data=dat_pan,aes(x=x,y=y),
            lty="dashed")+
  theme_bw()+
  theme(panel.grid = element_blank())+
  labs(x="Number of genomes",y="Pan-enome size\n(gene number)") -> p2

p2

library(patchwork)
p1+p2

image.png

這里有個問題是nls()函數(shù)擬合的時候會有一個start參數(shù)，這個參數(shù)里的值怎么確定预烙，暫時沒有想明白

示例數(shù)據(jù)和代碼可以給推文點贊墨微，然后點擊在看，最后留言獲取

歡迎大家關(guān)注我的公眾號

小明的數(shù)據(jù)分析筆記本

小明的數(shù)據(jù)分析筆記本 公眾號 主要分享：1扁掸、R語言和python做數(shù)據(jù)分析和數(shù)據(jù)可視化的簡單小例子翘县；2、園藝植物相關(guān)轉(zhuǎn)錄組學(xué)谴分、基因組學(xué)锈麸、群體遺傳學(xué)文獻(xiàn)閱讀筆記；3牺蹄、生物信息學(xué)入門學(xué)習(xí)資料及自己的學(xué)習(xí)筆記忘伞！

微信公眾號好像又有改動，如果沒有將這個公眾號設(shè)為星標(biāo)的話，會經(jīng)常錯過公眾號的推文氓奈，個人建議將 小明的數(shù)據(jù)分析筆記本 公眾號添加星標(biāo)翘魄，添加方法是

點開公眾號的頁面，右上角有三個點

image.png

點擊三個點舀奶，會跳出界面

image.png

直接點擊 設(shè)為星標(biāo) 就可以了