<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Goal2&nbsp;&nbsp;Var3&nbsp;&nbsp;Init4&nbsp;&nbsp;Load4.1&nbsp;&nbsp;Format4.2&nbsp;&nbsp;checks5&nbsp;&nbsp;Host tree subsample6&nbsp;&nbsp;PGLS6.1&nbsp;&nbsp;Digestive system6.1.1&nbsp;&nbsp;Saving output6.1.2&nbsp;&nbsp;Summary6.2&nbsp;&nbsp;Digestive system (simplified)6.2.1&nbsp;&nbsp;Saving output6.2.2&nbsp;&nbsp;Summary6.3&nbsp;&nbsp;Diet6.3.1&nbsp;&nbsp;Saving output6.3.2&nbsp;&nbsp;Summary6.3.3&nbsp;&nbsp;significant values6.3.3.1&nbsp;&nbsp;Saving results6.3.3.2&nbsp;&nbsp;Plotting7&nbsp;&nbsp;sessionInfo  Goal* PGLS to assess alpha diversity associations with diet or digestive systems  Var<jupyter_code>work_dir = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/modulating/PGLS/'

# microbiome
## 16S-arch phyloseq object (one sample per individual)
physeq_arch_IndD_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/merged/phyloseq/physeq_arch_IndD_clr-genus.RDS'

# diversity
alpha_div_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/merged/diversity-qiime2/alpha_div_genus.RDS'
alpha_div_hill_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S/LLA/diversity//alpha_div_iNEXT.RDS'

# host tree
host_tree_IndD_file = '/ebio/abt3_projects/Georg_animal_feces/data/animal/phylogeny/TT/unified_metadata_complete_200429_16S-arch_IndD.nwk'

# misc params
threads = 8
my_seed = 8394<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>library(dplyr)
library(tidyr)
library(ggplot2)

library(ape)
library(ade4)
library(caper)
library(phyloseq)

library(doParallel)
library(LeyLabRMisc)
df.dims()
set.seed(my_seed)
make_dir(work_dir)
RhpcBLASctl::blas_set_num_threads(1)<jupyter_output>Directory already exists: /ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/modulating/PGLS/<jupyter_text>Load<jupyter_code># phyloseq obj.
physeq = readRDS(physeq_arch_IndD_file)
physeq
# otu table
otu = physeq %>%
    phyloseq2df(otu_table) %>%
    as.data.frame
rownames(otu) = otu$OTU
otu$OTU = NULL
otu = otu %>% t %>% as.data.frame
otu
# metadata
metadata = physeq %>% 
    phyloseq2df(sample_data) %>%
    mutate(scientific_name = gsub(' ', '_', scientific_name),
           digestive_system = gsub('fermentation activity ', '', digestive_system),
           dig_sys = gsub('simple ', '', digestive_system)) %>%
    as.data.frame
rownames(metadata) = metadata$SampleID
metadata 
# alpha diversity calculated by QIIME2 (rarefaction)
alpha_div = readRDS(alpha_div_file) %>%
    do.call(cbind, .) %>%
    mutate(Sample = rownames(.))
rownames(alpha_div) = alpha_div$Sample
alpha_div
# # alpha diversity estimated by iNEXT
# alpha_div_hill = readRDS(alpha_div_hill_file)
# alpha_div_hill %>% names

# alpha_div_hill = alpha_div_hill$AsyEst
# alpha_div_hill
# # joining alpha div
# overlap(alpha_div$Sample, alpha_div_hill$Site)

# alpha_div = alpha_div %>%
#     inner_join(alpha_div_hill %>% 
#                    filter(Diversity == 'Shannon diversity') %>%
#                    dplyr::select(Site, Estimator) %>%
#                    rename('shannon_hill' = Estimator),
#                c('Sample'='Site')) %>%
#     as.data.frame

# rownames(alpha_div) = alpha_div$Sample
# alpha_div$Sample = NULL
# alpha_div
# host phylogeny
host_tree = read.tree(host_tree_IndD_file)
host_tree<jupyter_output><empty_output><jupyter_text>Format<jupyter_code>diet_vars = c(
'diet',
'diet_vertebrate',
'diet_invertebrate',
'diet_fruit',
'diet_flowers_nectar_pollen',
'diet_leaves_branches_bark_buds',
'diet_seeds',
'diet_grass_waterplants',
'diet_roots_tubers',
'digestive_system',
'dig_sys'
)
alpha_div = cbind(metadata[,diet_vars], alpha_div)
alpha_div
# encoding diet vars as numeric
alpha_div$diet = case_when(alpha_div$diet == 'herbivore' ~ 1,
                           alpha_div$diet == 'omnivore' ~ 0.5,
                           alpha_div$diet == 'carnivore' ~ 0)

alpha_div$digestive_system = case_when(alpha_div$digestive_system == 'hindgut' ~ 0,
                                 alpha_div$digestive_system == 'simple hindgut' ~ 0.25,
                                 alpha_div$digestive_system == 'simple' ~ 0.5,
                                 alpha_div$digestive_system == 'simple foregut' ~ 0.75,
                                 alpha_div$digestive_system == 'foregut' ~1)


alpha_div$dig_sys = case_when(alpha_div$dig_sys == 'hindgut' ~ 0,
                        alpha_div$dig_sys == 'simple' ~ 0.5,
                        alpha_div$dig_sys == 'foregut' ~1)

for(x in diet_vars[2:(length(diet_vars)-2)]){
    alpha_div[,x] = gsub(' .+', '', alpha_div[,x])
    alpha_div[,x] = ifelse(alpha_div[,x] == 'yes', 1, 0) %>% as.numeric
}

alpha_div$taxa = rownames(alpha_div)
alpha_div 
# checking values
colSums(alpha_div[,1:11])<jupyter_output><empty_output><jupyter_text>checks<jupyter_code>overlap(alpha_div$Sample, host_tree$tip.label)
# ordering
alpha_div = alpha_div[host_tree$tip.label,]
alpha_div<jupyter_output><empty_output><jupyter_text>Host tree subsample<jupyter_code># replicate analyses with different subsets of the host phylogeny
n_reps = 100
#' randomly selecting one per group
tree_subsample = function(L, df, tree){
    # get subsample (note: subsampling within each species)
    to_keep = df %>% 
        group_by(scientific_name) %>% 
        sample_n(1) %>%
        .$sample
    # subsampling tree
    to_rm = setdiff(tree$tip.label, to_keep)
    tree = drop.tip(tree, to_rm)
    return(tree)
}
# permuted trees
df = metadata %>%
    mutate(sample = SampleID) %>%
    dplyr::select(sample, scientific_name) 

doParallel::registerDoParallel(threads)
host_tree_l = plyr::llply(as.list(1:n_reps), 
                          function(x) tree_subsample(x, df, host_tree),
                          .parallel=TRUE)
host_tree_l %>% length
# check tips per tree
lapply(host_tree_l, function(x) x$tip.label %>% length) %>% unlist %>% summary<jupyter_output><empty_output><jupyter_text>PGLS<jupyter_code>pgls_each_var = function(y, x, tree, data){
    cdat = comparative.data(data=data, phy=tree, names.col="taxa")
    
    f = as.formula(paste(c(y, x), collapse=' ~ '))
    mod = pgls(f, cdat) %>% summary
    
    df = mod$coefficients %>% as.data.frame
    colnames(df) = c('Estimate', 'std_err', 't_value', 'p_value')
    df$R2 = mod$adj.r.squared %>% as.vector %>% .[1]
    df$x = x
    df$y = y
    df = df[2,] %>% as.data.frame
    return(df)
}
# #-- debug --#
# rep = 1
# host_tree = host_tree_l[[rep]]
# df = otu[host_tree$tip.label,] 

# pgls_each_var(y = '7efd5f90bad3210ac9182784e5882716',
#               x = 'diet',
#               tree = host_tree,
#               data = otu)
#' PGLS on each tree
pgls_each_tree = function(rep, y, x, tree_list, data, threads=1){
    tree = tree_list[[rep]]
    data = data[tree$tip.label,] 
    
    doParallel::registerDoParallel(threads)
    ret = plyr::llply(as.list(y),
                      pgls_each_var,
                      x = x,
                      tree = tree,
                      data = data,
                      .parallel = TRUE)
    ret = do.call(rbind, ret)
    ret$rep = rep
    return(ret)
}
#-- debug --#
y = colnames(alpha_div)[1:2]
ret = pgls_each_tree(1,
                     y = y,
                     x = 'digestive_system',
                     tree_list = host_tree_l,
                     data = alpha_div,
                     threads = threads)

ret = do.call(rbind, ret) %>%
    as.data.frame
colnames(ret) = y

ret 
# variables
alpha_div_labels = c('pielou_e', 'observed_otus', 'shannon', 'faith_pd')
alpha_div_labels<jupyter_output><empty_output><jupyter_text>Digestive system<jupyter_code># summary of covariate
metadata %>% .$digestive_system %>% table
ret = plyr::llply(1:length(host_tree_l),
                  pgls_each_tree,
                  y = alpha_div_labels,
                  x = 'digestive_system',
                  tree_list = host_tree_l,
                  data = alpha_div,
                  threads = threads)

ret = do.call(rbind, ret) %>%
    as.data.frame
ret 
send_email('http://127.0.0.1:9705/notebooks/notebooks/16S_arch/MiSeq-Runs116-122-126-189-190/03_modulating_factors/04d_PGLS_IndD_alpha-div.ipynb')
# formatting results
pgls_res = ret %>%
    group_by(rep) %>%
    mutate(p_value = p_value %>% as.Num,
           q_value = p.adjust(p_value, method='BH')) %>%
    group_by(y) %>%
    summarize(Estimate = mean(Estimate),
              std_err = mean(std_err),
              t_value = mean(t_value),
              p_value = mean(p_value),
              q_value = (length(q_value) - sum(q_value < 0.05)) / length(q_value),
              R2 = mean(R2),
              x = first(x),
              .groups='drop')

pgls_res$p_value %>% summary
pgls_res$q_value %>% summary
pgls_res<jupyter_output><empty_output><jupyter_text>Saving output<jupyter_code># output file name
pgls_res_file = file.path(work_dir, 'pgls_IndD_alpha-div_digSys.RDS')
pgls_res_file
# saving
saveRDS(pgls_res, file=pgls_res_file)
cat('File written:', pgls_res_file, '\n')
# re-load 
pgls_res = readRDS(pgls_res_file) 
pgls_res %>% dfhead<jupyter_output>[1] 4 9<jupyter_text>Summary<jupyter_code># summary
pgls_res$p_value %>% summary_x('p-value')
pgls_res$q_value %>% summary_x('adj p-value')
# p-value distribution
p = ggplot(pgls_res, aes(R2, p_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw()

p.dims(5,3.5)
plot(p)
# p-adj-value distribution
p = ggplot(pgls_res, aes(R2, q_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw() 

p.dims(5,3.5)
plot(p)<jupyter_output><empty_output><jupyter_text>Digestive system (simplified)<jupyter_code>ret = plyr::llply(1:length(host_tree_l),
                  pgls_each_tree,
                  y = alpha_div_labels,
                  x = 'dig_sys',
                  tree_list = host_tree_l,
                  data = alpha_div,
                  threads = threads)

ret = do.call(rbind, ret) %>%
    as.data.frame
ret 
# formatting results
pgls_res = ret %>%
    group_by(rep) %>%
    mutate(p_value = p_value %>% as.Num,
           q_value = p.adjust(p_value, method='BH')) %>%
    group_by(y) %>%
    summarize(Estimate = mean(Estimate),
              std_err = mean(std_err),
              t_value = mean(t_value),
              p_value = mean(p_value),
              q_value = (length(q_value) - sum(q_value < 0.05)) / length(q_value),
              R2 = mean(R2),
              x = first(x),
              .groups='drop')

pgls_res$p_value %>% summary
pgls_res$q_value %>% summary
pgls_res<jupyter_output><empty_output><jupyter_text>Saving output<jupyter_code># output file name
pgls_res_file = file.path(work_dir, 'pgls_IndD_alpha-div_digSysSimp.RDS')
pgls_res_file
# saving
saveRDS(pgls_res, file=pgls_res_file)
cat('File written:', pgls_res_file, '\n')
# re-load 
pgls_res = readRDS(pgls_res_file) 
pgls_res %>% dfhead<jupyter_output>[1] 4 9<jupyter_text>Summary<jupyter_code># summary
pgls_res$p_value %>% summary_x('p-value')
pgls_res$q_value %>% summary_x('adj p-value')
# p-value distribution
p = ggplot(pgls_res, aes(R2, p_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw()

p.dims(5,3.5)
plot(p)
# p-adj-value distribution
p = ggplot(pgls_res, aes(R2, q_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw() 

p.dims(5,3.5)
plot(p)<jupyter_output><empty_output><jupyter_text>Diet<jupyter_code># variables
alpha_div_labels = c('shannon', 'faith_pd')
alpha_div_labels
# summary of covariate
metadata %>% .$diet %>% table
# PGLS
ret = plyr::llply(1:length(host_tree_l),
                  pgls_each_tree,
                  y = alpha_div_labels,
                  x = 'diet',
                  tree_list = host_tree_l,
                  data = alpha_div,
                  threads = threads)

ret = do.call(rbind, ret) %>%
    as.data.frame
ret 
# formatting results
pgls_res = ret %>%
    group_by(rep) %>%
    mutate(p_value = p_value %>% as.Num,
           q_value = p.adjust(p_value, method='BH')) %>%
    group_by(y) %>%
    summarize(Estimate = mean(Estimate),
              std_err = mean(std_err),
              t_value = mean(t_value),
              p_value = mean(p_value),
              q_value = (length(q_value) - sum(q_value < 0.05)) / length(q_value),
              R2 = mean(R2),
              x = first(x))

pgls_res$p_value %>% summary
pgls_res$q_value %>% summary
pgls_res<jupyter_output><empty_output><jupyter_text>Saving output<jupyter_code># output file name
pgls_res_file = file.path(work_dir, 'pgls_IndD_alpha-div_diet.RDS')
pgls_res_file
# saving
saveRDS(pgls_res, file=pgls_res_file)
cat('File written:', pgls_res_file, '\n')
# re-load 
pgls_res = readRDS(pgls_res_file) 
pgls_res %>% dfhead
send_email('http://127.0.0.1:9705/notebooks/notebooks/16S_arch/MiSeq-Runs116-122-126-189-190/03_modulating_factors/04d_PGLS_IndD_alpha-div.ipynb')<jupyter_output>Sending email to: <EMAIL><jupyter_text>Summary<jupyter_code># summary
pgls_res$p_value %>% summary_x('p-value')
pgls_res$q_value %>% summary_x('adj p-value')
# p-value distribution
p = ggplot(pgls_res, aes(R2, p_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw()

p.dims(5,3.5)
plot(p)
# p-adj-value distribution
p = ggplot(pgls_res, aes(R2, q_value)) +
    geom_point(alpha=0.3) +
    scale_y_log10() +
    theme_bw() 

p.dims(5,3.5)
plot(p)<jupyter_output><empty_output><jupyter_text>significant values<jupyter_code># how many?
pgls_res %>%
    filter(q_value < 0.05) %>%
    .$y %>% unique_n('No. of OTUs')
# filtering & adding taxonomy
pgls_res_j = pgls_res %>%
    filter(q_value < 0.05) %>%
    inner_join(physeq %>% 
                   phyloseq2df(tax_table) %>%
                   mutate(OTU = gsub('^([0-9])', 'X\\1', OTU)),
               c('y' = 'OTU'))
pgls_res_j 
# viewing taxonomic distribution
df.dims(30)
pgls_res_j %>%
    dplyr::select(Estimate, p_value, q_value, R2,
                  Phylum, Class, Order, Family, Genus, Species) %>%
    arrange(-R2) 
df.dims()<jupyter_output><empty_output><jupyter_text>Saving results<jupyter_code># F = file.path(work_dir, 'PGLS_IndD_alpha-div_diet.RDS')
# saveRDS(pgls_res_j, file=F)
# cat('File written:', F, '\n')<jupyter_output><empty_output><jupyter_text>sessionInfo<jupyter_code>sessionInfo()
# adding taxonomy
pgls_res = pgls_res %>%
    inner_join(physeq %>% 
                   phyloseq2df(tax_table) %>%
                   mutate(OTU = gsub('^([0-9])', 'X\\1', OTU)),
               c('y' = 'OTU'))
pgls_res
p = pgls_res %>%
    mutate(y = y %>% as.factor %>% as.numeric %>% gsub('^', 'ASV-', .)) %>%
    mutate(Taxonomy = paste(Order, Family, Genus, Species, y, sep='; '),
           Taxonomy = stringr::str_wrap(Taxonomy, 50)) %>%
    ggplot(aes(Taxonomy, R2, color=Class)) +
    geom_boxplot() +
    labs(y='R^2') +
    theme_bw() 

p.dims(5,3)
plot(p)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14437751004016064, 'fraction_numerical': 0.02423025435073628, 'mean_word_length': 3.41910677314404, 'pattern_counts': {'":': 0, '<': 79, '<?xml version=': 0, '>': 160, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 4, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '221440', 'license': 'MIT', 'n_tokens_mistral': 6224, 'n_tokens_neox': 5797, 'n_words': 847, 'path': '03_modulating_factors/04a_PGLS_IndD_alpha-div.ipynb', 'repo_name': 'leylabmpi/16S-arc_vertebrate_paper', 'size': '14960'}