暂无图片
暂无图片
暂无图片
暂无图片
暂无图片

使用R语言画花瓣图(flower plot)

陈桂森 2021-12-29
1155

国际惯例先放图

在表示多个集合之间的关系,常常使用venn图进行可视化,但对于5个以上的集合,venn图在视觉感受中就很难提取有用的信息,这也是很多R包最多只提供画4~5个集合的原因。

对于基因组功能注释中,常常使用多个数据库,近期使用Pfam、interPro、GO、KEGG、Nr和SwissProt数据库对东星斑的基因组进行了功能注释,正好使用这个数据画花瓣图。

准备:
#install.packages(plotrix)
library(plotrix) #画图椭圆

#颜色范围
ellipse_col <- c('#6181BD4E','#F348004E','#64A10E4E','#9300264E','#464E044E','#049a0b4E','#4E0C664E','#D000004E','#FF6C004E','#FF00FF4E','#c7475b4E','#00F5FF4E','#BDA5004E','#A5CFED4E','#f0301c4E','#2B8BC34E','#FDA1004E','#54adf54E','#CDD7E24E','#9295C14E')

#定义画flower_plot的函数:引用自:https://www.maimengkong.com/kyjc/585.html
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))
    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),
            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))
}


画图:
anno = read.csv('function_annotion/anno.xls',sep='\t',header = T)

#提取各个数据库注释到的基因集
Nr = as.vector(na.omit(anno[,c('gene_id','nr')])[,c('gene_id')])
SwissProt = as.vector(na.omit(anno[,c('gene_id','swissprot')])[,c('gene_id')])
interPro = as.vector(na.omit(anno[,c('gene_id','interprosacn')])[,c('gene_id')])
KEGG = as.vector(na.omit(anno[,c('gene_id','KEGG')])[,c('gene_id')])
GO = as.vector(na.omit(anno[,c('gene_id','GO')])[,c('gene_id')])
Pfam = as.vector(na.omit(anno[,c('gene_id','pfam')])[,c('gene_id')])

#计算各个数据库基因集的交集元素的个数
core = length(Reduce(intersect,list(Nr,SwissProt,KEGG,GO,interPro,Pfam)))

#计算各个数据库基因集的元素个数
Nr=length(Nr)
SwissProt=length(SwissProt)
interPro=length(interPro)
KEGG=length(KEGG)
GO=length(GO)
Pfam=length(Pfam)
data = t(data.frame(Nr,SwissProt,interPro,KEGG,GO,Pfam)) #整合到数据框

#画花瓣图
flower_plot(sample =row.names(data), otu_num = data[,1], core_otu = core, start = 90, a = 1.5, b = 2, r = 0, ellipse_col = ellipse_col, circle_col = 'white')
#core_otu 交集元素个数
#a和b椭圆的半径
#r中心圆的半径
#circle_col中心圆的元素
#ellipse_col花瓣的颜色

文章转载自陈桂森,如果涉嫌侵权,请发送邮件至:contact@modb.pro进行举报,并提供相关证据,一经查实,墨天轮将立刻删除相关内容。

评论