Running SIPSim pipeline to simulate priming_exp gradient dataset

  • Basing simulation params off of priming_exp dataset
    • Using 'defaults' for other parameters

Setting variables


In [29]:
workDir = '/home/nick/notebook/SIPSim/dev/priming_exp/validation_sample/X12C.700.14.05_default/'
genomeDir = '/home/nick/notebook/SIPSim/dev/priming_exp/genomes/'
allAmpFrags = '/home/nick/notebook/SIPSim/dev/bac_genome1210/validation/ampFrags.pkl'
otuTableFile = '/var/seq_data/priming_exp/data/otu_table.txt'
metaDataFile = '/var/seq_data/priming_exp/data/allsample_metadata_nomock.txt'
primerFile = '/home/nick/notebook/SIPSim/dev/515F-806R.fna'

cdhit_dir = '/home/nick/notebook/SIPSim/dev/priming_exp/CD-HIT/'
R_dir = '/home/nick/notebook/SIPSim/lib/R/'
figureDir = '/home/nick/notebook/SIPSim/figures/'


# simulation params
comm_richness = 2170
seq_per_fraction = ['lognormal', 10.096, 1.116]

# for making genome_map file for genome fragment simulation
taxonMapFile = os.path.join(cdhit_dir, 'target_taxa.txt')
genomeFilterFile = os.path.join(cdhit_dir, 'genomeFile_seqID_filt.txt')
abundFile = os.path.join(workDir, 'X12C.700.14.05.NA_OTU.txt')

# misc
nprocs = 20

Init


In [7]:
import glob
import cPickle as pickle
import copy
from IPython.display import Image

In [8]:
%load_ext rpy2.ipython


The rpy2.ipython extension is already loaded. To reload it, use:
  %reload_ext rpy2.ipython

In [9]:
%%R
library(ggplot2)
library(dplyr)
library(tidyr)
library(gridExtra)

In [10]:
if not os.path.isdir(workDir):
    os.makedirs(workDir)

Creating a community file from OTU table


In [11]:
%%R -i otuTableFile 

tbl.otu = read.delim(otuTableFile, sep='\t')
tbl.otu[1:3,1:5]


     OTUId X12C.700.45.01.24 X12C.700.14.06.14 X12C.100.14.05.18
1 OTU.4776               157                34                13
2 OTU.2864                 0                 1                 2
3 OTU.8170                 0                 0                 0
  X12C.700.14.06.05
1               142
2                16
3                 0

In [12]:
%%R
tbl.otu.w = tbl.otu %>%
    gather('sample', 'count', 2:ncol(tbl.otu)) %>%
    group_by(sample) %>%
    mutate(rel_abund_perc = count / sum(count) * 100) %>%
    ungroup %>%
    filter(count > 0, sample == 'X12C.700.14.05.NA') %>%
    mutate(library = '1', rank = row_number(-rel_abund_perc)) %>%
    select(library, 'taxon_name' = OTUId, rel_abund_perc, rank) %>%
    arrange(rank)
tbl.otu.w %>% head(n=5)


Source: local data frame [5 x 4]

  library taxon_name rel_abund_perc rank
1       1     OTU.22       2.431002    1
2       1      OTU.1       2.202202    2
3       1     OTU.20       1.892369    3
4       1      OTU.2       1.682635    4
5       1    OTU.108       1.577768    5

In [14]:
%%R -i comm_richness
# number of OTUs
n.OTUs = tbl.otu.w$taxon_name %>% unique %>% length
cat('Number of OTUs: ', n.OTUs, '\n')

# assertion
cat('Community richness = number of OTUs?  ', comm_richness == n.OTUs, '\n')


Number of OTUs:  2170 
Community richness = number of OTUs?   TRUE 

In [15]:
%%R -i workDir

commFile = paste(c(workDir, 'comm.txt'), collapse='/')
write.table(tbl.otu.w, commFile, sep='\t', quote=F, row.names=F)

Plotting community distribution


In [16]:
%%R -i workDir

commFile = paste(c(workDir, 'comm.txt'), collapse='/')
comm = read.delim(commFile, sep='\t')
comm %>% head


  library taxon_name rel_abund_perc rank
1       1     OTU.22       2.431002    1
2       1      OTU.1       2.202202    2
3       1     OTU.20       1.892369    3
4       1      OTU.2       1.682635    4
5       1    OTU.108       1.577768    5
6       1      OTU.8       1.201201    6

In [17]:
%%R -w 900 -h 350

ggplot(comm, aes(rank, rel_abund_perc)) +
    geom_point() +
    labs(x='Rank', y='% relative abundance', title='Priming experiment community abundance distribution') +
    theme_bw() +
    theme(
        text = element_text(size=16)
        )


Simulating fragments

Making a genome index file to map genome fasta files to OTUs

  • Will be used for community simulation
  • Just OTUs with association to genomes

In [19]:
%%R -i taxonMapFile -i genomeFilterFile 

taxonMap = read.delim(taxonMapFile, sep='\t') %>%
    select(target_genome, OTU) %>%
    distinct()
taxonMap %>% nrow %>% print
taxonMap %>% head(n=3) %>% print

breaker = '----------------\n'
cat(breaker)

genomeFilter = read.delim(genomeFilterFile, sep='\t', header=F) 
genomeFilter %>% nrow %>% print
genomeFilter %>% head(n=3) %>% print

cat(breaker)

comm = read.delim(commFile, sep='\t') 
comm %>% nrow %>% print
comm %>% head(n=3) %>% print


[1] 236
                               target_genome      OTU
1 CP001738_Thermomonospora_curvata_DSM_43183 OTU.8540
2 CP001738_Thermomonospora_curvata_DSM_43183 OTU.9267
3 CP001738_Thermomonospora_curvata_DSM_43183 OTU.1457
----------------
[1] 187
                                                         V1
1                  CP003093_Pseudoxanthomonas_spadix_BD_a59
2                  CP000511_Mycobacterium_vanbaalenii_PYR_1
3 CP003344_Desulfitobacterium_dichloroeliminans_LMG_P_21439
                                                                                                            V2
1                  /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Pseudoxanthomonas_spadix_BD-a59.fasta
2                  /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Mycobacterium_vanbaalenii_PYR-1.fasta
3 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Desulfitobacterium_dichloroeliminans_LMG_P-21439.fasta
----------------
[1] 2170
  library taxon_name rel_abund_perc rank
1       1     OTU.22       2.431002    1
2       1      OTU.1       2.202202    2
3       1     OTU.20       1.892369    3

In [20]:
%%R
taxonMap$OTU %>% table %>% sort(decreasing=T) %>% head


.
    OTU.1    OTU.10   OTU.101   OTU.102 OTU.10237  OTU.1035 
        1         1         1         1         1         1 

In [21]:
%%R

tbl.j = inner_join(taxonMap, genomeFilter, c('target_genome' = 'V1')) %>%
     rename('fasta_file' = V2) %>%
     select(OTU, fasta_file, target_genome)

tbl.j %>% head(n=3)


       OTU
1 OTU.8540
2 OTU.9267
3 OTU.1457
                                                                                     fasta_file
1 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Thermomonospora_curvata_DSM_43183.fasta
2 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Thermomonospora_curvata_DSM_43183.fasta
3 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Thermomonospora_curvata_DSM_43183.fasta
                               target_genome
1 CP001738_Thermomonospora_curvata_DSM_43183
2 CP001738_Thermomonospora_curvata_DSM_43183
3 CP001738_Thermomonospora_curvata_DSM_43183

In [22]:
%%R
tbl.j$OTU %>% table %>% sort(decreasing=T) %>% head


.
    OTU.1    OTU.10   OTU.101   OTU.102 OTU.10237  OTU.1035 
        1         1         1         1         1         1 

In [25]:
%%R
tbl.j2 = inner_join(tbl.j, comm, c('OTU' = 'taxon_name')) 

n.target.genomes = tbl.j2$OTU %>% unique %>% length
cat('Number of target OTUs: ', n.target.genomes, '\n')
cat('--------', '\n')
tbl.j2 %>% head(n=3)


Number of target OTUs:  108 
-------- 
       OTU
1 OTU.1457
2   OTU.77
3 OTU.2347
                                                                                        fasta_file
1    /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Thermomonospora_curvata_DSM_43183.fasta
2      /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Pseudoxanthomonas_spadix_BD-a59.fasta
3 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Micromonospora_aurantiaca_ATCC_27029.fasta
                                  target_genome library rel_abund_perc rank
1    CP001738_Thermomonospora_curvata_DSM_43183       1    0.009533343 1038
2      CP003093_Pseudoxanthomonas_spadix_BD_a59       1    0.019066686  700
3 CP002162_Micromonospora_aurantiaca_ATCC_27029       1    0.004766671 1380

In [26]:
%%R -i workDir

outFile = paste(c(workDir, 'target_genome_index.txt'), collapse='/')
write.table(tbl.j2, outFile, sep='\t', quote=F, row.names=F, col.names=F)

Plotting community abundance distribution of target genomes


In [27]:
%%R -w 900 -h 350

ggplot(tbl.j2, aes(rank, rel_abund_perc)) +
    geom_point(size=3, shape='O', color='red') +
    labs(x='Rank', y='% relative abundance', title='Priming experiment community abundance distribution') +
    theme_bw() +
    theme(
        text = element_text(size=16)
        )


Simulating fragments of genomes that match priming_exp bulk OTUs


In [30]:
!cd $workDir; \
    SIPSim fragments \
    target_genome_index.txt \
    --fp $genomeDir \
    --fr $primerFile \
    --fld skewed-normal,9000,2500,-5 \
    --flr None,None \
    --nf 10000 \
    --np $nprocs \
    2> ampFrags.log \
    > ampFrags.pkl

Appending fragments from randomly selected genomes of total dataset (n=1210)

  • This is to obtain the richness of the bulk soil community
  • Random OTUs will be named after non-target OTUs

Making list of non-target OTUs


In [284]:
%%R -i workDir
# loading files

## target genome index (just OTUs with associated genome)
inFile = paste(c(workDir, 'target_genome_index.txt'), collapse='/')
tbl.target = read.delim(inFile, sep='\t', header=F)
colnames(tbl.target) = c('OTUId', 'fasta_file', 'genome_name')

## comm file of total community OTUs 
commFile = paste(c(workDir, 'comm.txt'), collapse='/')
tbl.comm = read.delim(commFile, sep='\t')

In [285]:
%%R
# just OTUs w/out an associated genome
tbl.j = anti_join(tbl.comm, tbl.target, c('taxon_name' = 'OTUId'))
n.nontarget.genomes = tbl.j$taxon_name %>% length
message('Number of non-target genomes: ', n.nontarget.genomes)
message('---------')
tbl.j %>% head(n=5)


Number of non-target genomes: 2062
---------
  library taxon_name rel_abund_perc rank
1       1  OTU.11779    0.004766671 2170
2       1   OTU.9498    0.004766671 2169
3       1   OTU.3673    0.004766671 2168
4       1    OTU.493    0.004766671 2167
5       1   OTU.2500    0.004766671 2166

In [286]:
%%R -i comm_richness
# checking assumptions
message('Target + nonTarget richness = total community richness?: ',
         n.target.genomes + n.nontarget.genomes == comm_richness)


Target + nonTarget richness = total community richness?: TRUE

In [287]:
%%R -i workDir
# writing out non-target OTU file
outFile = paste(c(workDir, 'comm_nonTarget.txt'), collapse='/')
write.table(tbl.j, outFile, sep='\t', quote=F, row.names=F)

Randomly selecting amplicon fragment length-GC KDEs from total genome pool


In [288]:
# List of non-target OTUs
inFile = os.path.join(workDir, 'comm_nonTarget.txt')
nonTarget = pd.read_csv(inFile, sep='\t')['taxon_name'].tolist()

print 'Number of non-target OTUs: {}'.format(len(nonTarget))
nonTarget[:4]


Number of non-target OTUs: 2062
Out[288]:
['OTU.11779', 'OTU.9498', 'OTU.3673', 'OTU.493']

In [289]:
# loading amplicon fragments from full genome KDE dataset
inFile = os.path.join(workDir, 'ampFrags.pkl')
ampFrag_target = []
with open(inFile, 'rb') as iFH:
    ampFrag_target = pickle.load(iFH)
print 'Target OTU richness: {}'.format(len(ampFrag_target))


Target OTU richness: 108

In [290]:
# loading amplicon fragments from full genome KDE dataset
ampFrag_all = []
with open(allAmpFrags, 'rb') as iFH:
    ampFrag_all = pickle.load(iFH)
print 'Count of frag-GC KDEs for all genomes: {}'.format(len(ampFrag_all))


Count of frag-GC KDEs for all genomes: 1210

In [291]:
# random selection from list
#target_richness = len(ampFrag_target)

target_richness = len(ampFrag_target)
richness_needed = comm_richness - target_richness
print 'Number of random taxa needed to reach richness: {}'.format(richness_needed)

if richness_needed > 0:
    index = range(target_richness)
    index = np.random.choice(index, richness_needed)
    
    ampFrag_rand = []
    for i in index:
        sys.stderr.write('{},'.format(i))
        ampFrag_rand.append(copy.deepcopy(ampFrag_all[i]))
else:
    ampFrag_rand = []


23,52,49,82,41,1,61,34,57,64,99,80,61,105,97,75,60,68,78,32,58,87,37,92,32,107,74,65,10,71,98,56,69,67,70,5,37,11,13,12,36,35,8,4,104,78,3,0,86,34,36,12,53,51,50,71,36,71,75,27,42,35,50,17,77,95,46,30,59,68,67,102,78,103,51,19,24,90,50,91,66,100,63,3,17,18,92,96,88,67,54,64,31,25,26,8,77,39,86,83,74,107,85,31,10,83,63,65,26,88,107,106,46,88,72,40,15,57,47,10,39,53,55,89,57,35,72,40,5,90,100,65,82,75,13,68,90,7,51,15,84,102,1,89,53,21,69,14,10,45,37,101,63,103,95,93,94,21,65,66,16,15,37,78,21,12,11,73,50,21,79,95,94,94,77,74,1,24,56,84,75,9,20,69,41,25,38,96,49,23,38,91,5,58,80,3,32,64,40,100,86,9,58,42,69,100,17,61,59,99,84,4,89,84,33,41,75,0,98,25,23,39,5,107,74,64,16,9,28,79,95,17,90,35,53,68,25,11,30,85,1,75,97,102,0,41,11,79,46,1,40,74,26,7,65,62,79,80,79,41,38,57,62,83,77,4,62,57,25,69,103,61,106,86,65,36,92,33,73,54,64,76,59,32,8,95,75,23,23,49,66,57,101,92,74,36,26,95,106,80,27,72,61,85,74,40,100,92,62,85,46,80,21,53,88,20,41,35,3,36,32,102,100,45,28,71,42,75,60,102,48,40,105,9,80,107,21,21,37,31,45,40,5,62,70,11,94,71,89,50,12,33,90,26,2,100,78,64,16,56,0,34,5,98,84,71,73,21,65,45,2,7,21,92,38,103,47,102,101,61,86,36,60,25,81,2,102,0,35,41,12,42,78,95,21,46,83,0,3,13,15,59,55,36,95,73,107,63,47,79,71,75,93,82,51,49,90,15,5,70,1,22,0,35,9,74,101,83,19,107,47,83,38,100,53,29,21,27,74,85,107,101,80,95,8,60,17,73,105,80,82,10,80,90,22,64,45,45,11,1,16,93,99,102,46,6,49,58,17,50,75,89,107,0,71,44,34,15,82,18,17,39,34,44,107,27,31,49,84,80,13,89,59,67,2,3,107,57,62,101,39,74,1,70,23,27,10,7,2,82,33,99,3,50,92,98,78,21,3,59,62,30,81,38,30,44,105,4,98,56,92,90,49,52,91,2,21,106,105,93,91,47,92,78,91,23,0,88,83,68,38,31,52,18,76,59,16,67,60,56,75,43,52,14,74,68,102,44,33,97,98,58,54,104,87,76,35,32,25,70,21,75,62,99,50,73,5,78,61,48,86,74,44,80,72,4,41,86,105,32,103,79,101,44,30,27,31,30,95,54,26,19,19,66,45,41,10,0,84,36,73,15,21,70,89,95,87,91,103,29,59,28,76,57,2,22,53,59,54,82,75,8,20,89,3,30,34,33,107,69,55,100,2,51,30,102,89,60,89,107,78,92,3,65,106,6,96,23,30,59,20,72,107,3,70,14,103,62,24,26,2,96,101,73,28,55,87,78,89,42,79,32,33,37,56,17,59,25,0,32,103,105,35,46,7,20,0,106,105,33,10,107,3,49,88,32,13,58,12,67,45,96,30,54,12,17,27,36,60,86,10,30,65,71,17,24,0,1,49,45,3,14,45,47,34,97,65,62,40,96,41,75,25,53,25,104,102,100,89,0,10,57,72,58,79,76,103,60,85,65,92,34,50,92,85,70,62,67,52,71,91,67,97,8,61,2,50,86,30,38,24,64,17,10,57,0,13,68,46,92,65,88,66,44,63,46,20,67,49,11,77,45,75,16,98,52,78,20,0,48,73,71,39,77,100,55,66,49,27,5,55,62,24,4,16,81,88,104,46,27,66,13,23,74,81,62,107,106,14,75,48,72,31,59,56,63,49,0,69,105,20,48,37,61,36,72,15,31,10,33,57,105,81,31,32,107,68,61,76,84,41,20,105,13,41,45,103,66,73,105,79,100,79,75,17,74,45,69,2,88,20,10,100,23,75,65,82,53,17,98,65,14,64,73,3,74,5,85,35,68,28,64,31,20,73,37,89,8,37,39,73,102,91,16,99,42,40,52,91,35,88,82,19,102,55,54,19,96,36,24,2,2,65,33,78,12,77,102,47,57,26,48,99,21,94,96,46,104,40,67,53,78,101,4,3,79,79,26,62,81,98,17,2,100,91,102,24,101,61,101,107,34,65,55,38,26,101,1,73,84,82,80,60,96,62,92,103,23,52,79,49,74,7,101,8,10,81,82,101,14,59,87,19,60,100,82,44,63,79,90,83,93,99,24,62,96,9,26,1,38,67,24,75,75,36,98,9,82,30,88,21,106,53,92,103,101,93,68,50,97,107,15,4,74,1,81,51,16,57,57,84,98,3,3,78,33,43,47,93,80,62,9,48,9,32,40,37,53,42,58,57,82,15,95,9,63,50,95,75,18,54,36,64,71,50,97,82,17,30,0,23,19,24,41,88,76,18,61,24,97,11,56,68,80,3,72,91,85,32,93,91,43,7,24,103,74,106,10,83,17,74,24,98,46,41,34,2,89,28,75,61,19,22,67,14,77,87,86,101,79,61,83,58,78,13,31,77,1,36,3,37,43,24,36,51,82,20,31,47,47,1,18,64,79,94,74,59,31,48,28,57,107,40,93,105,74,87,27,52,6,96,71,90,11,85,11,2,41,87,25,43,50,33,98,66,71,101,12,60,26,53,46,33,58,23,97,82,87,99,80,25,63,61,77,59,64,50,51,77,81,57,0,32,42,50,72,39,73,91,104,106,79,35,46,33,64,90,83,75,17,20,44,94,66,35,35,95,12,45,54,69,59,27,20,83,65,78,98,58,46,36,76,83,45,84,14,96,17,0,67,68,34,107,78,84,0,54,83,98,21,13,64,55,53,90,96,79,107,31,82,61,73,45,46,85,92,41,100,52,96,93,104,21,2,58,3,32,29,89,40,0,96,82,31,88,12,107,87,14,82,9,18,11,25,93,72,73,8,91,102,40,60,50,39,101,13,70,95,34,74,94,51,65,4,93,39,77,66,25,41,7,35,79,86,27,11,58,68,84,2,6,25,7,84,42,2,80,78,84,102,80,102,98,80,1,90,61,102,52,29,84,74,60,23,98,75,84,37,102,80,72,86,99,30,69,22,91,62,4,73,88,50,3,75,34,71,27,19,52,52,76,61,31,96,103,27,51,30,36,33,53,77,26,95,13,41,78,41,96,13,21,73,55,86,11,46,53,20,7,0,81,87,62,92,100,20,26,53,25,26,85,43,20,107,53,60,31,62,65,38,44,77,22,97,73,26,39,11,3,32,72,90,0,26,29,45,64,70,44,51,36,52,30,93,12,15,22,78,21,6,55,28,52,102,68,72,26,48,35,10,80,101,46,10,9,69,99,12,68,67,30,43,18,38,70,40,22,60,12,62,77,63,44,15,3,41,11,17,0,7,26,107,61,12,76,42,74,92,74,70,67,38,45,104,107,88,65,31,74,35,65,17,29,6,72,52,74,52,65,14,76,43,79,34,91,97,103,41,62,97,60,104,68,46,71,69,105,93,56,60,80,72,103,54,44,94,59,97,18,0,5,46,33,103,83,106,51,101,84,79,4,94,73,73,48,101,14,25,49,23,25,9,91,69,11,21,29,23,81,24,67,74,15,75,31,24,90,25,99,48,30,18,10,28,78,77,43,50,21,25,20,105,22,63,22,103,85,44,46,53,41,14,29,8,35,40,15,32,86,71,29,8,74,29,43,42,89,40,6,58,94,57,60,91,38,101,42,79,12,75,81,25,31,34,87,71,106,69,2,78,89,10,47,59,53,81,64,105,33,56,65,2,57,0,21,18,51,96,82,64,37,29,107,2,105,61,76,70,46,76,34,0,25,30,39,48,107,35,81,107,26,43,12,18,94,2,80,5,18,62,82,94,29,3,66,20,39,71,46,30,46,65,104,59,2,26,47,82,12,54,37,63,61,35,22,15,1,38,32,13,73,46,98,69,45,15,57,97,9,49,0,62,76,5,47,18,91,48,32,43,51,10,22,64,48,24,42,74,76,70,7,90,47,48,4,74,2,50,31,44,9,77,27,100,28,11,10,76,88,20,22,27,48,13,60,89,66,70,23,88,79,47,8,61,4,19,54,62,106,95,27,107,48,73,93,91,73,101,101,97,7,54,34,24,51,59,15,57,28,77,7,93,45,26,59,99,74,65,21,44,97,4,57,52,44,2,42,48,88,15,22,5,35,23,62,72,84,87,0,104,86,72,94,105,21,57,27,93,55,60,74,31,68,104,28,81,105,12,81,38,52,76,77,11,106,26,6,82,52,27,93,25,44,1,13,74,66,97,98,87,76,5,37,50,32,82,33,50,64,106,84,28,71,80,96,2,70,97,101,59,105,41,45,99,103,28,44,2,72,83,21,74,9,0,89,41,45,63,50,45,85,8,13,88,100,34,44,56,33,39,42,94,38,76,103,67,15,84,37,98,91,39,95,35,30,27,35,93,7,54,94,85,34,50,25,86,85,73,82,26,105,103,53,40,85,75,9,8,9,41,66,81,83,38,51,76,74,78,9,22,63,89,42,27,47,99,106,28,16,85,44,88,
Number of random taxa needed to reach richness: 2062

In [292]:
# renaming randomly selected KDEs by non-target OTU-ID
for i in range(len(ampFrag_rand)):
    ampFrag_rand[i][0] = nonTarget[i]

In [293]:
# appending random taxa to target taxa and writing
outFile = os.path.join(workDir, 'ampFrags_wRand.pkl')

with open(outFile, 'wb') as oFH:
    x = ampFrag_target + ampFrag_rand
    print 'Number of taxa in output: {}'.format(len(x))
    pickle.dump(x, oFH)


Number of taxa in output: 2170

Converting fragments to kde object


In [ ]:
!cd $workDir; \
    SIPSim fragment_kde \
    ampFrags_wRand.pkl \
    > ampFrags_wRand_kde.pkl

Adding diffusion


In [ ]:
!cd $workDir; \
    SIPSim diffusion \
    ampFrags_wRand_kde.pkl \
    --np $nprocs \
    > ampFrags_wRand_kde_dif.pkl


Processing: OTU.5687
Processing: OTU.12484
Processing: OTU.1376
Processing: OTU.3442
Processing: OTU.3023
Processing: OTU.2015
Processing: OTU.247
Processing: OTU.148
Processing: OTU.8432
Processing: OTU.8949
Processing: OTU.197
Processing: OTU.758
Processing: OTU.2304
Processing: OTU.6164
Processing: OTU.5875
Processing: OTU.774
Processing: OTU.7536
Processing: OTU.4541
Processing: OTU.3020
Processing: OTU.6855
Processing: OTU.3152
Processing: OTU.13746
Processing: OTU.243
Processing: OTU.4901
Processing: OTU.4907
Processing: OTU.1624
Processing: OTU.8941
Processing: OTU.126
Processing: OTU.508
Processing: OTU.1087
Processing: OTU.5241
Processing: OTU.6994
Processing: OTU.4429
Processing: OTU.9497
Processing: OTU.6875
Processing: OTU.196
Processing: OTU.11076
Processing: OTU.5683
Processing: OTU.7221
Processing: OTU.7227
Processing: OTU.179
Processing: OTU.394
Processing: OTU.777
Processing: OTU.1155
Processing: OTU.12912
Processing: OTU.164
Processing: OTU.8946
Processing: OTU.9334
Processing: OTU.509
Processing: OTU.1086
Processing: OTU.6077
Processing: OTU.6658
Processing: OTU.2969
Processing: OTU.9324
Processing: OTU.8100
Processing: OTU.9314
Processing: OTU.1793
Processing: OTU.630
Processing: OTU.631
Processing: OTU.6300
Processing: OTU.1154
Processing: OTU.191
Processing: OTU.12094
Processing: OTU.5466
Processing: OTU.165
Processing: OTU.3307
Processing: OTU.1941
Processing: OTU.771
Processing: OTU.2190
Processing: OTU.1085
Processing: OTU.1908
Processing: OTU.6992
Processing: OTU.1141
Processing: OTU.1438
Processing: OTU.6871
Processing: OTU.5328
Processing: OTU.1792
Processing: OTU.633
Processing: OTU.6074
Processing: OTU.99
Processing: OTU.6303
Processing: OTU.7315
Processing: OTU.242
Processing: OTU.166
Processing: OTU.3308
Processing: OTU.2192
Processing: OTU.190
Processing: OTU.1083
Processing: OTU.1082
Processing: OTU.752
Processing: OTU.995
Processing: OTU.4152
Processing: OTU.1147
Processing: OTU.9454
Processing: OTU.82
Processing: OTU.770
Processing: OTU.1790
Processing: OTU.12917
Processing: OTU.634
Processing: OTU.7314
Processing: OTU.385
Processing: OTU.1150
Processing: OTU.6073
Processing: OTU.167
Processing: OTU.98
Processing: OTU.1524
Processing: OTU.2194
Processing: OTU.1081
Processing: OTU.265
Processing: OTU.9171
Processing: OTU.9498
Processing: OTU.83
Processing: OTU.193
Processing: OTU.192
Processing: OTU.1797
Processing: OTU.751
Processing: OTU.997
Processing: OTU.635
Processing: OTU.13304
Processing: OTU.1949
Processing: OTU.4432
Processing: OTU.9159
Processing: OTU.160
Processing: OTU.3231
Processing: OTU.8842
Processing: OTU.2195
Processing: OTU.1080
Processing: OTU.3002
Processing: OTU.599
Processing: OTU.1145
Processing: OTU.415
Processing: OTU.342
Processing: OTU.80
Processing: OTU.1796
Processing: OTU.636
Processing: OTU.12726
Processing: OTU.7823
Processing: OTU.1638
Processing: OTU.750
Processing: OTU.996
Processing: OTU.161
Processing: OTU.3236
Processing: OTU.779
Processing: OTU.2197
Processing: OTU.3810
Processing: OTU.3817
Processing: OTU.9158
Processing: OTU.1322
Processing: OTU.6496
Processing: OTU.8695
Processing: OTU.81
Processing: OTU.8846
Processing: OTU.1799
Processing: OTU.3799
Processing: OTU.637
Processing: OTU.91
Processing: OTU.2178
Processing: OTU.7827
Processing: OTU.162
Processing: OTU.8051
Processing: OTU.2200
Processing: OTU.609
Processing: OTU.5635
Processing: OTU.756
Processing: OTU.991
Processing: OTU.5338
Processing: OTU.10634
Processing: OTU.2987
Processing: OTU.6823
Processing: OTU.86
Processing: OTU.778
Processing: OTU.7040
Processing: OTU.8048
Processing: OTU.638
Processing: OTU.2176
Processing: OTU.2949
Processing: OTU.3798
Processing: OTU.163
Processing: OTU.90
Processing: OTU.1528
Processing: OTU.1999
Processing: OTU.6506
Processing: OTU.1321
Processing: OTU.4206
Processing: OTU.10638
Processing: OTU.13967
Processing: OTU.87
Processing: OTU.4683
Processing: OTU.359
Processing: OTU.4127
Processing: OTU.10804
Processing: OTU.755
Processing: OTU.993
Processing: OTU.639
Processing: OTU.1098
Processing: OTU.12238
Processing: OTU.2947
Processing: OTU.2370
Processing: OTU.519
Processing: OTU.2204
Processing: OTU.13614
Processing: OTU.7249
Processing: OTU.4204
Processing: OTU.1205
Processing: OTU.93
Processing: OTU.1509
Processing: OTU.1503
Processing: OTU.84
Processing: OTU.4680
Processing: OTU.4126
Processing: OTU.1555
Processing: OTU.643
Processing: OTU.12125
Processing: OTU.358
Processing: OTU.2945
Processing: OTU.754
Processing: OTU.992
Processing: OTU.2372
Processing: OTU.2542
Processing: OTU.1099
Processing: OTU.2205
Processing: OTU.12219
Processing: OTU.4209
Processing: OTU.7831
Processing: OTU.1500
Processing: OTU.85
Processing: OTU.2223
Processing: OTU.4129
Processing: OTU.8111
Processing: OTU.1556
Processing: OTU.3790
Processing: OTU.748
Processing: OTU.92
Processing: OTU.1976
Processing: OTU.2944
Processing: OTU.982
Processing: OTU.2541
Processing: OTU.2759
Processing: OTU.5529
Processing: OTU.357
Processing: OTU.12748
Processing: OTU.5440
Processing: OTU.2499
Processing: OTU.3165
Processing: OTU.1745
Processing: OTU.7834
Processing: OTU.1501
Processing: OTU.88
Processing: OTU.1095
Processing: OTU.5898
Processing: OTU.7816
Processing: OTU.7672
Processing: OTU.11380
Processing: OTU.759
Processing: OTU.1975
Processing: OTU.5262
Processing: OTU.4328
Processing: OTU.188
Processing: OTU.3173
Processing: OTU.95
Processing: OTU.144
Processing: OTU.6068
Processing: OTU.2753
Processing: OTU.7111
Processing: OTU.231
Processing: OTU.1200
Processing: OTU.1504
Processing: OTU.11478
Processing: OTU.8763
Processing: OTU.2971
Processing: OTU.604
Processing: OTU.8115
Processing: OTU.1558
Processing: OTU.2246
Processing: OTU.3499
Processing: OTU.3477
Processing: OTU.1097
Processing: OTU.9
Processing: OTU.554
Processing: OTU.981
Processing: OTU.6069
Processing: OTU.6217
Processing: OTU.13619
Processing: OTU.189
Processing: OTU.619
Processing: OTU.1203
Processing: OTU.97
Processing: OTU.8070
Processing: OTU.8004
Processing: OTU.131
Processing: OTU.8934
Processing: OTU.2553
Processing: OTU.740
Processing: OTU.8
Processing: OTU.607
Processing: OTU.9549
Processing: OTU.3497
Processing: OTU.1715
Processing: OTU.3175
Processing: OTU.511
Processing: OTU.1090
Processing: OTU.5611
Processing: OTU.6447
Processing: OTU.4256
Processing: OTU.9183
Processing: OTU.7837
Processing: OTU.89
Processing: OTU.7650
Processing: OTU.3707
Processing: OTU.8930
Processing: OTU.2006
Processing: OTU.182
Processing: OTU.741
Processing: OTU.96
Processing: OTU.1
Processing: OTU.14219
Processing: OTU.984
Processing: OTU.510
Processing: OTU.1633
Processing: OTU.6213
Processing: OTU.6443
Processing: OTU.3490
Processing: OTU.4524
Processing: OTU.8583
Processing: OTU.2956
Processing: OTU.733
Processing: OTU.4973
Processing: OTU.8427
Processing: OTU.4354
Processing: OTU.12047
Processing: OTU.529
Processing: OTU.742
Processing: OTU.3
Processing: OTU.11661
Processing: OTU.183
Processing: OTU.1315
Processing: OTU.3176
Processing: OTU.513
Processing: OTU.762
Processing: OTU.5742
Processing: OTU.345
Processing: OTU.2950
Processing: OTU.168
Processing: OTU.8424
Processing: OTU.4351
Processing: OTU.353
Processing: OTU.7322
Processing: OTU.298
Processing: OTU.3491
Processing: OTU.2143
Processing: OTU.743
Processing: OTU.744
Processing: OTU.1092
Processing: OTU.2
Processing: OTU.4818
Processing: OTU.988
Processing: OTU.2548
Processing: OTU.763
Processing: OTU.235
Processing: OTU.180
Processing: OTU.7661
Processing: OTU.6781
Processing: OTU.145
Processing: OTU.9408
Processing: OTU.9555
Processing: OTU.6049
Processing: OTU.785
Processing: OTU.2005
Processing: OTU.746
Processing: OTU.4
Processing: OTU.600
Processing: OTU.1308
Processing: OTU.2788
Processing: OTU.5797
Processing: OTU.2140
Processing: OTU.4663
Processing: OTU.3016
Processing: OTU.1093
Processing: OTU.760
Processing: OTU.7809
Processing: OTU.4117
Processing: OTU.6780
Processing: OTU.181
Processing: OTU.9954
Processing: OTU.9082
Processing: OTU.7275
Processing: OTU.296
Processing: OTU.747
Processing: OTU.127
Processing: OTU.5374
Processing: OTU.7
Processing: OTU.6853
Processing: OTU.2648
Processing: OTU.514
Processing: OTU.3015
Processing: OTU.603
Processing: OTU.1956
Processing: OTU.7808
Processing: OTU.12647
Processing: OTU.5794
Processing: OTU.1406
Processing: OTU.1435
Processing: OTU.396
Processing: OTU.6783
Processing: OTU.3827
Processing: OTU.186
Processing: OTU.9406
Processing: OTU.13773
Processing: OTU.378
Processing: OTU.297
Processing: OTU.2317
Processing: OTU.5090
Processing: OTU.1305
Processing: OTU.410
Processing: OTU.1537
Processing: OTU.506
Processing: OTU.507
Processing: OTU.7253
Processing: OTU.2688
Processing: OTU.11797
Processing: OTU.1433
Processing: OTU.10124
Processing: OTU.184
Processing: OTU.14205
Processing: OTU.3445
Processing: OTU.3368
Processing: OTU.1619
Processing: OTU.350
Processing: OTU.4921
Processing: OTU.2163
Processing: OTU.12397
Processing: OTU.2315
Processing: OTU.3823
Processing: OTU.5092
Processing: OTU.4227
Processing: OTU.504
Processing: OTU.505
Processing: OTU.199
Processing: OTU.767
Processing: OTU.596
Processing: OTU.1430
Processing: OTU.9211
Processing: OTU.10129
Processing: OTU.1536
Processing: OTU.185
Processing: OTU.1319
Processing: OTU.3332
Processing: OTU.371
Processing: OTU.831
Processing: OTU.647
Processing: OTU.3838
Processing: OTU.177
Processing: OTU.5865
Processing: OTU.2514
Processing: OTU.2320
Processing: OTU.502
Processing: OTU.734
Processing: OTU.198
Processing: OTU.764
Processing: OTU.7805
Processing: OTU.7669
Processing: OTU.4251
Processing: OTU.1628
Processing: OTU.137
Processing: OTU.3330
Processing: OTU.986
Processing: OTU.292
Processing: OTU.5470
Processing: OTU.411
Processing: OTU.1534
Processing: OTU.3837
Processing: OTU.5312
Processing: OTU.503
Processing: OTU.195
Processing: OTU.176
Processing: OTU.1952
Processing: OTU.393
Processing: OTU.1892
Processing: OTU.1969
Processing: OTU.340
Processing: OTU.1370
Processing: OTU.1372
Processing: OTU.5640
Processing: OTU.349
Processing: OTU.1314
Processing: OTU.3337
Processing: OTU.373
Processing: OTU.837
Processing: OTU.4783
Processing: OTU.5866
Processing: OTU.10401
Processing: OTU.1532
Processing: OTU.500
Processing: OTU.194
Processing: OTU.1612
Processing: OTU.7069
Processing: OTU.236
Processing: OTU.7788
Processing: OTU.344
Processing: OTU.6120
Processing: OTU.175
Processing: OTU.1313
Processing: OTU.372
Processing: OTU.290
Processing: OTU.1893
Processing: OTU.1965
Processing: OTU.570
Processing: OTU.5642
Processing: OTU.1762
Processing: OTU.4752
Processing: OTU.501
Processing: OTU.374
Processing: OTU.1089
Processing: OTU.140
Processing: OTU.1530
Processing: OTU.8067
Processing: OTU.1051
Processing: OTU.613
Processing: OTU.835
Processing: OTU.1311
Processing: OTU.576
Processing: OTU.174
Processing: OTU.1407
Processing: OTU.6514
Processing: OTU.1961
Processing: OTU.1420
Processing: OTU.1988
Processing: OTU.377
Processing: OTU.1088
Processing: OTU.3577
Processing: OTU.413
Processing: OTU.8065
Processing: OTU.1054
Processing: OTU.346
Processing: OTU.8204
Processing: OTU.574
Processing: OTU.173
Processing: OTU.2761
Processing: OTU.418
Processing: OTU.1144
Processing: OTU.1009
Processing: OTU.1981
Processing: OTU.566
Processing: OTU.1422
Processing: OTU.1611
Processing: OTU.6053
Processing: OTU.1413
Processing: OTU.3228
Processing: OTU.1055
Processing: OTU.1620
Processing: OTU.1994
Processing: OTU.12282
Processing: OTU.579
Processing: OTU.578
Processing: OTU.172
Processing: OTU.1898
Processing: OTU.419
Processing: OTU.13605
Processing: OTU.7866
Processing: OTU.707
Processing: OTU.12801
Processing: OTU.567
Processing: OTU.1425
Processing: OTU.151
Processing: OTU.1415
Processing: OTU.1538
Processing: OTU.3380
Processing: OTU.1056
Processing: OTU.4969
Processing: OTU.5320
Processing: OTU.294
Processing: OTU.171
Processing: OTU.2764
Processing: OTU.2354
Processing: OTU.1627
Processing: OTU.941
Processing: OTU.7102
Processing: OTU.317
Processing: OTU.293
Processing: OTU.10065
Processing: OTU.12754
Processing: OTU.568
Processing: OTU.10694
Processing: OTU.150
Processing: OTU.3841
Processing: OTU.1417
Processing: OTU.1057
Processing: OTU.12598
Processing: OTU.170
Processing: OTU.6201
Processing: OTU.149
Processing: OTU.11135
Processing: OTU.957
Processing: OTU.3771
Processing: OTU.291
Processing: OTU.7103
Processing: OTU.1186
Processing: OTU.569
Processing: OTU.1427
Processing: OTU.341
Processing: OTU.153
Processing: OTU.5612
Processing: OTU.597
Processing: OTU.316
Processing: OTU.2878
Processing: OTU.4138
Processing: OTU.2211
Processing: OTU.2300
Processing: OTU.3845
Processing: OTU.7770
Processing: OTU.956
Processing: OTU.146
Processing: OTU.1049
Processing: OTU.7104
Processing: OTU.1187
Processing: OTU.2510
Processing: OTU.1426
Processing: OTU.1458
Processing: OTU.7263
Processing: OTU.152
Processing: OTU.12311
Processing: OTU.1625
Processing: OTU.1334
Processing: OTU.6882
Processing: OTU.1419
Processing: OTU.315
Processing: OTU.5537
Processing: OTU.2302
Processing: OTU.775
Processing: OTU.1709
Processing: OTU.147
Processing: OTU.1048
Processing: OTU.910
Processing: OTU.1184
Processing: OTU.5666
Processing: OTU.1650
Processing: OTU.8462
Processing: OTU.155
Processing: OTU.848
Processing: OTU.1454
Processing: OTU.4212
Processing: OTU.11322
Processing: OTU.2216
Processing: OTU.559
Processing: OTU.416
Processing: OTU.2643
Processing: OTU.9054
Processing: OTU.314
Processing: OTU.1185
Processing: OTU.4398
Processing: OTU.7599
Processing: OTU.7612
Processing: OTU.154
Processing: OTU.4776
Processing: OTU.64
Processing: OTU.65
Processing: OTU.5346
Processing: OTU.4175
Processing: OTU.3779
Processing: OTU.295
Processing: OTU.558
Processing: OTU.417
Processing: OTU.2645
Processing: OTU.156
Processing: OTU.1043
Processing: OTU.205
Processing: OTU.1180
Processing: OTU.2502
Processing: OTU.8468
Processing: OTU.331
Processing: OTU.409
Processing: OTU.847
Processing: OTU.943
Processing: OTU.158
Processing: OTU.1457
Processing: OTU.66
Processing: OTU.366
Processing: OTU.889
Processing: OTU.2248
Processing: OTU.229
Processing: OTU.142
Processing: OTU.2646
Processing: OTU.3754
Processing: OTU.6739
Processing: OTU.1181
Processing: OTU.6450
Processing: OTU.319
Processing: OTU.115
Processing: OTU.8912
Processing: OTU.9444
Processing: OTU.330
Processing: OTU.5408
Processing: OTU.2658
Processing: OTU.1326
Processing: OTU.1338
Processing: OTU.940
Processing: OTU.367
Processing: OTU.1450
Processing: OTU.3581
Processing: OTU.118
Processing: OTU.10834
Processing: OTU.1846
Processing: OTU.228
Processing: OTU.143
Processing: OTU.2429
Processing: OTU.3287
Processing: OTU.1041
Processing: OTU.3975
Processing: OTU.318
Processing: OTU.221
Processing: OTU.5507
Processing: OTU.9046
Processing: OTU.159
Processing: OTU.4820
Processing: OTU.48
Processing: OTU.4493
Processing: OTU.364
Processing: OTU.7606
Processing: OTU.9479
Processing: OTU.881
Processing: OTU.5129
Processing: OTU.545
Processing: OTU.1849
Processing: OTU.227
Processing: OTU.2401
Processing: OTU.412
Processing: OTU.2539
Processing: OTU.3285
Processing: OTU.1040
Processing: OTU.1189
Processing: OTU.3342
Processing: OTU.2736
Processing: OTU.5294
Processing: OTU.1078
Processing: OTU.2366
Processing: OTU.845
Processing: OTU.49
Processing: OTU.365
Processing: OTU.5993
Processing: OTU.880
Processing: OTU.4912
Processing: OTU.2241
Processing: OTU.2506
Processing: OTU.141
Processing: OTU.7136
Processing: OTU.3282
Processing: OTU.3922
Processing: OTU.1046
Processing: OTU.2473
Processing: OTU.1676
Processing: OTU.486
Processing: OTU.3608
Processing: OTU.120
Processing: OTU.442
Processing: OTU.405
Processing: OTU.5602
Processing: OTU.11850
Processing: OTU.362
Processing: OTU.1594
Processing: OTU.887
Processing: OTU.1447
Processing: OTU.1445
Processing: OTU.2243
Processing: OTU.557
Processing: OTU.3117
Processing: OTU.1198
Processing: OTU.7443
Processing: OTU.1045
Processing: OTU.5755
Processing: OTU.1677
Processing: OTU.6238
Processing: OTU.5290
Processing: OTU.1076
Processing: OTU.6407
Processing: OTU.404
Processing: OTU.6551
Processing: OTU.47
Processing: OTU.121
Processing: OTU.363
Processing: OTU.1292
Processing: OTU.1296
Processing: OTU.885
Processing: OTU.4164
Processing: OTU.3620
Processing: OTU.1703
Processing: OTU.286
Processing: OTU.13145
Processing: OTU.3288
Processing: OTU.8277
Processing: OTU.460
Processing: OTU.3349
Processing: OTU.3603
Processing: OTU.8484
Processing: OTU.1075
Processing: OTU.402
Processing: OTU.2091
Processing: OTU.6558
Processing: OTU.44
Processing: OTU.1340
Processing: OTU.360
Processing: OTU.1297
Processing: OTU.884
Processing: OTU.1440
Processing: OTU.3312
Processing: OTU.3628
Processing: OTU.2084
Processing: OTU.551
Processing: OTU.122
Processing: OTU.966
Processing: OTU.13147
Processing: OTU.300
Processing: OTU.5356
Processing: OTU.1954
Processing: OTU.1678
Processing: OTU.899
Processing: OTU.289
Processing: OTU.4552
Processing: OTU.401
Processing: OTU.842
Processing: OTU.45
Processing: OTU.361
Processing: OTU.5285
Processing: OTU.7298
Processing: OTU.1284
Processing: OTU.6090
Processing: OTU.4957
Processing: OTU.967
Processing: OTU.1192
Processing: OTU.13547
Processing: OTU.527
Processing: OTU.301
Processing: OTU.1920
Processing: OTU.1679
Processing: OTU.891
Processing: OTU.220
Processing: OTU.123
Processing: OTU.1073
Processing: OTU.400
Processing: OTU.7146
Processing: OTU.42
Processing: OTU.368
Processing: OTU.6951
Processing: OTU.799
Processing: OTU.1281
Processing: OTU.436
Processing: OTU.1706
Processing: OTU.960
Processing: OTU.687
Processing: OTU.6884
Processing: OTU.302
Processing: OTU.12589
Processing: OTU.2345
Processing: OTU.895
Processing: OTU.10520
Processing: OTU.1070
Processing: OTU.7840
Processing: OTU.974
Processing: OTU.840
Processing: OTU.125
Processing: OTU.485
Processing: OTU.40
Processing: OTU.369
Processing: OTU.4070
Processing: OTU.798
Processing: OTU.1283
Processing: OTU.434
Processing: OTU.2500
Processing: OTU.961
Processing: OTU.839
Processing: OTU.6885
Processing: OTU.303
Processing: OTU.10985
Processing: OTU.428
Processing: OTU.896
Processing: OTU.7728
Processing: OTU.116
Processing: OTU.447
Processing: OTU.971
Processing: OTU.12406
Processing: OTU.41
Processing: OTU.7848
Processing: OTU.1601
Processing: OTU.942
Processing: OTU.2696
Processing: OTU.435
Processing: OTU.2408
Processing: OTU.1851
Processing: OTU.962
Processing: OTU.3908
Processing: OTU.59
Processing: OTU.1282
Processing: OTU.304
Processing: OTU.6983
Processing: OTU.139
Processing: OTU.897
Processing: OTU.7729
Processing: OTU.110
Processing: OTU.970
Processing: OTU.4481
Processing: OTU.10671
Processing: OTU.1600
Processing: OTU.7738
Processing: OTU.2695
Processing: OTU.432
Processing: OTU.1856
Processing: OTU.1061
Processing: OTU.963
Processing: OTU.3909
Processing: OTU.58
Processing: OTU.305
Processing: OTU.4377
Processing: OTU.138
Processing: OTU.788
Processing: OTU.1519
Processing: OTU.949
Processing: OTU.8395
Processing: OTU.2468
Processing: OTU.6194
Processing: OTU.3919
Processing: OTU.1882
Processing: OTU.1603
Processing: OTU.2559
Processing: OTU.794
Processing: OTU.433
Processing: OTU.5510
Processing: OTU.6187
Processing: OTU.528
Processing: OTU.55
Processing: OTU.306
Processing: OTU.4374
Processing: OTU.421
Processing: OTU.789
Processing: OTU.7282
Processing: OTU.1518
Processing: OTU.11632
Processing: OTU.1060
Processing: OTU.113
Processing: OTU.12027
Processing: OTU.3911
Processing: OTU.4947
Processing: OTU.3333
Processing: OTU.1602
Processing: OTU.3745
Processing: OTU.3742
Processing: OTU.793
Processing: OTU.431
Processing: OTU.2233
Processing: OTU.288
Processing: OTU.680
Processing: OTU.54
Processing: OTU.1058
Processing: OTU.420
Processing: OTU.307
Processing: OTU.786
Processing: OTU.1515
Processing: OTU.11633
Processing: OTU.112
Processing: OTU.12021
Processing: OTU.5770
Processing: OTU.11108
Processing: OTU.4543
Processing: OTU.1604
Processing: OTU.3291
Processing: OTU.2692
Processing: OTU.386
Processing: OTU.4247
Processing: OTU.438
Processing: OTU.2232
Processing: OTU.560
Processing: OTU.683
Processing: OTU.57
Processing: OTU.226
Processing: OTU.3120
Processing: OTU.308
Processing: OTU.787
Processing: OTU.1514
Processing: OTU.10237
Processing: OTU.5366
Processing: OTU.6192
Processing: OTU.5772
Processing: OTU.1429
Processing: OTU.1606
Processing: OTU.203
Processing: OTU.791
Processing: OTU.439
Processing: OTU.2235
Processing: OTU.6715
Processing: OTU.561
Processing: OTU.5765
Processing: OTU.56
Processing: OTU.12196
Processing: OTU.422
Processing: OTU.309
Processing: OTU.780
Processing: OTU.11920
Processing: OTU.2530
Processing: OTU.9264
Processing: OTU.5367
Processing: OTU.572
Processing: OTU.13222
Processing: OTU.5372
Processing: OTU.1428
Processing: OTU.1609
Processing: OTU.7456
Processing: OTU.4598
Processing: OTU.313
Processing: OTU.2727
Processing: OTU.5763
Processing: OTU.562
Processing: OTU.51
Processing: OTU.425
Processing: OTU.3355
Processing: OTU.781
Processing: OTU.11922
Processing: OTU.1065
Processing: OTU.7928
Processing: OTU.4659
Processing: OTU.571
Processing: OTU.7852
Processing: OTU.10056
Processing: OTU.5378
Processing: OTU.2531
Processing: OTU.2757
Processing: OTU.8712
Processing: OTU.7457
Processing: OTU.2725
Processing: OTU.312
Processing: OTU.832
Processing: OTU.563
Processing: OTU.50
Processing: OTU.53
Processing: OTU.424
Processing: OTU.1666
Processing: OTU.3517
Processing: OTU.8003
Processing: OTU.2428
Processing: OTU.5160
Processing: OTU.7851
Processing: OTU.3648
Processing: OTU.79
Processing: OTU.78
Processing: OTU.1064
Processing: OTU.1067
Processing: OTU.3370
Processing: OTU.311
Processing: OTU.333
Processing: OTU.238
Processing: OTU.3631
Processing: OTU.1172
Processing: OTU.833
Processing: OTU.564
Processing: OTU.52
Processing: OTU.9612
Processing: OTU.135
Processing: OTU.1724
Processing: OTU.1665
Processing: OTU.5215
Processing: OTU.1371
Processing: OTU.7858
Processing: OTU.6036
Processing: OTU.5172
Processing: OTU.225
Processing: OTU.2421
Processing: OTU.310
Processing: OTU.862
Processing: OTU.332
Processing: OTU.3680
Processing: OTU.1066
Processing: OTU.12128
Processing: OTU.830
Processing: OTU.565
Processing: OTU.426
Processing: OTU.202
Processing: OTU.3356
Processing: OTU.239
Processing: OTU.448
Processing: OTU.13682
Processing: OTU.9286
Processing: OTU.4922
Processing: OTU.11804
Processing: OTU.10036
Processing: OTU.2351
Processing: OTU.11695
Processing: OTU.2092
Processing: OTU.5955
Processing: OTU.6089
Processing: OTU.2666
Processing: OTU.2009
Processing: OTU.2845
Processing: OTU.5614
Processing: OTU.950
Processing: OTU.2529
Processing: OTU.204
Processing: OTU.1661
Processing: OTU.5015
Processing: OTU.5165
Processing: OTU.4395
Processing: OTU.924
Processing: OTU.5178
Processing: OTU.542
Processing: OTU.11786
Processing: OTU.6416
Processing: OTU.765
Processing: OTU.9923
Processing: OTU.855
Processing: OTU.906
Processing: OTU.7868
Processing: OTU.8021
Processing: OTU.2842
Processing: OTU.223
Processing: OTU.953
Processing: OTU.207
Processing: OTU.3358
Processing: OTU.1738
Processing: OTU.12171
Processing: OTU.60
Processing: OTU.133
Processing: OTU.429
Processing: OTU.697
Processing: OTU.1461
Processing: OTU.109
Processing: OTU.1492
Processing: OTU.4025
Processing: OTU.1712
Processing: OTU.2663
Processing: OTU.6419
Processing: OTU.5714
Processing: OTU.1358
Processing: OTU.2843
Processing: OTU.954
Processing: OTU.206
Processing: OTU.1668
Processing: OTU.61
Processing: OTU.2267
Processing: OTU.1725
Processing: OTU.592
Processing: OTU.550
Processing: OTU.132
Processing: OTU.694
Processing: OTU.1462
Processing: OTU.676
Processing: OTU.6601
Processing: OTU.1493
Processing: OTU.8038
Processing: OTU.4923
Processing: OTU.2394
Processing: OTU.11122
Processing: OTU.8280
Processing: OTU.2708
Processing: OTU.128
Processing: OTU.540
Processing: OTU.62
Processing: OTU.5712
Processing: OTU.2264
Processing: OTU.526
Processing: OTU.423
Processing: OTU.7776
Processing: OTU.695
Processing: OTU.4707
Processing: OTU.443
Processing: OTU.933
Processing: OTU.1469
Processing: OTU.4020
Processing: OTU.1138
Processing: OTU.2395
Processing: OTU.2841
Processing: OTU.2344
Processing: OTU.2258
Processing: OTU.129
Processing: OTU.944
Processing: OTU.237
Processing: OTU.6176
Processing: OTU.63
Processing: OTU.2265
Processing: OTU.130
Processing: OTU.616
Processing: OTU.13207
Processing: OTU.928
Processing: OTU.670
Processing: OTU.684
Processing: OTU.7719
Processing: OTU.7774
Processing: OTU.6607
Processing: OTU.7155
Processing: OTU.3730
Processing: OTU.7934
Processing: OTU.114
Processing: OTU.1468
Processing: OTU.2705
Processing: OTU.2402
Processing: OTU.230
Processing: OTU.3419
Processing: OTU.4410
Processing: OTU.3126
Processing: OTU.9017
Processing: OTU.544
Processing: OTU.671
Processing: OTU.867
Processing: OTU.931
Processing: OTU.1385
Processing: OTU.9355
Processing: OTU.13201
Processing: OTU.1133
Processing: OTU.11361
Processing: OTU.77
Processing: OTU.2256
Processing: OTU.7773
Processing: OTU.4940
Processing: OTU.263
Processing: OTU.482
Processing: OTU.7088
Processing: OTU.117
Processing: OTU.1263
Processing: OTU.9818
Processing: OTU.698
Processing: OTU.6763
Processing: OTU.6770
Processing: OTU.2605
Processing: OTU.866
Processing: OTU.930
Processing: OTU.1386
Processing: OTU.11368
Processing: OTU.8254
Processing: OTU.1131
Processing: OTU.4450
Processing: OTU.5575
Processing: OTU.4943
Processing: OTU.262
Processing: OTU.2488
Processing: OTU.858
Processing: OTU.136
Processing: OTU.157
Processing: OTU.481
Processing: OTU.440
Processing: OTU.2600
Processing: OTU.914
Processing: OTU.11913
Processing: OTU.936
Processing: OTU.12077
Processing: OTU.699
Processing: OTU.1387
Processing: OTU.9449
Processing: OTU.1007
Processing: OTU.328
Processing: OTU.75
Processing: OTU.74
Processing: OTU.4457
Processing: OTU.2254
Processing: OTU.1714
Processing: OTU.261
Processing: OTU.1126
Processing: OTU.1472
Processing: OTU.427
Processing: OTU.1474
Processing: OTU.2489
Processing: OTU.6567
Processing: OTU.905
Processing: OTU.469
Processing: OTU.935
Processing: OTU.11860
Processing: OTU.1380
Processing: OTU.489
Processing: OTU.5264
Processing: OTU.111
Processing: OTU.73
Processing: OTU.1136
Processing: OTU.1260

Making an incorp config file


In [ ]:
!cd $workDir; \
    SIPSim incorpConfigExample \
    --percTaxa 0 \
    --percIncorpUnif 100 \
    > PT0_PI100.config

Adding isotope incorporation to BD distribution


In [ ]:
!cd $workDir; \
    SIPSim isotope_incorp \
    ampFrags_wRand_kde_dif.pkl \
    PT0_PI100.config \
    --comm comm.txt \
    --np $nprocs \
    > ampFrags_wRand_kde_dif_incorp.pkl

Calculating BD shift from isotope incorporation


In [ ]:
!cd $workDir; \
    SIPSim BD_shift \
    ampFrags_wRand_kde_dif.pkl \
    ampFrags_wRand_kde_dif_incorp.pkl \
    --np $nprocs \
    > ampFrags_wRand_kde_dif_incorp_BD-shift.txt

Simulating gradient fractions


In [ ]:
!cd $workDir; \
    SIPSim gradient_fractions \
    comm.txt \
    > fracs.txt

Simulating an OTU table


In [ ]:
!cd $workDir; \
    SIPSim OTU_table \
    ampFrags_wRand_kde_dif_incorp.pkl \
    comm.txt \
    fracs.txt \
    --abs 1e9 \
    --np $nprocs \
    > OTU_abs1e9.txt

Plotting taxon abundances


In [311]:
%%R -i workDir
setwd(workDir)

# loading file
tbl = read.delim('OTU_abs1e9.txt', sep='\t')

In [45]:
%%R
## BD for G+C of 0 or 100
BD.GCp0 = 0 * 0.098 + 1.66
BD.GCp100 = 1 * 0.098 + 1.66

In [313]:
%%R -w 800 -h 300
# plotting absolute abundances

tbl.s = tbl %>%
    group_by(library, BD_mid) %>%
    summarize(total_count = sum(count))

## plot
p = ggplot(tbl.s, aes(BD_mid, total_count)) +
    geom_area(stat='identity', alpha=0.3, position='dodge') +
    geom_histogram(stat='identity') +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16) 
    )
p



In [314]:
%%R -w 800 -h 300
# plotting number of taxa at each BD

tbl.nt = tbl %>%
    filter(count > 0) %>%
    group_by(library, BD_mid) %>%
    summarize(n_taxa = n())

## plot
p = ggplot(tbl.nt, aes(BD_mid, n_taxa)) +
    geom_area(stat='identity', alpha=0.3, position='dodge') +
    geom_histogram(stat='identity') +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density', y='Number of taxa') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p



In [315]:
%%R -w 800 -h 250
# plotting relative abundances

## plot
p = ggplot(tbl, aes(BD_mid, count, fill=taxon)) +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p + geom_area(stat='identity', position='dodge', alpha=0.5)



In [316]:
%%R -w 800 -h 250

p + geom_area(stat='identity', position='fill')


Subsampling from the OTU table


In [37]:
dist,loc,scale = seq_per_fraction

!cd $workDir; \
    SIPSim OTU_subsample \
    --dist $dist \
    --dist_params mean:$loc,sigma:$scale \
    --walk 1 \
    OTU_abs1e9.txt \
    > OTU_abs1e9_sub.txt

Testing/Plotting seq count distribution of subsampled fraction samples


In [38]:
%%R -h 300
setwd(workDir)

tbl = read.csv('OTU_abs1e9_sub.txt', sep='\t') 

tbl.s = tbl %>% 
    group_by(library, fraction) %>%
    summarize(total_count = sum(count)) %>%
    ungroup() %>%
    mutate(library = as.character(library))

ggplot(tbl.s, aes(total_count)) +
    geom_density(fill='blue')



In [39]:
%%R -h 300 -w 600
setwd(workDir)

tbl.s = tbl %>%
    group_by(fraction, BD_min, BD_mid, BD_max) %>%
    summarize(total_count = sum(count)) 

ggplot(tbl.s, aes(BD_mid, total_count)) +
    geom_point() +
    geom_line() +
    labs(x='Buoyant density', y='Total sequences') +
    theme_bw() +
    theme(
        text = element_text(size=16)
        )


Getting list of target taxa


In [41]:
%%R -i workDir

inFile = paste(c(workDir, 'target_genome_index.txt'), collapse='/')

tbl.target = read.delim(inFile, sep='\t', header=F)
colnames(tbl.target) = c('OTUId', 'genome_file', 'genome_ID', 'X', 'Y', 'Z')
tbl.target = tbl.target %>% distinct(OTUId)


cat('Number of target OTUs: ', tbl.target$OTUId %>% unique %>% length, '\n')
cat('----------\n')
tbl.target %>% head(n=3)


Number of target OTUs:  108 
----------
     OTUId
1 OTU.1457
2   OTU.77
3 OTU.2347
                                                                                       genome_file
1    /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Thermomonospora_curvata_DSM_43183.fasta
2      /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Pseudoxanthomonas_spadix_BD-a59.fasta
3 /home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/Micromonospora_aurantiaca_ATCC_27029.fasta
                                      genome_ID X           Y    Z
1    CP001738_Thermomonospora_curvata_DSM_43183 1 0.009533343 1038
2      CP003093_Pseudoxanthomonas_spadix_BD_a59 1 0.019066686  700
3 CP002162_Micromonospora_aurantiaca_ATCC_27029 1 0.004766671 1380

Plotting abundance distributions


In [46]:
%%R -w 800 -h 250
# plotting relative abundances

tbl = tbl %>% 
    group_by(fraction) %>%
    mutate(rel_abund = count / sum(count))


## plot
p = ggplot(tbl, aes(BD_mid, count, fill=taxon)) +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p + geom_area(stat='identity', position='dodge', alpha=0.5)



In [47]:
%%R -w 800 -h 250

p = ggplot(tbl, aes(BD_mid, rel_abund, fill=taxon)) +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p + geom_area(stat='identity')


Abundance distribution of just target taxa


In [48]:
%%R

targets = tbl.target$OTUId %>% as.vector %>% unique 

tbl.f = tbl %>%
    filter(taxon %in% targets)

tbl.f %>% head


Source: local data frame [6 x 8]
Groups: fraction

  library    fraction taxon BD_min BD_mid BD_max count rel_abund
1       1  -inf-1.660 OTU.1   -Inf  1.660  1.660     0         0
2       1 1.660-1.664 OTU.1  1.660  1.662  1.664     0         0
3       1 1.664-1.666 OTU.1  1.664  1.665  1.666     0         0
4       1 1.666-1.670 OTU.1  1.666  1.668  1.670     0         0
5       1 1.670-1.676 OTU.1  1.670  1.673  1.676     0         0
6       1 1.676-1.679 OTU.1  1.676  1.677  1.679     0         0

In [49]:
%%R -w 800 -h 250
# plotting absolute abundances

## plot
p = ggplot(tbl.f, aes(BD_mid, count, fill=taxon)) +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p + geom_area(stat='identity', position='dodge', alpha=0.5)



In [50]:
%%R -w 800 -h 250
# plotting relative abundances

p = ggplot(tbl.f, aes(BD_mid, rel_abund, fill=taxon)) +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p + geom_area(stat='identity')


Plotting 'true' taxon abundance distribution (from priming exp dataset)


In [51]:
%%R -i metaDataFile
# loading priming_exp metadata file

meta = read.delim(metaDataFile, sep='\t')
meta %>% head(n=4)


            Sample FractionNum Bulk Control CC X100 X700 H2O Day Density rep
1 12C.000.28.03.07           7    0       1  1    0    0   0  28  1.7646    
2 12C.000.28.03.08           8    0       1  1    0    0   0  28  1.7614    
3 12C.000.28.03.09           9    0       1  1    0    0   0  28  1.7537    
4 12C.000.28.03.10          10    0       1  1    0    0   0  28  1.7483    
  contolVlabel Treatment
1      control    12C000
2      control    12C000
3      control    12C000
4      control    12C000

In [52]:
%%R -i otuTableFile
# loading priming_exp OTU table 

tbl.otu.true = read.delim(otuTableFile, sep='\t') %>%
    select(OTUId, starts_with('X12C.700.14')) 
tbl.otu.true %>% head(n=3)


     OTUId X12C.700.14.06.14 X12C.700.14.06.05 X12C.700.14.06.16
1 OTU.4776                34               142                11
2 OTU.2864                 1                16                 1
3 OTU.8170                 0                 0                 0
  X12C.700.14.06.11 X12C.700.14.06.13 X12C.700.14.06.09 X12C.700.14.06.04
1                18                30                27                97
2                 2                 0                 1                 3
3                 0                 0                 0                 0
  X12C.700.14.06.20 X12C.700.14.06.10 X12C.700.14.06.15 X12C.700.14.06.08
1               226                23                19                35
2                 1                 0                 1                 0
3                 0                 1                 0                 0
  X12C.700.14.06.12 X12C.700.14.06.18 X12C.700.14.06.17 X12C.700.14.06.19
1                25                61                27               136
2                 0                 2                 2                 0
3                 0                 0                 0                 0
  X12C.700.14.06.23 X12C.700.14.06.21 X12C.700.14.06.06 X12C.700.14.06.07
1                85               283                50                 0
2                 0                 1                 0                 0
3                 0                 0                 0                 0
  X12C.700.14.06.22 X12C.700.14.05.NA X12C.700.14.06.NA
1               208                51                26
2                 2                 0                 0
3                 0                 0                 0

In [53]:
%%R
# editing table
tbl.otu.true.w = tbl.otu.true %>%
    gather('sample', 'count', 2:ncol(tbl.otu.true)) %>%
    mutate(sample = gsub('^X', '', sample)) %>%
    group_by(sample) %>%
    mutate(rel_abund = count / sum(count)) %>%
    ungroup() %>%
    filter(count > 0)
tbl.otu.true.w %>% head(n=5)


Source: local data frame [5 x 4]

     OTUId           sample count    rel_abund
1 OTU.4776 12C.700.14.06.14    34 7.263716e-04
2 OTU.2864 12C.700.14.06.14     1 2.136387e-05
3 OTU.5406 12C.700.14.06.14     2 4.272774e-05
4 OTU.6001 12C.700.14.06.14     1 2.136387e-05
5  OTU.611 12C.700.14.06.14    23 4.913690e-04

In [54]:
%%R
tbl.true.j = inner_join(tbl.otu.true.w, meta, c('sample' = 'Sample'))
tbl.true.j %>% as.data.frame %>% head(n=3)


     OTUId           sample count    rel_abund FractionNum Bulk Control CC X100
1 OTU.4776 12C.700.14.06.14    34 7.263716e-04          14    0       1  0    0
2 OTU.2864 12C.700.14.06.14     1 2.136387e-05          14    0       1  0    0
3 OTU.5406 12C.700.14.06.14     2 4.272774e-05          14    0       1  0    0
  X700 H2O Day Density rep contolVlabel Treatment
1    1   0  14  1.7122          control    12C700
2    1   0  14  1.7122          control    12C700
3    1   0  14  1.7122          control    12C700

In [55]:
%%R -w 800 -h 300
# plotting number of taxa at each BD

tbl.nt = tbl %>%
    filter(count > 0) %>%
    group_by(library, BD_mid) %>%
    summarize(n_taxa = n())

## plot
p = ggplot(tbl.nt, aes(BD_mid, n_taxa)) +
    geom_area(stat='identity', alpha=0.5) +
    geom_point() +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density', y='Number of taxa') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )
p



In [56]:
%%R -w 700 -h 350

tbl.true.j.s = tbl.true.j %>%
    filter(count > 0) %>%
    group_by(sample, Density) %>%
    summarize(n_taxa = sum(count > 0))

ggplot(tbl.true.j.s, aes(Density, n_taxa)) +
    geom_area(stat='identity', alpha=0.5) +
    geom_point() +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density', y='Number of taxa') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )


Plotting total counts for each sample


In [57]:
%%R -h 300 -w 600
tbl.true.j.s = tbl.true.j %>%
    group_by(sample, Density) %>%
    summarize(total_count = sum(count)) 

ggplot(tbl.true.j.s, aes(Density, total_count)) +
    geom_point() +
    geom_line() +
    labs(x='Buoyant density', y='Total sequences') +
    theme_bw() +
    theme(
        text = element_text(size=16)
        )


Plotting abundance distribution of target OTUs


In [58]:
%%R
tbl.true.j.f = tbl.true.j %>%
    filter(OTUId %in% targets) %>%
    arrange(OTUId, Density) %>%
    group_by(sample)
tbl.true.j.f %>% head(n=3) %>% as.data.frame


  OTUId           sample count  rel_abund FractionNum Bulk Control CC X100 X700
1 OTU.1 12C.700.14.06.23   736 0.02455871          23    0       1  0    0    1
2 OTU.1 12C.700.14.06.22  1229 0.02870556          22    0       1  0    0    1
3 OTU.1 12C.700.14.06.21  1560 0.03320067          21    0       1  0    0    1
  H2O Day Density rep contolVlabel Treatment
1   0  14  1.6772          control    12C700
2   0  14  1.6805          control    12C700
3   0  14  1.6849          control    12C700

In [59]:
%%R -w 800 -h 250
# plotting relative abundances

## plot
ggplot(tbl.true.j.f, aes(Density, rel_abund, fill=OTUId)) +
    geom_area(stat='identity') +
    geom_vline(xintercept=c(BD.GCp0, BD.GCp100), linetype='dashed', alpha=0.5) +
    labs(x='Buoyant density') +
    theme_bw() +
    theme( 
        text = element_text(size=16),
        legend.position = 'none'
    )


Combining true and simulated OTU tables for target taxa


In [60]:
%%R
tbl.f.e = tbl.f %>%
    mutate(library = 'simulation') %>%
    rename('density' = BD_mid) %>%
    select(-BD_min, -BD_max)

tbl.true.e = tbl.true.j.f %>% 
    select('taxon' = OTUId,
           'fraction' = sample,
           'density' = Density,
           count, rel_abund) %>%
    mutate(library = 'true') 
           
           
tbl.sim.true = rbind(tbl.f.e, tbl.true.e)
tbl.f.e = data.frame()
tbl.true.e = data.frame()

tbl.sim.true %>% head(n=3)


Source: local data frame [3 x 6]
Groups: fraction

     library    fraction taxon density count rel_abund
1 simulation  -inf-1.660 OTU.1   1.660     0         0
2 simulation 1.660-1.664 OTU.1   1.662     0         0
3 simulation 1.664-1.666 OTU.1   1.665     0         0

In [62]:
%%R
# check
cat('Number of target taxa: ', tbl.sim.true$taxon %>% unique %>% length, '\n')


Number of target taxa:  108 

Abundance distributions of each target taxon


In [ ]:
%%R -w 900 -h 3000

tbl.sim.true = tbl.sim.true %>%
    ungroup() %>%
    group_by(taxon) %>%
    mutate(mean_rel_abund = mean(rel_abund)) %>%
    ungroup()

tbl.sim.true$taxon = reorder(tbl.sim.true$taxon, -tbl.sim.true$mean_rel_abund)

ggplot(tbl.sim.true, aes(density, rel_abund, color=library)) +
    geom_point() +
    geom_line() +
    theme_bw() +
    facet_wrap(~ taxon, ncol=4, scales='free_y')



In [ ]:


In [ ]:


In [ ]:

sandbox


In [504]:
def partial_shuffle(x, n, frac=False):
    # if frac is True, n assumed to be fraction of len(x)
    if frac is True:
        n = int(round(len(x) * frac,0))
    else:
        n = int(n)
    
    # x = list
    try:
        x = list(x)
    except TypeError:
        msg = 'x must be a list-like object'
        raise TypeError, msg
        
        
    # making list of random indices
    rand = lambda x: random.randrange(0,len(x),1)
    randints = lambda x,n: [rand(x) for i in xrange(n)]
    rand_i = zip(randints(x,n), randints(x,n))
    
    # suffle
    for i,j in rand_i:
        tmp = x[i]
        x[i] = x[j]
        x[j] = tmp
        
    return x

In [506]:
partial_shuffle(range(20), 2.0)


Out[506]:
[0, 8, 2, 3, 4, 5, 18, 7, 1, 9, 10, 11, 12, 13, 14, 15, 16, 17, 6, 19]

In [507]:
import scipy.stats as ss

In [582]:
def shuffle_walk(x, max_walk):
    # x = list
    try:
        x = list(x)
    except TypeError:
        msg = 'x must be a list-like object'
        raise TypeError, msg
        
    x_len = len(x)
    
    # ranks of values
    ## (index, value, value_rank)
    ranks = zip(range(x_len), x, ss.rankdata(x))
 
    # starting range
    start_i = random.randrange(0,x_len,1)
    cur_rank = ranks[start_i]

    # filtering cur_rank from ranks 
    ranks = [x for x in ranks if x[0] != cur_rank[0]]
    
    # moving through ranks
    x_new_order = [cur_rank[1]]
    for i in xrange(x_len-1):       
        # select max walk distance
        max_range = random.randrange(1,max_walk+1)
    
        # filter to just ranks w/in rank distance
        filt_ranks = [x for x in ranks 
                      if abs(x[2] - cur_rank[2]) <= max_range]
    
        # selecting randomly from filtered ranks
        cur_rank = random.sample(filt_ranks, k=1)[0]
        
        # adding value to new list
        x_new_order.append(cur_rank[1])
        
        # filtering cur_rank from ranks 
        ranks = [x for x in ranks if x[0] != cur_rank[0]]

        # re-ranking remaining values
        rank_vals = [x[2] for x in ranks]
        new_ranks = ss.rankdata(rank_vals)
        for i,v in enumerate(new_ranks):
            ranks[i] = (ranks[i][0], ranks[i][1], v)
                    
    return x_new_order

In [637]:
vals = range(1, 100)
sw = shuffle_walk(vals, 3)
random.shuffle(vals)

In [638]:
%%R -w 700 -h 300 -i sw -i vals

tbl = data.frame('shuffle_walk' = sw,
                 'shuffle' = vals)

#ggplot(tbl, aes(shuffle_walk, shuffle)) +
#    geom_point() 

tbl = tbl %>% 
    mutate(order = 1:nrow(tbl)) %>%
    gather('algorithim', 'value', 1:2)

ggplot(tbl, aes(order, value, color=algorithim)) +
    geom_point() +
    geom_line() +
    theme_bw() +
    theme(
        text = element_text(size=16)
        )



In [ ]: