数据分析:“散点图:菌群与代谢物的相关性”

2020-10-23  本文已影响0人  生信学习者2

散点图:菌群与代谢物的相关性

散点图能较好地展示两坐标的相关趋势,再通过线性回归可计算两指标的相关关系,最终结合散点图以及相关关系结果共同在图上展示其结果。更多知识分享请到 https://zouhua.top/

导入所需R包

library(dplyr)
library(tibble)
library(ggplot2)
library(patchwork)

grp <- c("Case", "Control")
grp.col <- c("#283891", "#ED1C24")

读取数据

所需数据已经上传至百度网盘: https://pan.baidu.com/s/1nu1AoJd2YOfKEPWLxr3wpg 提取码:请邮件 zouhua1@outlook.com

phen <- read.csv("phenotype.csv")
metabolite <- read.csv("metabolite.csv") %>% column_to_rownames("X")
phylum <- read.csv("phylum.csv") %>% column_to_rownames("phylum")

函数

cor_plot_fun <- function(meta="Putrescine", 
                         taxo="p__Bacteroidetes"){
 
  mdat <- phen %>% inner_join(metabolite %>% 
                               rownames_to_column("SampleID"), 
                             by="SampleID") %>%
          tidyr::gather(key="metabolism", value="MetaScore", 
                        -c("SampleID", "Group")) %>%
          inner_join(phylum %>% t() %>% data.frame() %>% 
                       select(taxo) %>%
                       setNames("TaxScore") %>% 
                       rownames_to_column("SampleID"), 
                     by="SampleID")
  
  xlab <- paste(meta, "activity score")
  ylab <- paste(taxo, "Abundance")
  mdat_cln <- mdat %>% filter(metabolism%in%meta) %>%
              mutate(Group=factor(Group, levels = grp))

  dat_lab <- c()
  fr <- levels(mdat_cln$Group)
  for(i in 1:length(fr)){
    mdat_cln_fr <- mdat_cln %>% filter(Group%in%fr[i])
    fit <- lm(formula = TaxScore ~ MetaScore, data = mdat_cln_fr)
    lm.lst <- list(a = as.numeric(format(coef(fit)[1], digits = 4)),
              b = as.numeric(format(coef(fit)[2], digits = 4)),
              r2 = format(summary(fit)$r.squared, digits = 4),
              p = format(summary(fit)$coefficients[2,4], digits = 4))
    eq <- substitute(italic(y) == a + b %.% italic(x)~","~italic(R)^2~"="~r2~","~italic(P)~"="~p, lm.lst)
    lab <- substitute(~italic(R)^2~"="~r2~","~italic(P)~"="~p, lm.lst)
    
    MScore <- with(mdat_cln_fr, (max(MetaScore)-min(MetaScore))/2+min(MetaScore))
    TScore <- with(mdat_cln_fr, max(TaxScore)-(max(TaxScore)-min(TaxScore))/8)
    
    temp <- data.frame(MetaScore=MScore,
                       TaxScore=TScore,
                       Group=fr[i],
                       label=as.character(as.expression(lab)))
    dat_lab <- rbind(dat_lab, temp)
  }
  
    
  pl <- ggplot(mdat_cln, aes(x=MetaScore, y=TaxScore))+ 
    geom_point()+
    stat_smooth(method='lm', formula = y~x, color="red")+
    geom_text(data=dat_lab, aes(x=MetaScore, y=TaxScore, label = label),
                size = 4, parse = TRUE) +
    labs(x=xlab, y=ylab)+
    facet_wrap(~Group, scales = "free")+  
    theme_bw()+
    theme(axis.title = element_text(face = 'bold',color = 'black',size = 10),
          axis.text.y = element_text(color = 'black',size = 8),
          axis.text.x = element_text(color = 'black',size = 8),
          strip.text = element_text(size = 10, face = "bold"))
  
  return(pl)  
}

运行函数

Putrescine_plot <- cor_plot_fun(meta = "Putrescine", taxo="p__Bacteroidetes")
Eicosenoic__plot <- cor_plot_fun(meta = "Eicosenoic_acid", taxo="p__Bacteroidetes")
Glyceric_plot <- cor_plot_fun(meta = "Glyceric_acid", taxo="p__Bacteroidetes")

(Putrescine_plot / Eicosenoic__plot / Glyceric_plot) +
  plot_layout(ncol = 1)

结果显示,p__Bacteroidetes与三类代谢物均不存在显著相关关系。

sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    
system code page: 936

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] patchwork_1.0.1 ggplot2_3.3.2   tibble_3.0.3    dplyr_1.0.0    

loaded via a namespace (and not attached):
 [1] rstudioapi_0.11  knitr_1.29       magrittr_1.5     tidyselect_1.1.0 munsell_0.5.0   
 [6] colorspace_1.4-1 R6_2.4.1         rlang_0.4.7      tools_4.0.2      grid_4.0.2      
[11] gtable_0.3.0     xfun_0.18        withr_2.2.0      htmltools_0.5.0  ellipsis_0.3.1  
[16] yaml_2.2.1       digest_0.6.25    lifecycle_0.2.0  crayon_1.3.4     purrr_0.3.4     
[21] vctrs_0.3.2      glue_1.4.1       evaluate_0.14    rmarkdown_2.4    compiler_4.0.2  
[26] pillar_1.4.6     generics_0.0.2   scales_1.1.1     pkgconfig_2.0.3 

引用

参考文章如引起任何侵权问题,可以与我联系,谢谢。

上一篇 下一篇

猜你喜欢

热点阅读