%%R
otu.tbl.file1 = '/home/nick/notebook/SIPSim/dev/bac_genome1210/atomIncorp_taxaIncorp/0/10/1/OTU_n2_abs1e9_sub-norm_filt.physeq'
otu.tbl.file2 = '/home/nick/notebook/SIPSim/dev/bac_genome1210/atomIncorp_taxaIncorp/100/10/1/OTU_n2_abs1e9_sub-norm_filt.physeq'

physeq1 = readRDS(otu.tbl.file1)
physeq2 = readRDS(otu.tbl.file2)

%%R

ord1 = ordinate(physeq1, method='NMDS', distance='bray')
ord2 = ordinate(physeq2, method='NMDS', distance='bray')

ord1 %>% scores %>% head %>% print
ord2 %>% scores %>% head %>% print

%%R

get.fracs = function(ord){
    fracs = gsub('.+__', '', rownames(ord %>% scores)) %>% as.data.frame()
    colnames(fracs) = c('fractions')
    fracs = fracs %>% 
        separate(fractions, c('start','end'), sep='-', convert=T) %>%
        mutate(start = start * 1000,
               end = end * 1000)
    return(fracs)
    }

ord1.f = get.fracs(ord1)
ord2.f = get.fracs(ord2)

%%R
library(IRanges)

%%R

ord1.r = IRanges(start=ord1.f$start, end=ord1.f$end)
ord2.r = IRanges(start=ord2.f$start, end=ord2.f$end)

%%R

ov = findOverlaps(ord1.r, ord2.r, select='first')
ov 

%%R

ov = findOverlaps(ord1.r, ord2.r)
ov 

"""
Explanation: Mapping fractions between gradient communities in order to perform procrustes
End of explanation
"""


%%R
otu.tbl.file1 = '/home/nick/notebook/SIPSim/dev/bac_genome1210/atomIncorp_taxaIncorp/0/10/1/OTU_n2_abs1e9_sub-norm_filt.physeq'
otu.tbl.file2 = '/home/nick/notebook/SIPSim/dev/bac_genome1210/atomIncorp_taxaIncorp/100/10/1/OTU_n2_abs1e9_sub-norm_filt.physeq'

physeq1 = readRDS(otu.tbl.file1)
physeq2 = readRDS(otu.tbl.file2)

%%R

ord1 = ordinate(physeq1, method='NMDS', distance='bray')
ord2 = ordinate(physeq2, method='NMDS', distance='bray')

%%R

grps = as.character(rep(seq(1,nrow(ord1$points) / 2), 2))
grps = append(grps, '2')

plot(ord1, type = "p", display='sites')
elps = ordiellipse(ord1, grps, kind="se", conf=0.95, lwd=2, col="blue")

elps = elps %>% summary %>% t %>% as.data.frame
elps

%%R

ggplot(elps, aes(NMDS1, NMDS2)) +
    geom_point()

%%R

get.ellipse = function(ord){
    grps = as.character(rep(seq(1,nrow(ord$points) / 2), 2))
    grps = append(grps, '2')
    
    plot(ord, type = "p", display='sites')
    elps = ordiellipse(ord, grps, kind="se", conf=0.95, lwd=2, col="blue")

    elps = elps %>% summary %>% t %>% as.data.frame
    return(elps)
    }


get.ellipse(ord1)

%%R

BD.range = seq(1.6, 1.9, 0.004)

mid = function(x, y){ (x + y)/2 }

get.BD.range = function(tbl, BD.range){
    tbl = as.data.frame(tbl)
    tbl$lib = gsub('__.+', '', rownames(tbl)) %>% as.character
    tbl$BD.start = gsub('.+__([0-9.]+)-.+', '\\1', rownames(tbl)) %>% as.numeric
    tbl$BD.end = gsub('.+-', '', rownames(tbl)) %>% as.numeric
    tbl$BD.mid = mapply(mid, tbl$BD.start, tbl$BD.end)
    tbl$BD.bin = cut(tbl$BD.mid, breaks=BD.range)
    
    tbl = tbl %>%
        unite(group, lib, BD.bin, sep='__')
    
    return(tbl)
    }

ord.BD = get.BD.range(ord1 %>% scores, BD.range=BD.range)
ord.BD %>% head

%%R
?tidyr::unite

%%R
# making fixed BD-range & binning by BD.mid

BD.range = seq(1.6, 1.9, 0.004)
BD.range

%%R
ord.BD$BD.bin = cut(ord.BD$BD.mid, breaks=BD.range)
ord.BD %>% head

"""
Explanation: Calculating centroid of binned fraction samples

centroid of all 20 replicates for fraction samples that fall into the BD-range bin
trying oriellipse() function from vegan package
End of explanation
"""


