In [1]:
workDir = '/home/nick/notebook/SIPSim/dev/bac_genome10/'
SIPSimExe = '/home/nick/notebook/SIPSim/SIPSim'
In [2]:
import os
import sys
import numpy as np
import pandas as pd
import subprocess
In [3]:
%load_ext rpy2.ipython
%pylab inline
In [4]:
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
In [8]:
! cd $workDir; \
$SIPSimExe indexGenomes genomes/genomes10.txt \
--fp ./genomes/ --np 10 > index_log.txt
In [285]:
%%bash -s "$workDir" "$SIPSimExe"
# amplicon fragments
cd $1
$2 fragGC genomes/genomes10.txt \
--fp ./genomes/ --fr 515Fm-927Rm.fna \
--np 10 > genome10_ampFragGC.txt
In [286]:
%%bash -s "$workDir" "$SIPSimExe"
# shotgun fragments
cd $1
$2 fragGC genomes/genomes10.txt \
--fp ./genomes/ --np 10 \
> genome10_shotFragGC.txt
In [287]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 gradientComms genomes/genomes10.txt \
--fp ./genomes/ --pf grinder_profile \
> genome10_comm_n3.txt
In [288]:
import os
# making config file
config = """
[library 1]
# baseline: no incorp
[[intraPopDist 1]]
distribution = uniform
weight = 1
[[[start]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 0
end = 0
[[[end]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 0
end = 0
[library 2]
# split intra-populations
## to get some taxa with split; use inter-pop mixture for 2nd intra-pop mu, where mixture is highly uneven
[[intraPopDist 1]]
distribution = normal
weight = 0.5
[[[mu]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 90
sigma = 2
[[[sigma]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 5
sigma = 2
[[intraPopDist 2]]
distribution = normal
weight = 0.5
[[[mu]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 5
sigma = 2
[[[sigma]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 5
sigma = 2
[library 3]
# split inter-pop distribution (some approx. full; others none)
[[intraPopDist 1]]
distribution = normal
weight = 1
[[[mu]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 90
sigma = 2
# these taxa in the community get no incorp
[[[[interPopDist 2]]]]
distribution = uniform
start = 0
end = 0
[[[sigma]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 5
sigma = 2
"""
outfile = os.path.join(workDir, 'genome10_n3.config')
outf = open(outfile, 'wb')
outf.write(config)
outf.close()
In [289]:
import os
# making config file
config = """
[library 1]
# baseline: no incorp
[[intraPopDist 1]]
distribution = uniform
weight = 1
[[[start]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 0
end = 0
[[[end]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 0
end = 0
[library 2]
# full incorp
[[intraPopDist 1]]
distribution = uniform
weight = 1
[[[start]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 100
end = 100
[[[end]]]
[[[[interPopDist 1]]]]
distribution = uniform
start = 100
end = 100
[library 3]
# split inter-pop distribution (some approx. full; others none)
[[intraPopDist 1]]
distribution = normal
weight = 1
[[[mu]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 90
sigma = 2
# these taxa in the community get no incorp
[[[[interPopDist 2]]]]
distribution = uniform
start = 0
end = 0
[[[sigma]]]
[[[[interPopDist 1]]]]
distribution = normal
mu = 5
sigma = 2
"""
outfile = os.path.join(workDir, 'genome10_n3.config')
outf = open(outfile, 'wb')
outf.write(config)
outf.close()
In [290]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 isoIncorp genome10_comm_n3.txt genome10_n3.config > genome10_comm_n3_incorp.txt
In [291]:
%%R
library(ggplot2)
library(dplyr)
In [292]:
%%R -i workDir
infile = paste(c(workDir, 'genome10_comm_n3_incorp.txt'), collapse='/')
tbl = read.csv(infile, sep='\t')
In [293]:
%%R -w 700
tbl$taxon_name = reorder(tbl$taxon_name, tbl$param_value, max)
ggplot(tbl, aes(taxon_name, param_value, color=param)) +
geom_point() +
facet_grid(param ~ library, scales='free_y') +
theme(
text = element_text(size=16),
axis.text.x = element_text(angle=90, hjust=1),
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
In [294]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 fractions genome10_comm_n3.txt > genome10_comm_n3_fracs.txt
In [295]:
%%R -i workDir
infile = paste(c(workDir, 'genome10_comm_n3_fracs.txt'), collapse='/')
tbl = read.csv(infile, sep='\t')
In [296]:
%%R -w 700
#tbl$taxon_name = reorder(tbl$taxon_name, tbl$param_value, max)
tbl$library = as.character(tbl$library)
ggplot(tbl, aes(library, fraction_size)) +
geom_boxplot()
In [297]:
%%R -h 300
tbl_sum = group_by(tbl, library) %>%
summarize( n_fracs = n())
ggplot(tbl_sum, aes(library, n_fracs)) +
geom_bar(stat='identity')
In [298]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 OTU_table genome10_ampFragGC.txt genome10_comm_n3.txt \
genome10_comm_n3_incorp.txt genome10_comm_n3_fracs.txt \
--abs_abund 1e4 --log genome10_OTU_abnd1e4_log.txt \
> genome10_OTU_abnd1e4.txt
In [299]:
%%R -i workDir
# loading file
infile = paste(c(workDir, 'genome10_OTU_abnd1e4_log.txt'), collapse='/')
tbl = read.csv(infile, sep='\t')
# reformat
tbl$taxon = gsub("(.*?_.*?)_(.+)", "\\1\n\\2", tbl$taxon)
tbl$fragment_length = tbl$fragment_length / 1000
In [300]:
%%R -w 600 -h 1000
# plot
ggplot(tbl, aes(fragment_length, fragment_GC)) +
geom_density2d() +
facet_grid(taxon ~ library) +
labs(x="fragment length [kb]", y="fragment G+C") +
theme(
text = element_text(size=18),
axis.text.x = element_text(angle=90)
)
In [301]:
%%R -w 1000 -h 700
# plot
tbl.l1 = tbl %>% filter(library == 1)
ggplot(tbl.l1, aes(fragment_length, fragment_GC)) +
geom_density2d() +
facet_wrap( ~ taxon) +
labs(x="fragment length [kb]", y="fragment G+C") +
theme(
text = element_text(size=18),
axis.text.x = element_text(angle=90)
)
In [302]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 OTU_table genome10_shotFragGC.txt genome10_comm_n3.txt \
genome10_comm_n3_incorp.txt genome10_comm_n3_fracs.txt \
--abs_abund 1e5 > genome10_OTU_abnd1e6.txt
In [5]:
%%R -i workDir
# loading file
infile = paste(c(workDir, 'genome10_OTU_abnd1e6.txt'), collapse='/')
tbl = read.csv(infile, sep='\t')
In [6]:
%%R
# formatting table
tbl$BD_min = gsub('-.+', '', tbl$fractions)
tbl$BD_min = as.numeric(tbl$BD_min)
tbl$BD_max = gsub('.+-', '', tbl$fractions)
tbl$BD_max = as.numeric(tbl$BD_max)
In [7]:
%%R
# summarizing counts (should be approx. total abundance)
tbl %>%
group_by(library) %>%
summarize(sum(count))
In [8]:
%%R -w 800
# plotting absolute abundances
## BD for G+C of 0 or 100
BD.GCp0 = 0 * 0.098 + 1.66
BD.GCp100 = 1 * 0.098 + 1.66
## plot
ggplot(tbl, aes(BD_min, count, fill=taxon, group=taxon)) +
geom_area(stat='identity', alpha=0.5, position='dodge') +
geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
facet_grid(library ~ .) +
labs(x='Buoyant density') +
theme( text = element_text(size=16) )
In [9]:
%%R -w 800
# plotting relative abundances
ggplot(tbl, aes(BD_min, count, fill=taxon, group=taxon)) +
geom_area(stat='identity', alpha=0.8, position='fill') +
geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
facet_grid(library ~ .) +
labs(x='Buoyant density') +
theme( text = element_text(size=16) )
Notes:
In [48]:
%%R -w 500 -h 400
taxa.sel = unique(tbl$taxon)[1:5]
tbl.tmp = filter(tbl, library == 1, taxon %in% taxa.sel, count > 0)
x = as.character(tbl.tmp$BD_min)
tbl.tmp$BD_min = factor(x, levels=x[order(tbl.tmp$BD_min)])
# plotting relative abundances
ggplot(tbl.tmp, aes(BD_min, count, fill=taxon, group=taxon)) +
geom_bar(stat='identity', position='fill') +
labs(x='Buoyant density', y='Relative abundance') +
coord_flip() +
theme_bw() +
theme(
text = element_text(size=16),
axis.text.y = element_blank()
)
In [359]:
!cd $workDir; \
$SIPSimExe fragGC genomes/genomes10.txt \
--fp ./genomes/ \
--fr 515Fm-927Rm.fna \
--fld skewed-normal,9000,2500,-5 \
--flr 500,None \
--np 10 \
> genome10_ampFragGC_difSkew.txt
In [360]:
%%R -i workDir
setwd(workDir)
infile = 'genome10_ampFragGC_difSkew.txt'
tbl = read.delim(infile, sep='\t')
tbl$frag_len = abs(tbl$fragEnd - tbl$fragStart) + 1
In [361]:
%%R -h 350 -w 650
# stats
print(length(tbl[tbl$frag_len < 4000, 'frag_len']) / nrow(tbl))
# plot
ggplot(tbl, aes(frag_len)) +
geom_histogram(binwidth=100) +
geom_vline(xintercept=4000, linetype='dashed', alpha=0.5) +
theme(
text = element_text(size=18)
)
In [336]:
%%bash -s "$workDir" "$SIPSimExe"
cd $1
$2 OTU_table genome10_ampFragGC_difSkew.txt genome10_comm_n3.txt \
genome10_comm_n3_incorp.txt genome10_comm_n3_fracs.txt \
--abs_abund 1e8 > genome10_OTU_abnd1e8.txt
In [342]:
%%R -i workDir
# loading file
infile = paste(c(workDir, 'genome10_OTU_abnd1e8.txt'), collapse='/')
tbl = read.csv(infile, sep='\t')
In [343]:
%%R
# formatting table
tbl$BD_min = gsub('-.+', '', tbl$fractions)
tbl$BD_min = as.numeric(tbl$BD_min)
tbl$BD_max = gsub('.+-', '', tbl$fractions)
tbl$BD_max = as.numeric(tbl$BD_max)
In [344]:
%%R
# summarizing counts (should be approx. total abundance)
tbl %>%
group_by(library) %>%
summarize(sum(count))
In [345]:
%%R -w 800
# plotting absolute abundances
## BD for G+C of 0 or 100
BD.GCp0 = 0 * 0.098 + 1.66
BD.GCp100 = 1 * 0.098 + 1.66
## plot
ggplot(tbl, aes(BD_min, count, fill=taxon, group=taxon)) +
geom_area(stat='identity', alpha=0.5, position='dodge') +
geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
facet_grid(library ~ .) +
labs(x='Buoyant density') +
theme( text = element_text(size=16) )
In [346]:
%%R -w 800
# plotting relative abundances
ggplot(tbl, aes(BD_min, count, fill=taxon, group=taxon)) +
geom_area(stat='identity', alpha=0.8, position='fill') +
geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
facet_grid(library ~ .) +
labs(x='Buoyant density') +
theme( text = element_text(size=16) )
In [ ]:
In [ ]:
In [367]:
%%R -i workDir
infile = paste(c(workDir, 'genome10_OTU_abnd1e6.txt'), collapse='/')
tbl = read.csv(infile, sep='\t', row.names=1)
tbl$taxon_name = rownames(tbl)
In [368]:
%%R
tbl.m = melt(tbl, id.var=c('taxon_name'))
colnames(tbl.m) = c('taxon_name', 'variable', 'abundance')
tbl.m$lib = gsub('X|\\..+', '', tbl.m$variable)
tbl.m$BD_min = gsub('X[0-9]+\\.([0-9]+\\.[0-9]+).+', '\\1', tbl.m$variable)
tbl.m$BD_min = as.numeric(tbl.m$BD_min)
tbl.m$BD_max = gsub('X[0-9]+\\.[0-9]+\\.[0-9]+\\.(.+)', '\\1', tbl.m$variable)
tbl.m$BD_max = as.numeric(tbl.m$BD_max)
In [369]:
%%R -w 800
ggplot(tbl.m, aes(BD_min, abundance, color=taxon_name, group=taxon_name)) +
geom_point(size=1.5) +
geom_line(alpha=0.5) +
facet_grid(lib ~ .) +
theme( text = element_text(size=16) )
In [370]:
%%R -w 800
ggplot(tbl.m, aes(BD_min, abundance, fill=taxon_name, group=taxon_name)) +
geom_area(stat='identity', alpha=0.5, position='dodge') +
facet_grid(lib ~ .) +
labs(x='Buoyant density') +
theme( text = element_text(size=16) )
In [373]:
%%R -w 800 -h 800
tbl.m2 = tbl.m
tbl.m2$taxon_name = gsub("_","\n", tbl.m2$taxon_name)
ggplot(tbl.m2, aes(BD_min, abundance, fill=taxon_name, group=taxon_name)) +
geom_area(stat='identity', alpha=0.5, position='dodge') +
facet_grid(taxon_name ~ lib) <