<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Goal2&nbsp;&nbsp;Var3&nbsp;&nbsp;Init4&nbsp;&nbsp;Load5&nbsp;&nbsp;Plotting distances6&nbsp;&nbsp;Testing for association via Mantel7&nbsp;&nbsp;Testing for association via Protest8&nbsp;&nbsp;Most similar/divergent via Procrustes residuals9&nbsp;&nbsp;sessionInfo  Goal* Determine whether bacterial and archaeal alpha diversity significantly correlate  Var<jupyter_code>work_dir = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/diversity-qiime2/'

# alpha diversity (calculated via qiime2)
beta_div_arc_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/diversity-qiime2_arc/beta_div_genus.RDS'
beta_div_bac_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/diversity-qiime2_bac/beta_div_genus.RDS'

# 16S-arch phyloseq object; rarefied
physeq_arc_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/physeq_arch_arc-bac_r250.RDS'
# 16S-bac phyloseq object; rarefied
physeq_bac_file = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/physeq_bac_arc-bac_r5k.RDS'

# params
threads = 8<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>library(dplyr)
library(tidyr)
library(ggplot2)
library(phyloseq)
library(ape)
library(vegan)
library(doParallel)
library(LeyLabRMisc)
df.dims()
make_dir(work_dir)<jupyter_output>Directory already exists: /ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/arch-bac/diversity-qiime2/<jupyter_text>Load<jupyter_code># loading phyloseq object
meta = readRDS(physeq_arc_file) %>%
    phyloseq2df(sample_data) %>%
    dplyr::select(SampleID, class, order, family, genus, scientific_name,
                  diet, habitat, sample_type)
meta 
# archaea beta div
beta_div_arc = readRDS(beta_div_arc_file)
beta_div_arc %>% names
# archaea beta div
beta_div_bac = readRDS(beta_div_bac_file)
beta_div_bac %>% names<jupyter_output><empty_output><jupyter_text>Plotting distances<jupyter_code>dist2df = function(D, L){
    m = L[[D]] %>% as.matrix 
    m[upper.tri(m, diag=TRUE)] = NA
    
    m %>% as.data.frame %>%
        mutate(X=rownames(.)) %>%
        gather(Y, value, -X) %>%
        filter(!is.na(value)) %>%
        mutate(Dist = D)
}

df = dist2df('jaccard', beta_div_arc)
df<jupyter_output><empty_output><jupyter_text>Testing for association via Mantel<jupyter_code>D = 'jaccard'
beta_div_arc[[D]] %>% as.matrix %>% dim
beta_div_bac[[D]] %>% as.matrix %>% dim
Mantel = function(D, L1, L2){
    res = vegan::mantel(L1[[D]], L2[[D]], method='spearman', permutations=999)
    return(c(D, res$statistic, res$signif))
}

doParallel::registerDoParallel(threads)
res = beta_div_arc %>%
    names %>% 
    plyr::llply(Mantel, L1=beta_div_arc, L2=beta_div_bac, .parallel=TRUE) %>%
    do.call(rbind, .) %>%
    as.data.frame
colnames(res) = c('measure', 'statistic', 'signif')
res = res %>%
    mutate(signif = signif %>% as.Num,
           padj = p.adjust(signif, method='BH'))
res
p = res %>%
    mutate(ss = statistic %>% as.Num) %>%
    mutate(measure = case_when(
        measure == 'bray_curtis' ~ 'Bray\nCurtis',
        measure == 'jaccard' ~ 'Jaccard',
        measure == 'weighted_unifrac' ~ 'Weighted\nUniFrac',
        measure == 'unweighted_unifrac' ~ 'Unweighted\nUniFrac')) %>%
    mutate(sig = padj < 0.05) %>%
    ggplot(aes(measure, ss, fill=sig)) +
    geom_bar(stat='identity') +
    scale_fill_discrete('Adj. P < 0.05') +
    labs(x='Diversity measure', y='Mantel rho') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1),
        legend.title = element_text(size=9)
    )

p.dims(4,3)
plot(p)<jupyter_output><empty_output><jupyter_text>Testing for association via Protest<jupyter_code>Protest = function(D, L1, L2){
    res = vegan::protest(L1[[D]], L2[[D]], permutations=999)
    return(c(D, res$signif, res$ss))
}

# protest in parallel
doParallel::registerDoParallel(threads)
proc_res = beta_div_arc %>%
    names %>% 
    plyr::llply(Protest, L1=beta_div_arc, L2=beta_div_bac, .parallel=TRUE) %>%
    do.call(rbind, .) %>%
    as.data.frame

# formatting
colnames(proc_res) = c('measure', 'signif', 'ss')
proc_res = proc_res %>%
    mutate(signif = signif %>% as.Num,
           padj = p.adjust(signif, method='BH'))
proc_res
# formatting
proc_res = proc_res %>%
    mutate(measure = case_when(measure == 'bray_curtis' ~ 'Bray-Curtis',
                               measure == 'jaccard' ~ 'Jaccard',
                               measure == 'unweighted_unifrac' ~ 'Unweighted UniFrac',
                               measure == 'weighted_unifrac' ~ 'Weighted UniFrac'))
proc_res
p = proc_res %>%
    mutate(ss = ss %>% as.Num) %>%
    mutate(sig = padj < 0.05) %>%
    ggplot(aes(measure, ss, fill=sig)) +
    geom_bar(stat='identity') +
    scale_fill_discrete('Adj. P < 0.05') +
    labs(x='Diversity measure', y='Protest: sum of squares') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(4,3.5)
plot(p)<jupyter_output><empty_output><jupyter_text>Most similar/divergent via Procrustes residuals<jupyter_code>Procrustes = function(D, L1, L2){
    res = vegan::procrustes(L2[[D]], L1[[D]])
    return(res)
}

proc_res1 = beta_div_arc %>%
    names %>% 
    lapply(Procrustes, L1=beta_div_arc, L2=beta_div_bac) 
names(proc_res1) = beta_div_arc %>% names
p.dims(6,5)
x = proc_res1 %>% lapply(plot)
# procrustes residuals
resids = proc_res1 %>% lapply(function(x) summary(x)$resid) %>%
    do.call(cbind, .) %>%
    as.data.frame %>%
    mutate(sample = rownames(.))
resids
# formatting
resids = resids %>%
    gather(measure, value, -sample) %>%
    inner_join(meta, c('sample'='SampleID')) %>%
    mutate(measure = case_when(measure == 'bray_curtis' ~ 'Bray-Curtis',
                               measure == 'jaccard' ~ 'Jaccard',
                               measure == 'unweighted_unifrac' ~ 'Unweighted UniFrac',
                               measure == 'weighted_unifrac' ~ 'Weighted UniFrac'))
resids
# plotting
p = resids %>%
    ggplot(aes(class, value, color=diet)) +
    geom_boxplot() +
    scale_color_discrete('Host diet') +
    labs(x='Host class', y='Procrustes residuals') +
    facet_wrap(~ measure, scales='free') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(7,5)
plot(p)<jupyter_output><empty_output><jupyter_text>sessionInfo<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1564002403846154, 'fraction_numerical': 0.029146634615384616, 'mean_word_length': 4.188620420888543, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '155226', 'license': 'MIT', 'n_tokens_mistral': 2809, 'n_tokens_neox': 2581, 'n_words': 374, 'path': '07_arc-bac/03b_compare_beta_div.ipynb', 'repo_name': 'leylabmpi/16S-arc_vertebrate_paper', 'size': '6656'}