当前位置： 首页 > news >正文

做搜狗手机网站优化,wordpress+discuz+seo,学编程有什么好处,可以做外贸私单的网站欢迎大家关注全网生信学习者系列#xff1a; WX公zhong号#xff1a;生信学习者Xiao hong书#xff1a;生信学习者知hu#xff1a;生信学习者CDSN#xff1a;生信学习者2 介绍 本教程使用基于R的函数来估计微生物群落的香农指数和丰富度#xff0c;使用MetaPhlAn prof…欢迎大家关注全网生信学习者系列 WX公zhong号生信学习者Xiao hong书生信学习者知hu生信学习者CDSN生信学习者2 介绍 本教程使用基于R的函数来估计微生物群落的香农指数和丰富度使用MetaPhlAn profile数据。估计结果进一步进行了可视化并与元数据关联以测试统计显著性。 数据 大家通过以下链接下载数据 百度网盘链接https://pan.baidu.com/s/1f1SyyvRfpNVO3sLYEblz1A提取码: 请关注WX公zhong号_生信学习者_后台发送 复现msm 获取提取码 R 包 SummarizedExperimentmiaggpubrggplot2lfe Alpha diversity estimation and visualization 使用alpha_diversity_funcs.R计算alpha多样性和可视化。 代码 SE_converter - function(md_rows, tax_starting_row, mpa_md) {# SE_converter function is to convery metadata-wedged mpa table into SummarisedExperiment structure.# md_rows: a vector specifying the range of rows indicating metadata.# tax_starting_row: an interger corresponding to the row where taxonomic abundances start.# mpa_md: a metaphlan table wedged with metadata, in the form of dataframe.md_df - mpa_md[md_rows,] # extract metadata part from mpa_md tabletax_df - mpa_md[tax_starting_row: nrow(mpa_md),] # extract taxonomic abundances part from mpa_md table### convert md_df to a form compatible with SummarisedExperiment ### SE_md_df - md_df[, -1]rownames(SE_md_df) - md_df[, 1]SE_md_df - t(SE_md_df)### convert md_df to a form compatible with SummarisedExperiment ###### prep relab values in a form compatible with SummarisedExperiment ###SE_relab_df - tax_df[, -1]rownames(SE_relab_df) - tax_df[, 1]col_names - colnames(SE_relab_df)SE_relab_df[, col_names] - apply(SE_relab_df[, col_names], 2, function(x) as.numeric(as.character(x)))### prep relab values in a form compatible with SummarisedExperiment ###SE_tax_df - tax_df[, 1:2]rownames(SE_tax_df) - tax_df[, 1]SE_tax_df - SE_tax_df[-2]colnames(SE_tax_df) - c(species)SE_data - SummarizedExperiment::SummarizedExperiment(assays list(relative_abundance SE_relab_df),colData SE_md_df,rowData SE_tax_df)SE_data }est_alpha_diversity - function(se_data) {# This function is to estimate alpha diversity (shannon index and richness) of a microbiome and output results in dataframe.# se_data: the SummarizedExperiment data structure containing metadata and abundance values.se_data - se_data |mia::estimateRichness(abund_values relative_abundance, index observed)se_data - se_data |mia::estimateDiversity(abund_values relative_abundance, index shannon)se_alpha_div - data.frame(SummarizedExperiment::colData(se_data))se_alpha_div }make_boxplot - function(df, xlabel, ylabel, font_size 11, jitter_width 0.2, dot_size 1, font_style Arial, stats TRUE, pal NULL) {# This function is to create a boxplot using categorical data.# df: The dataframe containing microbiome alpha diversities, e.g. shannon and observed with categorical metadata.# xlabel: the column name one will put along x-axis.# ylabel: the index estimate one will put along y-axis.# font_size: the font size, default: [11]# jitter_width: the jitter width, default: [0.2]# dot_size: the dot size inside the boxplot, default: [1]# font_style: the font style, default: Arial# pal: a list of color codes for pallete, e.g. c(#888888, #eb2525). The order corresponds the column order of boxplot.# stats: wilcox rank-sum test. default: TRUEif (stats) {nr_group length(unique(df[, xlabel])) # get the number of groupsif (nr_group 2) {group_pair list(unique(df[, xlabel]))ggpubr::ggboxplot(data df, x xlabel, y ylabel, color xlabel,palette pal, ylab ylabel, xlab xlabel,add jitter, add.params list(size dot_size, jitter jitter_width)) ggpubr::stat_compare_means(comparisons group_pair, exact T, alternative less) ggplot2::stat_summary(fun.data function(x) data.frame(y max(df[, ylabel]), label paste(Mean,mean(x))), geomtext) ggplot2::theme(text ggplot2::element_text(size font_size, family font_style))} else {group_pairs my_combn(unique((df[, xlabel])))ggpubr::ggboxplot(data df, x xlabel, y ylabel, color xlabel,palette pal, ylab ylabel, xlab xlabel,add jitter, add.params list(size dot_size, jitter jitter_width)) ggpubr::stat_compare_means() ggpubr::stat_compare_means(comparisons group_pairs, exact T, alternative greater) ggplot2::stat_summary(fun.data function(x) data.frame(y max(df[, ylabel]), label paste(Mean,mean(x))), geomtext) ggplot2::theme(text ggplot2::element_text(size font_size, family font_style))}} else {ggpubr::ggboxplot(data df, x xlabel, y ylabel, color xlabel,palette pal, ylab ylabel, xlab xlabel,add jitter, add.params list(size dot_size, jitter jitter_width)) ggplot2::theme(text ggplot2::element_text(size font_size, family font_style))} }my_combn - function(x) {combs - list()comb_matrix - combn(x, 2)for (i in 1: ncol(comb_matrix)) {combs[[i]] - comb_matrix[,i]}combs }felm_fixed - function(data_frame, f_factors, t_factor, measure) {# This function is to perform fixed effect linear modeling# data_frame: a dataframe containing measures and corresponding effects # f_factors: a vector of header names in the dataframe which represent fixed effects# t_factors: test factor name in the form of string# measure: the measured values in column, e.g., shannon or richness

all_factors - c(t_factor, f_factors)

for (i in all_factors) {

vars - unique(data_frame[, i])

lookup - setNames(seq_along(vars) -1, vars)

data_frame[, i] - lookup[data_frame[, i]]

}

View(data_frame)str1 - paste0(c(t_factor, paste0(f_factors, collapse )), collapse )str2 - paste0(c(measure, str1), collapse ~ )felm_stats - lfe::felm(eval(parse(text str2)), data data_frame)felm_stats

}加载一个包含元数据和分类群丰度的合并MetaPhlAn profile文件 mpa_df - data.frame(read.csv(./data/merged_abundance_table_species_sgb_md.tsv, header TRUE, sep \t))sampleP057P054P052…P049sexual_orientationMSMMSMMSM…Non-MSMHIV_statusnegativepositivepositive…negativeethnicityCaucasianCaucasianCaucasian…Caucasianantibiotics_6monthYesNoNo…NoBMI_kg_m2_WHOObeseClassIOverweightNormal…OverweightMethanomassiliicoccales_archaeon0.00.00.0…0.01322………………Methanobrevibacter_smithii0.00.00.0…0.19154 接下来我们将数据框转换为SummarizedExperiment数据结构以便使用SE_converter函数继续分析该函数需要指定三个参数 md_rows: a vector specifying the range of rows indicating metadata. Note: 1-based.tax_starting_row: an interger corresponding to the row where taxonomic abundances start.mpa_md: a metaphlan table wedged with metadata, in the form of dataframe. SE - SE_converter(md_rows 1:5,tax_starting_row 6, mpa_md mpa_df)SE class: SummarizedExperiment dim: 1676 139 metadata(0): assays(1): relative_abundance rownames(1676): Methanomassiliicoccales_archaeon|tSGB376GGB1567_SGB2154|tSGB2154 … Entamoeba_dispar|tEUK46681Blastocystis_sp_subtype_1|tEUK944036 rowData names(1): species colnames(139): P057 P054 … KHK16 KHK11 colData names(5): sexual_orientation HIV_status ethnicityantibiotics_6month BMI_kg_m2_WHO接下来我们可以使用est_alpha_diversity函数来估计每个宏基因组样本的香农指数和丰富度。 alpha_df - est_alpha_diversity(se_data SE) alpha_dfsexual_orientationHIV_statusethnicityantibiotics_6monthBMI_kg_m2_WHOobservedshannonP057MSMnegativeCaucasianYesObeseClassI1343.1847P054MSMpositiveCaucasianNoOverweight1412.1197……………………P052MSMpositiveCaucasianNoNormal1522.5273 为了比较不同组之间的alpha多样性差异我们可以使用make_boxplot函数并使用参数 df: The dataframe containing microbiome alpha diversities, e.g. shannon and observed with categorical metadata.xlabel: the column name one will put along x-axis.ylabel: the index estimate one will put along y-axis.font_size: the font size, default: [11]jitter_width: the jitter width, default: [0.2]dot_size: the dot size inside the boxplot, default: [1]font_style: the font style, default: Arialpal: a list of color codes for pallete, e.g. c(#888888, #eb2525). The order corresponds the column order of boxplot.stats: wilcox rank-sum test. default: TRUE shannon - make_boxplot(df alpha_df,xlabel sexual_orientation,ylabel shannon,stats TRUE,pal c(#888888, #eb2525))richness - make_boxplot(df alpha_df,xlabel sexual_orientation, ylabel observed,stats TRUE,pal c(#888888, #eb2525)) multi_plot - ggpubr::ggarrange(shannon, richness, ncol 2) ggplot2::ggsave(file shannon_richness.svg, plot multi_plot, width 4, height 5)通过固定效应线性模型估计关联的显著性 在宏基因组分析中除了感兴趣的变量例如性取向之外通常还需要处理多个变量例如HIV感染和抗生素使用。因此在测试微生物群落矩阵例如香农指数或丰富度与感兴趣的变量例如性取向之间的关联时控制这些混杂效应非常重要。在这里我们使用基于固定效应线性模型的felm_fixed函数该函数实现在R包lfe 中以估计微生物群落与感兴趣变量之间的关联显著性同时控制其他变量的混杂效应。 data_frame: The dataframe containing microbiome alpha diversities, e.g. shannon and observed with multiple variables.f_factors: A vector of variables representing fixed effects.t_factor: The variable of interest for testing.measure: The header indicating microbiome measure, e.g. shannon or richness lfe_stats - felm_fixed(data_frame alpha_df,f_factors c(c(HIV_status, antibiotics_6month)),t_factor sexual_orientation,measure shannon) summary(lfe_stats) Residuals:Min 1Q Median 3Q Max -2.3112 -0.4666 0.1412 0.5200 1.4137 Coefficients:Estimate Std. Error t value Pr(|t|)
(Intercept) 3.62027 0.70476 5.137 9.64e-07 *** sexual_orientationMSM 0.29175 0.13733 2.125 0.0355 *
HIV_statuspositive -0.28400 0.14658 -1.937 0.0548 .
antibiotics_6monthNo -0.10405 0.67931 -0.153 0.8785

antibiotics_6monthYes 0.01197 0.68483 0.017 0.9861

Signif. codes: 0 ‘’ 0.001 ‘’ 0.01 ‘’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.6745 on 134 degrees of freedom Multiple R-squared(full model): 0.07784 Adjusted R-squared: 0.05032 Multiple R-squared(proj model): 0.07784 Adjusted R-squared: 0.05032 F-statistic(full model):2.828 on 4 and 134 DF, p-value: 0.02725 F-statistic(proj model): 2.828 on 4 and 134 DF, p-value: 0.02725