%load_ext rpy2.ipython

%%R

workDir = '/home/nick/notebook/SIPSim/dev/fullCyc/'
physeqDir = '/home/nick/notebook/fullCyc/data/MiSeq_16S/515f-806r/V4_Lib1-7/phyloseq/'
physeqBulkCore = 'bulk-core'
physeqSIP = 'SIP-core_unk'

"""
Explanation: General info on the fullCyc dataset (as it pertains to SIPSim validation)

Simulating 12C gradients
Determining if simulated taxon abundance distributions resemble the true distributions
Simulation parameters to infer from dataset:
Infer total richness of bulk soil community 
richness of starting community


Infer abundance distribution of bulk soil community
NO: distribution fit
INSTEAD: using relative abundances of bulk soil community


Get distribution of total OTU abundances per fraction
Number of sequences per sample



User variables
End of explanation
"""


%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(phyloseq)
library(fitdistrplus)
library(sads)

%%R

dir.create(workDir, showWarnings=FALSE)

"""
Explanation: Init
End of explanation
"""


%%R
# bulk core samples
F = file.path(physeqDir, physeqBulkCore)
physeq.bulk = readRDS(F)
physeq.bulk.m = physeq.bulk %>% sample_data
physeq.bulk

%%R
# SIP core samples
F = file.path(physeqDir, physeqSIP)
physeq.SIP = readRDS(F)
physeq.SIP.m = physeq.SIP %>% sample_data
physeq.SIP

"""
Explanation: Loading phyloseq datasets
End of explanation
"""


%%R
physeq.bulk.12C = prune_samples(physeq.bulk.m$Substrate == '12C-Con', physeq.bulk) %>%
    filter_taxa(function(x) sum(x) > 0, TRUE)
physeq.bulk.12C.m = physeq.bulk.12C %>% sample_data
physeq.bulk.12C

%%R
physeq.SIP.12C = prune_samples(physeq.SIP.m$Substrate == '12C-Con', physeq.SIP) %>%
    transform_sample_counts(function(x) x/sum(x))
physeq.SIP.12C

%%R
# number of fractions per 12C day
as.Num = function(x) as.numeric(as.character(x))

physeq.SIP.12C %>% sample_data %>% 
    as.matrix %>% as.data.frame %>% 
    group_by(Day) %>%
    summarize(n = n()) %>%
    mutate(Day = as.Num(Day)) %>%
    arrange(Day)

"""
Explanation: Selecting gradients to simulate
End of explanation
"""


%%R -w 700

p = plot_richness(physeq.bulk.12C, measures=c('Observed', 'Chao1'))
p = p + theme(text = element_text(size=18))
p

%%R
ggplot_build(p)$data[[2]] %>% arrange(PANEL, x)

"""
Explanation: Infer total richness of bulk soil community
End of explanation
"""


%%R
df.OTU = physeq.bulk.12C %>%
    transform_sample_counts(function(x) x/sum(x)) %>%
    otu_table %>%
    as.matrix %>% as.data.frame 

df.OTU$OTU = rownames(df.OTU)
df.OTU = df.OTU %>% 
    gather('sample', 'abundance', 1:(ncol(df.OTU)-1))

df.OTU %>% head(n=3)

%%R -w 700 -h 350


days = physeq.bulk.12C.m$Day %>% unique %>% as.character
days = reorder(days, days %>% as.Num)

df.OTU.l = list()
for(day in days %>% sort){
    day %>% print
    rx = paste(c('\\.D', day, '\\.'), collapse='')
    df.OTU.l[[day]] = df.OTU %>% 
                filter(grepl(rx, sample), abundance > 0)
    plotdist(df.OTU.l[[day]]$abundance)
}


%%R -w 450 -h 400

lapply(df.OTU.l, function(x) descdist(x$abundance, boot=1000))

%%R
tmp1 = df.OTU.l[['1']]$abundance

f.norm = fitdist(tmp1, 'norm')
f.exp = fitdist(tmp1, 'exp')
f.logn = fitdist(tmp1, 'lnorm')
f.gamma = fitdist(tmp1, 'gamma')
f.beta = fitdist(tmp1, 'beta')

f.list = list(f.norm, f.exp, f.logn, f.gamma, f.beta)

plot.legend = c('normal', 'exponential', 'lognormal', 'gamma', 'beta')

par(mfrow = c(2,1))
denscomp(f.list, legendtext=plot.legend)
qqcomp(f.list, legendtext=plot.legend)

%%R
gofstat(f.list, fitnames=plot.legend)

%%R
summary(f.beta)

%%R
tmp1 %>% summary %>% print
rbeta(10000, 0.605, 1206.999) %>% summary

"""
Explanation: Infer abundance distribution of bulk soil community

distribution fit
End of explanation
"""


%%R
physeq.SIP.12C = prune_samples(physeq.SIP.m$Substrate == '12C-Con', physeq.SIP)
physeq.SIP.12C

%%R
df.OTU = physeq.SIP.12C %>% 
    otu_table %>%
    as.matrix %>% as.data.frame 

df.OTU$OTU = rownames(df.OTU)
df.OTU = df.OTU %>% 
    gather('sample', 'abundance', 1:(ncol(df.OTU)-1))

df.OTU %>% head(n=3)

%%R -h 300

df.OTU.s = df.OTU %>%
    group_by(sample) %>%
    summarize(total_abundance = sum(abundance))

cat('Number of samples (fractions): ', nrow(df.OTU.s), '\n')

p = ggplot(df.OTU.s, aes(total_abundance)) +
    geom_histogram(binwidth=100, aes(y = ..density..)) +
    geom_density() + 
    theme_bw() +
    theme(
        text = element_text(size=16)
        )
p

"""
Explanation: Notes

Instead of trying to simulate OTU abundances based on a distribution, I will use the observed abundances of each OTU

Get distribution of total OTU abundances per fraction

This provides the parameters for simulating OTU table subsampling (ie., sequencing simulation)
End of explanation
"""


%%R -w 700 -h 350

plotdist(df.OTU.s$total_abundance)

%%R -w 450 -h 400
descdist(df.OTU.s$total_abundance, boot=1000)

%%R
f.n = fitdist(df.OTU.s$total_abundance, 'norm')
f.ln = fitdist(df.OTU.s$total_abundance, 'lnorm')
f.ll = fitdist(df.OTU.s$total_abundance, 'logis')
f.list = list(f.n, f.ln, f.ll)

plot.legend = c('normal', 'log-normal', 'logistic')

par(mfrow = c(2,1))
denscomp(f.list, legendtext=plot.legend)
qqcomp(f.list, legendtext=plot.legend)

%%R

gofstat(list(f.n, f.ln, f.ll), fitnames=plot.legend)

%%R
summary(f.ln)

%%R -h 300
rlnorm(1000, 9.432, 0.5) %>% hist(breaks=30)

"""
Explanation: Notes

Total abundances of each sample
End of explanation
"""


