<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Goal2&nbsp;&nbsp;Var3&nbsp;&nbsp;Init4&nbsp;&nbsp;Load5&nbsp;&nbsp;Filtering dataset6&nbsp;&nbsp;Checks7&nbsp;&nbsp;Subsampling tree8&nbsp;&nbsp;Creating datasets9&nbsp;&nbsp;parafit: global test9.1&nbsp;&nbsp;Saving results9.2&nbsp;&nbsp;Summarizing10&nbsp;&nbsp;parafit: individual links10.1&nbsp;&nbsp;Saving results10.2&nbsp;&nbsp;Summarizing10.3&nbsp;&nbsp;Links per host11&nbsp;&nbsp;sessionInfo  Goal* parafit analysis to support PACo results  Var<jupyter_code>work_dir = '/ebio/abt3_projects/Georg_animal_feces/data/16S_arch/MiSeq-Runs-116-122-126-189-190/LLA/cophylo/'

# 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_rare250.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
threads = 12
my_seed = 8394<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(tidytable)

library(phyloseq)
library(ape)
library(ggpubr)

library(LeyLabRMisc)
clustermq_setup() 
library(clustermq)
df.dims()
make_dir(work_dir)
set.seed(my_seed)
rescale_dist_mtx = function(m){
    m = m %>% as.matrix
    labs = m %>% colnames
    n_row = m %>% nrow
    n_col = m %>% ncol
    x = m %>% as.vector 
    x = scales::rescale(x) 
    m = matrix(x, nrow=n_row, ncol=n_col)
    colnames(m) = labs
    rownames(m) = labs
    m = m %>% as.dist
    return(m)
}<jupyter_output><empty_output><jupyter_text>Load<jupyter_code># phyloseq
physeq_arch_IndD = readRDS(physeq_arch_IndD_file)
physeq_arch_IndD
# host tree
host_tree_IndD = read.tree(host_tree_IndD_file)
host_tree_IndD<jupyter_output><empty_output><jupyter_text>Filtering dataset<jupyter_code># filtering
physeq_arch_IndD = physeq_arch_IndD %>%
    filter_taxa(function(x) sum(x > 0) >= 2, TRUE) %>%
    prune_samples(sample_sums(.)>0, .)
physeq_arch_IndD
# metadata
metadata = physeq_arch_IndD %>%
    phyloseq2df(sample_data) %>%
    dplyr::select(SampleID, scientific_name) %>%
    as.data.frame

rownames(metadata) = metadata$SampleID
metadata[sample_names(physeq_arch_IndD),]
metadata<jupyter_output><empty_output><jupyter_text>Checks<jupyter_code>overlap(sample_names(physeq_arch_IndD), host_tree_IndD$tip.label)<jupyter_output>intersect(x,y): 185 
setdiff(x,y): 0 
setdiff(y,x): 0 
union(x,y): 185<jupyter_text>Subsampling tree* multiple subsamples: one sample per species<jupyter_code># number of subsampled trees
ntrees = 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)
}
# subsampling trees
df = metadata %>%
    mutate(sample = rownames(.)) %>%
    dplyr::select(sample, scientific_name) 

doParallel::registerDoParallel(threads)
host_tree_l = plyr::llply(as.list(1:ntrees), 
                          function(x) tree_subsample(x, df, host_tree_IndD),
                          .parallel=TRUE)

# tree lengths
host_tree_l %>% length %>% print
lapply(host_tree_l, function(x) x$tip.label %>% length) %>% unlist %>% summary
# checking overlap of labels
setdiff(host_tree_l[[1]]$tip.label, rownames(metadata)) %>% length %>% print<jupyter_output>[1] 0<jupyter_text>Creating datasets<jupyter_code>make_parafit_input = function(host_tree, physeq){
    # subsampled phyloseq object
    physeq_f = prune_samples(sample_names(physeq) %in% host_tree$tip.label, 
                             physeq) 

    # microbial tree distance matrix
    micro_D = physeq_f %>%
        phy_tree %>% cophenetic %>% 
        rescale_dist_mtx %>% as.matrix
    # host tree distance matrix
    host_D = host_tree %>% cophenetic %>% 
             rescale_dist_mtx %>% as.matrix
    
    # otu pres/abs matrix
    otu = physeq_f %>% 
        otu_table %>%
        as.data.frame %>%
        as.matrix %>%
        t %>%
        apply(2, function(x) ifelse(x > 0, 1, 0)) %>% 
        as.matrix                
              
              
    # checking overlap
    x = length(intersect(rownames(otu), rownames(host_D)))
    y = length(union(rownames(otu), rownames(host_D)))
    stopifnot(x == y)
              
    # preparing paco data
    return(list(host.D=host_D, para.D = micro_D, HP=otu))
}
# creating list of parafit input objects
doParallel::registerDoParallel(threads)
PF_l = plyr::llply(host_tree_l, make_parafit_input, physeq=physeq_arch_IndD, .parallel=TRUE)
PF_l %>% length<jupyter_output><empty_output><jupyter_text>parafit: global test<jupyter_code># Function for running parafit on each dataset subsample
parafit_each = function(L, nperm=99, seed=3843, correction="cailliez", test.links=TRUE){
    parafit(L$host.D, L$para.D, L$HP, 
            nperm=nperm, 
            seed=seed,
            correction=correction,
            test.links=test.links)
}

## debug (WARNING: Slow)
#res = lapply(PF_l, parafit_each, nperm=9, test.links=FALSE)
# cluster job submission
tmpl = list(job_time = '00:59:00', job_mem = '10G', 
            conda = '/ebio/abt3_projects/Georg_animal_feces/envs/phyloseq',
            log_file = clustermq_logfile())  
params = list(nperm = 999, test.links = FALSE)
PF_res_l = Q(parafit_each, PF_l, const=params, pkgs=c('ape'),
             n_jobs=100, job_size=1, template=tmpl)
# # Running on cluster
# ## cluster resources
# resources = list(h_rt = '00:59:00',
#                  h_vmem = '12G',
#                  threads = '1',
#                  conda.env = '/ebio/abt3_projects/Georg_animal_feces/envs/phyloseq', 
#                  conda.path = '/ebio/abt3_projects/software/miniconda3_gt4.4/bin')
# plan(batchtools_sge, resources=resources, workers=100)
# ## job submission
# PF_res_l = future_lapply(PF_l, parafit_each, 
#                          nperm=9, 
#                          test.links=FALSE,
#                          future.packages=c('ape'))
# results
PF_res_l %>% length<jupyter_output><empty_output><jupyter_text>Saving results<jupyter_code># for loading results instead of re-running
PF_file = file.path(work_dir, 'physeq_IndD-sens_parafit.RDS')
PF_file
# saving
saveRDS(PF_res_l, file=PF_file)
cat('File written:', PF_file, '\n')
# re-loading
PF_res_l = readRDS(PF_file)
PF_res_l %>% length
send_email('http://127.0.0.1:9805/notebooks/notebooks/16S_arch/MiSeq-Runs116-122-126-189-190/04_cophylo/03_parafit_rarefy.ipynb')<jupyter_output>Sending email to: <EMAIL><jupyter_text>Summarizing<jupyter_code># distribution of effect size
lapply(PF_res_l, function(x) x$ParaFitGlobal) %>% unlist %>% summary %>% print
# distribution of p-values
lapply(PF_res_l, function(x) x$p.global) %>% unlist %>% summary %>% print<jupyter_output>Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.00100 0.00100 0.00100 0.00106 0.00100 0.00500<jupyter_text>parafit: individual links<jupyter_code>parafit_each = function(L, nperm=99, seed=3843, correction="cailliez", test.links=TRUE){
    parafit(L$host.D, L$para.D, L$HP, 
            nperm=nperm, 
            seed=seed,
            correction=correction,
            test.links=test.links)
}
# # cluster jobs
# ## cluster resources
# resources = list(h_rt = '200:00:00',
#                  h_vmem = '10G',
#                  threads = '1',
#                  conda.env = 'py3_physeq_eco', 
#                  conda.path = "/ebio/abt3_projects/software/miniconda3_gt4.4/bin")
# plan(batchtools_sge, resources=resources, workers=100)
# ## submission
# PF_res_links_l = future_lapply(PF_l, parafit_each, 
#                          nperm=999, 
#                          test.links=TRUE,
#                          future.packages=c('ape'))
# cluster job submission
tmpl = list(job_time = '72:00:00', job_mem = '10G', 
            conda = '/ebio/abt3_projects/Georg_animal_feces/envs/phyloseq',
            log_file = clustermq_logfile())  
params = list(nperm = 999, test.links = TRUE)
PF_res_links_l = Q(parafit_each, PF_l, const=params, pkgs=c('ape'),
             n_jobs=100, job_size=1, template=tmpl)
PF_res_links_l = PF_res_l 
# results
PF_res_links_l %>% length<jupyter_output><empty_output><jupyter_text>Saving results<jupyter_code># for loading results instead of re-running
PF_file = file.path(work_dir, 'physeq_IndD-sens_parafit-links.RDS')
PF_file
# saving
saveRDS(PF_res_links_l, file=PF_file)
cat('File written:', PF_file, '\n')
# re-loading
PF_res_links_l = readRDS(PF_file)
PF_res_links_l %>% length
send_email('http://127.0.0.1:9805/notebooks/notebooks/16S_arch/MiSeq-Runs116-122-126-189-190/04_cophylo/03_parafit_rarefy.ipynb')<jupyter_output><empty_output><jupyter_text>Summarizing<jupyter_code># distribution of effect size
lapply(PF_res_links_l, function(x) x$ParaFitGlobal) %>% unlist %>% summary
# distribution of p-values
lapply(PF_res_links_l, function(x) x$p.global) %>% unlist %>% summary<jupyter_output><empty_output><jupyter_text>Links per host<jupyter_code># format
para_per_host = lapply(PF_res_links_l, function(x) x$para.per.host) %>%
    do.call(rbind, .) %>%
    as.data.frame %>%
    mutate(Dataset_rep = 1:nrow(.)) %>%
    gather(Sample, N_links, -Dataset_rep) 
               
para_per_host
# adding metadata
para_per_host = physeq_arch_IndD %>%
    sample_data %>%
    as.matrix %>%
    as.data.frame %>%
    dplyr::select(SampleID, class, order, family, genus, scientific_name, diet) %>%
    inner_join(para_per_host, c('SampleID'='Sample')) 

para_per_host 
# summarizing by species: check
para_per_host %>%
    group_by(scientific_name) %>%
    summarize(n=n()) %>%
    ungroup() %>%
    .$n %>% summary %>% print
# summarizing by species
alpha=0.95

para_per_host = para_per_host %>%
    group_by(scientific_name, class, order, family, genus, diet) %>%
    summarize(min_n_links = min(N_links),
              q1_n_links = quantile(N_links, 0.25),
              mean_n_links = mean(N_links),
              median_n_links = median(N_links),
              q3_n_links = quantile(N_links, 0.75),
              max_n_links = max(N_links),
              sd_n_links = sd(N_links),
              CI_low = quantile(N_links, (1-alpha)/2),
              CI_high = quantile(N_links, alpha+((1-alpha)/2))) %>%
    ungroup() %>%
    mutate(scientific_name = scientific_name %>% reorder(-median_n_links))

para_per_host$CI_low %>% summary %>% print
para_per_host 
# plotting distribution
p = ggplot(para_per_host, aes(scientific_name, median_n_links, color=class)) +
    geom_linerange(aes(ymin=median_n_links-sd_n_links,
                       ymax=median_n_links+sd_n_links)) +
    geom_point() +
    theme_bw() +
    theme(
        axis.text.x = element_blank()
    )

p.dims(10,4)
plot(p)
# plotting distribution
p = para_per_host %>%
    ggplot(aes(class, mean_n_links)) +
    geom_boxplot() +
    scale_y_log10() +
    labs(x='Host class', y='No. of host-symbiont links') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(4,3.5)
plot(p)
# plotting distribution
p = para_per_host %>%
    ggplot(aes(class, mean_n_links, color=diet)) +
    geom_boxplot() +
    scale_color_discrete('Diet') +
    labs(x='Class', y='No. of host-symbiont links') +
    theme_bw() +
    theme(
        axis.text.x = element_text(angle=45, hjust=1)
    )

p.dims(5,3.5)
plot(p)
# anova
aov_res = car::Anova(lm(mean_n_links ~ class + diet, data=para_per_host), type=2)
aov_res<jupyter_output><empty_output><jupyter_text>sessionInfo<jupyter_code>sessionInfo()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1345941820032228, 'fraction_numerical': 0.027563395810363836, 'mean_word_length': 3.8526748971193414, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 137, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, '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': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '81121', 'license': 'MIT', 'n_tokens_mistral': 4792, 'n_tokens_neox': 4422, 'n_words': 688, 'path': '04_cophylo/03_parafit.ipynb', 'repo_name': 'leylabmpi/16S-arc_vertebrate_paper', 'size': '11811'}