<jupyter_start><jupyter_text>03 .Ordination ==============  **Overview of relationships between samples and features**<jupyter_code>library(phyloseq)
suppressPackageStartupMessages(library(plyr))
suppressPackageStartupMessages(library(vegan))
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(ggplot2))<jupyter_output><empty_output><jupyter_text>**Low dimensional representation of relationships between samples/features**  We are going to follow closely phyloseq's ordination examples (https://joey711.github.io/phyloseq/plot_ordination-examples.html)  Read data:<jupyter_code>#rds <- readRDS("../../../../../data/tmp/microbiome/quant/qiime2/silva/physeq_filtered.rds")

rds <- readRDS(snakemake@input[[1]])
df <- data.frame(as(sample_data(rds), "matrix"))
group <- factor(df[,snakemake@params$condition])
cat("Condition: ", snakemake@params$condition, "\n")
cat("Summary:\n")
print(summary(group))
if (!is.null(snakemake@params$ref_level)){
    cat("Setting", snakemake@params$ref_level, "as reference level", "\n")
    group <- relevel(group, ref=snakemake@params$ref_level)
}
df[,"Sample_Group"] <- group
sample_data(rds) <- df

df <- data.frame(as(sample_data(rds), "matrix"))
batch <- NULL
if (!is.null(snakemake@params$batch)){
    print(cat("\nBatch factor given as: ", snakemake@params$batch, "\n"))
    batch <- factor(snakemake@params$batch)
    print(summary(batch))
    df[,snakemake@params$batch] <- batch
    sample_data(rds) <- df
}<jupyter_output>Condition:  Sample_Group 
Summary:
     HC     PPI PPI_FGP 
     11       9      16 
Setting HC as reference level 

Sample subset exracted from Sample_Biosource == saft 

We are now using 36 samples<jupyter_text>--------------- Next step is to agglomerate at a specific taxon level<jupyter_code>if (!is.null(snakemake@params$taxrank)){
    cat("Taxrank: ", snakemake@params$taxrank, "\n")
    if (require(speedyseq)){
        rds <- speedyseq::tax_glom(rds, taxrank=snakemake@params$taxrank)
    } else{
        rds <- tax_glom(rds, taxrank=snakemake@params$taxrank)
    } 
} else{
    cat("Taxrank: None, Data is not aggregated. Each feature is an ASV (specific sequence)", "\n")
}
rds<jupyter_output>Taxrank: None, Data is not aggregated. Each feature is an ASV (specific sequence)<jupyter_text>--------------- It is normal to filter a bit more aggresive in dimension reduction methods (as compared to diff. epxrs for instance) as we are looking for variance across all groups to identify global patterns. Also, we may work on a subset and thus prevalence thresholds may have changed.<jupyter_code>if (!is.null(snakemake@params$independent_filtering)){
    if (snakemake@params$independent_filtering == TRUE){
        wh0 = genefilter_sample(rds, filterfun_sample(function(x) x > 5), A=0.25*nsamples(rds))
        rds = prune_taxa(wh0, rds)
        phylum.sum = tapply(taxa_sums(rds), tax_table(rds)[, "Phylum"], sum, na.rm=TRUE)
        top6phyla = names(sort(phylum.sum, TRUE))[1:6]
        rds = prune_taxa((tax_table(rds)[, "Phylum"] %in% top6phyla), rds)
    }
}

cat("After filtering:\n")
rds<jupyter_output>After filtering:<jupyter_text>We like to start out with equal impact of all samples and start out with a normalization to equal count (Counts per million)<jupyter_code>rds = transform_sample_counts(rds, function(x) 1E6 * x/sum(x))
#https://github.com/LieberInstitute/jaffelab
cleaningY <- function(y, mod, P) {
    stopifnot(P <= ncol(mod))
    Hat <- solve(t(mod) %*% mod) %*% t(mod)
    ## For dealing with NAs
    ## https://stackoverflow.com/questions/16535084/matrix-multiplication-with-scattered-na-values
    ty <- t(y)
    ty[is.na(ty)] <- 0
    beta <- (Hat %*% ty)
    ## Note that y might still have the NAs, and NA - a number = NA
    ## so there's no need to reset the NAs back on cleany
    cleany <- y - t(as.matrix(mod[, -c(seq_len(P))]) %*% beta[-seq_len(P), ])
    return(cleany)
    }

if (!is.null(snakemake@params$block)){
    df <- data.frame(as(sample_data(rds), "matrix"))
    mod <- with(df, model.matrix(~ snakemake@params$condition + snakemake@params$batch))
    Y <- t(as(otu_table(rds), "matrix"))
    Yclean <- cleaningY(Y, mod, P=ncol(mod))
    otu_table(rds) <- Yclean
    }<jupyter_output><empty_output><jupyter_text>-------------------- Our default ordination methods is using Bray-Curtis distance (phylogenetic distance) and Non-negative multidimensional scaling (NMDS) as the dimension reduction method.<jupyter_code>GP.ord <- ordinate(rds, "NMDS", "bray")<jupyter_output>Square root transformation
Wisconsin double standardization
Run 0 stress 0.1997037 
Run 1 stress 0.2328027 
Run 2 stress 0.1997037 
... Procrustes: rmse 3.727728e-05  max resid 0.0001223678 
... Similar to previous best
Run 3 stress 0.233551 
Run 4 stress 0.2387247 
Run 5 stress 0.2083979 
Run 6 stress 0.210115 
Run 7 stress 0.199656 
... New best solution
... Procrustes: rmse 0.01588642  max resid 0.0838528 
Run 8 stress 0.2085257 
Run 9 stress 0.2145015 
Run 10 stress 0.2141249 
Run 11 stress 0.1997037 
... Procrustes: rmse 0.01588218  max resid 0.08386982 
Run 12 stress 0.3964296 
Run 13 stress 0.2367641 
Run 14 stress 0.2237325 
Run 15 stress 0.1997037 
... Procrustes: rmse 0.01589291  max resid 0.08385033 
Run 16 stress 0.2368311 
Run 17 stress 0.249399 
Run 18 stress 0.2085265 
Run 19 stress 0.2145014 
Run 20 stress 0.253712 
*** No convergence -- monoMDS stopping criteria:
    20: stress ratio > sratmax<jupyter_text>We start out by looking at the relations between taxa to look for any irregular patterns. Each point is either an OTU/ASV or summarized to a specific taxonomic level given by model name.<jupyter_code>p1 = plot_ordination(rds, GP.ord, type="taxa", color="Phylum", title="taxa")
print(p1)
p1 + facet_wrap(~Phylum, 3)<jupyter_output><empty_output><jupyter_text>Do the same figures for samples. Each point is one sample.<jupyter_code>p2 = plot_ordination(rds, GP.ord, type="samples", color="Sample_Group", shape="Sample_Group") 
p2 + geom_point(size=7) + ggtitle("samples")
p4 = plot_ordination(rds, GP.ord, type="split", color="Phylum", shape="Sample_Group", label="Sample_Group", title="biplot") 
p4<jupyter_output><empty_output><jupyter_text>Lets do a summary of all the available decompositions in phyloseq:<jupyter_code>dist = "bray"
ord_meths = c("DCA", "CCA", "RDA", "DPCoA", "NMDS", "MDS", "PCoA")
plist = llply(as.list(ord_meths), function(i, physeq, dist){
        ordi = ordinate(rds, method=i, distance=dist)
        plot_ordination(rds, ordi, "samples", color="Sample_Group")
}, rds, dist)
names(plist) <- ord_meths

pdataframe = ldply(plist, function(x){
    df = x$data[, 1:2]
    colnames(df) = c("Axis_1", "Axis_2")
    return(cbind(df, x$data))
})
names(pdataframe)[1] = "method"
p = ggplot(pdataframe, aes(Axis_1, Axis_2, color=Sample_Group, shape=Sample_Group, fill=Sample_Group))
p = p + geom_point(size=4)
p = p + facet_wrap(~method, scales="free")
p = p + scale_fill_brewer(type="qual", palette="Set1")
p = p + scale_colour_brewer(type="qual", palette="Set1")
p<jupyter_output><empty_output><jupyter_text>An example with weighted Unifrac (can be a bit more unstable than Bray) using Principal Coordinate Analysis<jupyter_code>ordu = ordinate(rds, "PCoA", "unifrac", weighted=TRUE)
p5 <- plot_ordination(rds, ordu, color="Sample_Group", shape="Sample_Group")
p5 = p5 + geom_point(size=7,alpha=0.8)
p5 = p5 + scale_colour_brewer(type="qual", palette="Set1")
p5 + ggtitle("MDS/PCoA on weighted-UniFrac distance")<jupyter_output><empty_output><jupyter_text>---------------------- Lets check for significance of our *condition* in the low dimensional space defined bu the ordination. Are the groups signficant different from eachother? For this we are using permuation ANOVA from the *vegan* package (adonis). We will extract the two first dimensions (PC) and use this for signficance testing.<jupyter_code>PC <- ordu$vectors[,1:2]
PC <- GP.ord$points
X <- cbind(PC, df)
adonis(PC ~ Sample_Group, data = X, method='eu')<jupyter_output><empty_output><jupyter_text>Write the default dimred coordinates to file<jupyter_code>write.table(X, file=snakemake@output[[1]])
cat("Wrote:", snakemake@output[[1]])<jupyter_output>Wrote: data/tmp/microbiome/analysis/GCF-2019-658-juice/ordinations.tsv
{'fraction_non_alphanumeric': 0.14123199806365727, 'fraction_numerical': 0.04199443301464359, 'mean_word_length': 4.919770773638969, 'pattern_counts': {'":': 0, '<': 76, '<?xml version=': 0, '>': 46, 'https://': 3, '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': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '245290', 'license': 'BSD-2-Clause', 'n_tokens_mistral': 3245, 'n_tokens_neox': 2920, 'n_words': 733, 'path': 'microbiome/rules/analysis/notebooks/03_phyloseq_ordination.ipynb', 'repo_name': 'gcfntnu/gcf-workflows', 'size': '8263'}