数据准备
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))#设置中心圆名字
}
png('flower.png', w=500, h=500)
flower_plot(sample = st_u$strain, otu_num = st_u$Freq, core_otu =nrow(df_summary %>% filter(type=='core')),
start = 90, a = 0.7, b = 2, r = 0.8, ellipse_col = ellipse_col, circle_col = 'white')
dev.off()
当样本数少于5时，可以选择韦恩图表示，大于等于5时，建议使用upset图
library(UpSetR)
upsetdf <- as.data.frame(lapply(df, as.numeric))
upsetdf [upsetdf != 0] <- 1
upset(upsetdf, nset = 9,
order.by = 'freq')
