orthofinder分析泛基因组及可视化

数据准备

GeneCount文件包含了所有同源基因组数目不为1的所有基因组，在各个样本中的数目；UnassignedGenes文件则包含了仅包含一个基因的同源基因组的样本来源信息。首先，我们要先将两个数据合并，结果文件是一个包含所有基因组的表格。

library(tidyverse)
df <- read.csv('Orthogroups.GeneCount.tsv', sep='\t',row.names = 1)
df 

df_uniq <- read.csv('Orthogroups_UnassignedGenes.tsv',sep = '\t', row.names = 1)
df <- df[,1:(ncol(df)-1)]
df_uniq[df_uniq != ''] <- 1
df_uniq[df_uniq == ''] <- 0
df <- rbind(df, df_uniq)

核心基因和特有基因统计

核心基因的筛选思路为：使用apply函数按行分析，若一行中所有值都不为0，即表示该基因组在全部样本中出现，对应的行名为核心基因。
特有基因的筛选思路为：使用apply函数按行分析，若一行中仅有1个值不为0，则为特有基因，并提取不为0列的列名为拥有该基因的个体。
附属基因：除核心基因、特有基因外者，为附属基因。

df_summary <- data.frame(gene=rownames(df))
df_summary$type <- apply(df, 1, function(x) if (length(x[x!=0]) == 1) colnames(df)[which(x != 0)] else NA)
df_summary$type[apply(df, 1, function(x) all(x != 0))] <- 'core'
df_summary[is.na(df_summary)] <- 'accessory'

st_u <- data.frame(t(table(df_summary$type))) %>% rename(strain=Var2) %>%select(strain, Freq) %>%
  merge(., data.frame(strain = colnames(df)), all.y=T)
st_u[is.na(st_u)] <- 0

整理好的数据结构，df_summary包含每个基因组类型，若为特有基因，则标记所属样本。st_u进一步统计了每个样本包含的特有基因数目。

饼图展示核心基因、附属基因、特有基因各自的数目和比例。我们使用ggpubr来绘制。

library(ggpubr) 
pdata1 <- data.frame(group = factor(c('core','unique','accessory'), levels = c('core','unique','accessory')), 
                     freq = c(nrow(df_summary %>% filter(type=='core')),sum(st_u$Freq), nrow(df_summary %>% filter(type=='accessory'))))
piedata1 <- pdata1%>% 
  mutate(labs=paste0(freq,'\n(',100*round(freq/sum(freq),4),'%)')) 
ggpie(piedata1, 'freq', #绘图，只用写频数就行，切记不用再写分组 
      fill = 'group', palette = 'jco',#颜色板为jco. 
      label = 'labs', lab.pos = 'in', lab.font = c(4, 'white'))+ #设置标签，标签的位置在图的内部，标签的大小为4， 颜色为白色. 
  theme(legend.position = "top")+ labs(fill="")
ggsave('result/pie.png')

花瓣图包含了对核心基因和每个样本的特有基因的数目统计。我们基于plotrix包来完成。首先先导入r包和绘图功能参数，并指定颜色。
需要注意的是，指定的颜色是8位字符的哦，比起6位字符颜色，最后两位表示其透明度。

ellipse_col <- c('#6181BD4E','#F348004E','#64A10E4E','#9300264E','#464E044E','#049a0b4E','#daa5204E','#1e90ff4e','#ba55d34E')# 椭圆的颜色设置
flower_plot <- function(sample, otu_num, core_otu, start, a, b, r, ellipse_col, circle_col) {
  par( bty = 'n', ann = F, xaxt = 'n', yaxt = 'n', mar = c(1,1,1,1))#start为椭圆的起始位置, a为椭圆的短轴, b为椭圆的长轴, r为中心圆的半径
  plot(c(0,10),c(0,10),type='n')
  n   <- length(sample)
  deg <- 360 / n
  res <- lapply(1:n, function(t){
    draw.ellipse(x = 5 + cos((start + deg * (t - 1)) * pi / 180), #绘制椭圆
                 y = 5 + sin((start + deg * (t - 1)) * pi / 180), 
                 col = ellipse_col[t],
                 border = ellipse_col[t],
                 a = a, b = b, angle = deg * (t - 1))
    text(x = 5 + 2.5 * cos((start + deg * (t - 1)) * pi / 180),#在图上设置每个样品的otu_num数目
         y = 5 + 2.5 * sin((start + deg * (t - 1)) * pi / 180),
         otu_num[t])
    
    if (deg * (t - 1) < 180 && deg * (t - 1) > 0 ) {
      text(x = 5 + 3.3 * cos((start + deg * (t - 1)) * pi / 180),#设置每个样品名
           y = 5 + 3.3 * sin((start + deg * (t - 1)) * pi / 180),
           sample[t],
           srt = deg * (t - 1) - start,
           adj = 1,
           cex = 1
      )
    } else {
      text(x = 5 + 3.3 * cos((start + deg * (t - 1)) * pi / 180),
           y = 5 + 3.3 * sin((start + deg * (t - 1)) * pi / 180),
           sample[t],
           srt = deg * (t - 1) + start,
           adj = 0,
           cex = 1
      )
    }
  })
  draw.circle(x = 5, y = 5, r = r, col = circle_col, border = NA)#绘制中心圆
  text(x = 5, y = 5, paste('Core:', core_otu))#设置中心圆名字
}  

当样本数少于5时，可以选择韦恩图表示，大于等于5时，建议使用upset图

library(UpSetR)
upsetdf <- as.data.frame(lapply(df, as.numeric))
upsetdf [upsetdf != 0] <- 1
upset(upsetdf, nset = 9,
      order.by = 'freq')