最近堰氓，眾多學術(shù)期刊上開始出現(xiàn)一些讓人眼前一亮的高顏值小提琴圖适荣。在網(wǎng)上搜羅一番现柠，發(fā)現(xiàn)，規(guī)范的稱呼應該是split violin弛矛。直譯有點太抽象了够吩，暫且譯為對半小提琴吧。
該功能來自網(wǎng)址https://stackoverflow.com/questions/35717353/split-violin-plot-with-ggplot2丈氓。其中列舉了兩種方法周循，一種是建立了一個功能geo_split_violin直接嵌入ggplot中，無需過多設(shè)置万俗，無腦操作即可湾笛；另一種是手動打造，先分組求密度分布曲線闰歪，然后通過調(diào)整分組軸上的坐標位置來實現(xiàn)錯位嚎研，但是需要對數(shù)據(jù)進行轉(zhuǎn)換，而且面對更多分組的時候库倘，就需要特別的設(shè)置临扮，門檻略高。
以下示例基于自己的數(shù)據(jù)
加載工具包
library(reshape2)
library(dplyr)
library(rfPermute)
library(phyloseq)
library(lmerTest)
library(lme4)
library(car)
library(piecewiseSEM)
加載整理數(shù)據(jù)
load('../01_filter_and_rarefy_Bacteria/.RData')
remove(physeq)
sample_variables(rarefy)
rarefy = subset_samples(rarefy, NP<100&EC<2000&CN<150& (!Sequencing_ID_16s_18s%in% c('M24','M12','M30')))
map = data.frame(sample_data(rarefy))
map = mutate(map,Latitude = abs(Latitude),
group = case_when(Vegetation == 'Forest'~ 'Forest', TRUE ~ 'NonForest'))
map$group = as.factor(map$group)
新建向量存儲分組變量于樟，以便后文直接引用
envs = c('Latitude', "Elevation", "Slope", # geologic
"MAT", "AI", # climate
"Plant_cover", # plant
"pH", "EC", "Clay_silt", 'WHC','TP') # soil
funs = c("BG","PHOS","NAG","Rb","PO4","NO3", "NH4", 'NPP', 'ORC')
計算物種多樣性
richness = estimate_richness(rarefy, measures="Observed")
rownames(richness) == map$Global_Atlas_Order
map$Bacteria_richness = richness$Observed
map$Global_Atlas_Order = NULL
混合效應模型
forest_scale = mutate_if(forest, is.numeric, scale)
nonforest_scale = mutate_if(nonforest, is.numeric, scale)
# additive formular
formular00 =as.formula(
paste("Bacteria_richness~",
paste(c(envs, 'Forest','(1|CLIMA2)','(1|LandUse)'), collapse="+")))
formular01 =as.formula(
paste("Plant_richness~",
paste(c(envs, 'Forest','(1|CLIMA2)','(1|LandUse)'), collapse="+")))
lmer01 = lmer(formular00, data = mutate(rbind(forest_scale, nonforest_scale),
Forest = case_when(group == 'Forest' ~1, TRUE~0)))
lmer02 = lmer(formular01, data = mutate(rbind(forest_scale, nonforest_scale),
Forest = case_when(group == 'Forest' ~1, TRUE~0)))
混合效應模型的結(jié)果評估
anova(lmer01,type="I")
summary(lmer01)
vif(lmer01)
rsquared(lmer01)
混合效應模型的bootstrap 檢驗
lmer_boot01 = bootMer(lmer01,fixef,nsim = 1000,parallel = "snow")
lmer_boot02 = bootMer(lmer02,fixef,nsim = 1000,parallel = "snow")
lmer_boot0 = rbind(data.frame(melt(lmer_boot01$t), group = 'Bacteria'),
data.frame(melt(lmer_boot02$t), group = 'Plant'))
結(jié)果整理公条，因子排序，并設(shè)置出圖順序
lmer_boot0 = subset(lmer_boot0, Var2 != '(Intercept)', select = c(Var2, value, group))
lmer_boot0$Var2 = gsub('_',' ',lmer_boot0$Var2)
test0 = aggregate(value~group*Var2, lmer_boot0, median) %>% arrange(group, desc(value))
lmer_boot0$Var2 = factor(lmer_boot0$Var2, levels = test0$Var2[1:12])
兩個小花樣
#方便后文繪圖時迂曲，設(shè)置間隔性陰影背景來增加多個分組的區(qū)分度
odds0 <- seq(1, 12, 2)
rect0 = lmer_boot0[odds0,]
#文章開頭提到的那個***split_violin***函數(shù)內(nèi)容有點多靶橱，直接放在本空間里有點臃腫，選擇新建一個外掛直接調(diào)用路捧。
source('utiles.R')
繪圖
lmer_boot0_plt = ggplot(lmer_boot0, aes(Var2, value, fill = group)) +
geom_rect(data = rect0,
xmin = odds0 - 0.5,
xmax = odds0 + 0.5,
ymin = -Inf, ymax = +Inf,
fill = 'grey', alpha = 0.5, inherit.aes = F) +
geom_hline(yintercept = 0, lty = 2, color = 'black') +
geom_split_violin() +
theme_classic()+
labs(x = NULL, y = 'standardized effect on diversity',fill = NULL) +
theme(legend.position = c(0.8, 0.9),
legend.text = element_text(color = 'black'),
legend.background = element_blank(),
legend.key.size = unit(0.3,'cm'),
axis.text = element_text(color = 'black'),
axis.title = element_text(color = 'black'),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
lmer0.jpg
下面分頁箱線圖插入間隔行的陰影背景請查考這里
lmers.jpg
