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 [64]:
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


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

In [66]:
%load_ext rpy2.ipython

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

In [67]:

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

Creating a community file from OTU table

In [69]:
%%R -i otuTableFile 

tbl.otu = read.delim(otuTableFile, sep='\t')

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

In [70]:
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) %>%
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 [71]:
%%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 [72]:
%%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() +
        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) %>%
taxonMap %>% nrow %>% print
taxonMap %>% head(n=3) %>% print

breaker = '----------------\n'

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


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
1                  CP003093_Pseudoxanthomonas_spadix_BD_a59
2                  CP000511_Mycobacterium_vanbaalenii_PYR_1
3 CP003344_Desulfitobacterium_dichloroeliminans_LMG_P_21439
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]:
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]:

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)

1 OTU.8540
2 OTU.9267
3 OTU.1457
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
1 CP001738_Thermomonospora_curvata_DSM_43183
2 CP001738_Thermomonospora_curvata_DSM_43183
3 CP001738_Thermomonospora_curvata_DSM_43183

In [22]:
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]:
tbl.j2 = inner_join(tbl.j, comm, c('OTU' = 'taxon_name')) = tbl.j2$OTU %>% unique %>% length
cat('Number of target OTUs: ',, '\n')
cat('--------', '\n')
tbl.j2 %>% head(n=3)

Number of target OTUs:  108 
1 OTU.1457
2   OTU.77
3 OTU.2347
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() +
        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='/') = read.delim(inFile, sep='\t', header=F)
colnames( = 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]:
# just OTUs w/out an associated genome
tbl.j = anti_join(tbl.comm,, c('taxon_name' = 'OTUId'))
n.nontarget.genomes = tbl.j$taxon_name %>% length
message('Number of non-target genomes: ', n.nontarget.genomes)
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.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))

Number of non-target OTUs: 2062
['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:
    ampFrag_rand = []

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

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

In [45]:
## 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() +
        text = element_text(size=16) 

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() +
        text = element_text(size=16),
        legend.position = 'none'

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() +
        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

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)) +

In [39]:
%%R -h 300 -w 600

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() +
        text = element_text(size=16)

Getting list of target taxa

In [41]:
%%R -i workDir

inFile = paste(c(workDir, 'target_genome_index.txt'), collapse='/') = read.delim(inFile, sep='\t', header=F)
colnames( = c('OTUId', 'genome_file', 'genome_ID', 'X', 'Y', 'Z') = %>% distinct(OTUId)

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

Number of target OTUs:  108 
1 OTU.1457
2   OTU.77
3 OTU.2347
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() +
        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() +
        text = element_text(size=16),
        legend.position = 'none'
p + geom_area(stat='identity')

Abundance distribution of just target taxa

In [48]:

targets =$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() +
        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() +
        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.           7    0       1  1    0    0   0  28  1.7646    
2 12C.           8    0       1  1    0    0   0  28  1.7614    
3 12C.           9    0       1  1    0    0   0  28  1.7537    
4 12C.          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]:
# 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]:
tbl.true.j = inner_join(tbl.otu.true.w, meta, c('sample' = 'Sample'))
tbl.true.j %>% %>% 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() +
        text = element_text(size=16),
        legend.position = 'none'

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() +
        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() +
        text = element_text(size=16)

Plotting abundance distribution of target OTUs

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

  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() +
        text = element_text(size=16),
        legend.position = 'none'

Combining true and simulated OTU tables for target taxa

In [60]:
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]:
# 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 [63]:
%%R -w 900 -h 3000

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

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')