【绘图进阶】之lefse定制化绘图(五)

原创生信团上海天昊生物

lefse分析在微生物多样性中经常使用，但是由于数据库分类问题，我们会对其中一些没有意义的菌进行过滤，或挑选LDA值较高且P值显著的菌，重新绘制条形图。这里，我们展示了几种lefse的条形图。

0一

数据预处理

In [1]:
df = read.table('taxon.lefse.result.xls',sep = ' ',header = T)
head(df,4)

Out[1]:
Taxonomy Group LDA Pvalue
k__Archaea A 6.879678 0.03899022
k__Archaea.p__Euryarchaeota B 5.548066 0.03899022
k__Archaea.p__Euryarchaeota.c__Thermoplasmata B 5.497864 0.03899022
k__Archaea.p__Euryarchaeota.c__Thermoplasmata.o__Methanomassiliicoccales B 5.514124 0.02732372

如上图所示，我们拿到的lefse的文件格式一般是这个样子的，为了减少文字长度，我们通常只显示最后一个分类层级的名称，需要用“.” 作为分隔符，取最后一列。linux上使用awk -F '.' '{print $NF}' taxon.lefse.result.xls 获取。结果如下所示：

数据读取

In [2]:
df = read.table('lefse.result.xls',sep = ' ',header = T)
head(df)

Out[2]:
Taxonomy Group LDA Pvalue

k__Archaea A 6.879678 0.0389902

								2

p__Euryarchaeota B 5.548066 0.03899022
c__Thermoplasmata B 5.497864 0.03899022
o__Methanomassiliicoccales B 5.514124 0.02732372
f__Methanomassiliicoccaceae B 5.501945 0.02732372
g__Methanomassiliicoccus B 5.520169 0.02732372

二

Lefse绘图

In [3]:
library(ggplot2)

1.按照LDA值进行排序

In [4]:
ggplot(df,aes(x = reorder(Taxonomy,LDA),y=LDA,fill = Group)) +
geom_bar(stat = 'identity',colour = 'black',width = 0.9,position = position_dodge(0.7))+
xlab('') + ylab('LDA') + coord_flip() +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +theme(axis.text.y = element_blank(), axis.ticks = element_blank(), axis.text.x = element_text(face = 'bold'),axis.title.x = element_text(face='bold')) +
theme(panel.border = element_blank())+
geom_text(aes(y = ifelse(df$LDA >0,-0.1,0.1),label=Taxonomy),fontface=4,size=1,hjust = ifelse(df$LDA>0,1,0))

Out[4]:

2. 按照分组和LDA值进行排序

In [5]:
df=df[order(df[,2],df[,3]),]
head(df)

Out[5]:
Taxonomy Group LDA Pvalue
161 g__Oceanibaculum A 3.942752 0.02413390
147 s__Methylobacterium_oxalidis A 4.063120 0.03794154
71 s__Actinocorallia_longicatena A 4.134189 0.04593813
116 s__Paenibacillus_edaphicus A 4.138317 0.03899022
26 s__Aquihabitans_daechungensis A 4.149702 0.04964780
88 g__Segetibacter A 4.174600 0.03473526

In [6]:
df$Taxonomy = factor(df$Taxonomy,levels = as.character(df$Taxonomy))

In [7]:
df$Taxonomy = factor(df$Taxonomy,levels = as.character(df$Taxonomy))
ggplot(df,aes(x = Taxonomy,y=LDA,fill = Group)) +
geom_bar(stat = 'identity',colour = 'black',width = 0.9,position = position_dodge(0.7))+
xlab('') + ylab('LDA') + coord_flip() +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
theme(axis.text.y = element_blank(),axis.ticks = element_blank(),
axis.text.x = element_text(face = 'bold'),axis.title.x = element_text(face='bold')) +
theme(panel.border = element_blank())+
geom_text(aes(y = ifelse(df$LDA >0,-0.1,0.1),label=Taxonomy),fontface=4,size=1,hjust = ifelse(df$LDA>0,1,0))

Out[7]:

3. 带正负值的绘图

`将B组的LDA值改成负数`

In [8]:
df[df$Group=='A',3] = 0 - df[df$Group=='A',3]
head(df)

Out[8]:
Taxonomy Group LDA Pvalue
161 g__Oceanibaculum A -3.942752 0.02413390
147 s__Methylobacterium_oxalidis A -4.063120 0.03794154
71 s__Actinocorallia_longicatena A -4.134189 0.04593813
116 s__Paenibacillus_edaphicus A -4.138317 0.03899022
26 s__Aquihabitans_daechungensis A -4.149702 0.04964780
88 g__Segetibacter A -4.174600 0.03473526

In [9]:
df$Taxonomy = factor(df$Taxonomy,levels = as.character(df$Taxonomy))

In [10]:
ggplot(df,aes(x = Taxonomy,y=LDA,fill = Group)) +
geom_bar(stat = 'identity',colour = 'black',width = 0.9,position = position_dodge(0.7))+
xlab('') + ylab('LDA') + coord_flip() +
theme_bw() +
theme(panel.grid.major.y = element_blank()) +
theme(axis.text.y = element_blank(),axis.ticks = element_blank(),
axis.text.x = element_text(face = 'bold'),axis.title.x = element_text(face='bold')) +
theme(panel.border = element_blank()) +
geom_text(aes(y = ifelse(df$LDA >0,-0.1,0.1),label=Taxonomy),fontface=4,size=1,hjust = ifelse(df$LDA>0,1,0))

Out[10]:

4.绘制top10的正负区分条形图

In [11]:
df = rbind(head(df[df$Group=='A',],10),head(df[df$Group=='B',],10),head(df[df$Group=='C',],10))head(df)

Out[11]:
Taxonomy Group LDA Pvalue
161 g__Oceanibaculum A -3.942752 0.02413390
147 s__Methylobacterium_oxalidis A -4.063120 0.03794154
71 s__Actinocorallia_longicatena A -4.134189 0.04593813
116 s__Paenibacillus_edaphicus A -4.138317 0.03899022
26 s__Aquihabitans_daechungensis A -4.149702 0.04964780
88 g__Segetibacter A -4.174600 0.03473526

通过geom_text(size=4) 调节字体大小

In [12]:
df$Taxonomy = factor(df$Taxonomy,levels = as.character(df$Taxonomy))
ggplot(df,aes(x = Taxonomy,y=LDA,fill = Group)) +
geom_bar(stat = 'identity',colour = 'black',width = 0.9,position = position_dodge(0.7))+
xlab('') + ylab('LDA') + coord_flip() +
theme_bw() +
theme(axis.text.y = element_blank(),axis.ticks = element_blank(),
axis.text.x = element_text(face = 'bold'),axis.title.x = element_text(face='bold')) +
theme(panel.border = element_blank()) +
geom_text(aes(y = ifelse(df$LDA >0,-0.1,0.1),label=Taxonomy),fontface=4,size=4,hjust = ifelse(df$LDA>0,1,0))

Out[12]: