<jupyter_start><jupyter_text>Validation against Wawrik et al. 2009 This notebook compares metaSIPSim simulation results against those of Wawrik et al. 2009. Comparisons are made to figure 3B from said paper.For the simulation, I will be generating fragments only. I will estimate the amount of DNA from each genome in each fraction by multiplying fragment size by simulated abundance then summing within each genome.<NAME>., <NAME>., & <NAME>. (2009). Use of Inorganic and Organic Nitrogen by Synechococcus spp. and Diatoms on the West Florida Shelf as Measured Using Stable Isotope Probing. Applied and Environmental Microbiology, 75(21), 6662–6670. DOI:10.1128/AEM.01002-09  Intitiation<jupyter_code>import os
import pandas as pd
import ConfigParser
import rpy2.ipython
%load_ext rpy2.ipython
%%R
library(dplyr)
library(ggplot2)
# What is your working directory?
workdir = '/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/'

# Where are your genomes?
genomedir = '/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/genomes/'


# making directories
## working directory
if not os.path.isdir(workdir):
    os.makedirs(workdir)
%cd $workdir

## genome directory
if not os.path.isdir(genomedir):
    os.makedirs(genomedir)
print(genomedir)<jupyter_output>/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009
/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/genomes/<jupyter_text>Part 1: Get genomes<jupyter_code>%cd $genomedir
# The genomes you want are:
## Synechococcus_sp_WH7803: GCF_000063505.1


!wget -O Synechococcus_sp_WH7803.fna.gz ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/063/505/GCF_000063505.1_ASM6350v1/GCF_000063505.1_ASM6350v1_genomic.fna.gz

!gunzip Synechococcus_sp_WH7803.fna.gz<jupyter_output>/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/genomes
--2019-05-22 13:09:36--  ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/063/505/GCF_000063505.1_ASM6350v1/GCF_000063505.1_ASM6350v1_genomic.fna.gz
           => ‘Synechococcus_sp_WH7803.fna.gz’
Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 172.16.58.3, fdf8:f53e:61e4::18
Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|172.16.58.3|:21... connected.
Logging in as anonymous ... Logged in!
==> SYST ... done.    ==> PWD ... done.
==> TYPE I ... done.  ==> CWD (1) /genomes/all/GCF/000/063/505/GCF_000063505.1_ASM6350v1 ... done.
==> SIZE GCF_000063505.1_ASM6350v1_genomic.fna.gz ... 694327
==> PASV ... done.    ==> RETR GCF_000063505.1_ASM6350v1_genomic.fna.gz ... done.
Length: 694327 (678K) (unauthoritative)

GCF_000063505.1_ASM 100%[===================>] 678.05K  --.-KB/s    in 0.07s   

2019-05-22 13:09:37 (8.92 MB/s) - ‘Synechococcus_sp_WH7803.fna.gz’ saved [694327]<jupyter_text>Part 2: Simulation  A) Get genome indexes<jupyter_code>%cd $workdir
# make index file of genomes 
!ls $genomedir | grep ".fna" > genome_files.txt
!tail genome_files.txt
with open("genome_files.txt", 'r') as genome_files, open ("genome_index.txt", 'w') as genome_index:
    genomes = genome_files.readlines()
    for line in genomes:
        fasta = line
        taxa = fasta.replace('.fna\n', '\t')
        genome_index.write(taxa + fasta)
!cat genome_index.txt<jupyter_output>Synechococcus_sp_WH7803	Synechococcus_sp_WH7803.fna<jupyter_text>B) Set up experimental design  List of incorporatorsIn this case Synechococcus_sp_WH7803 will be 100% incorporating in library 2<jupyter_code>library = [2]

taxon_name = ['Synechococcus_sp_WH7803'] 

percent_incorporation = [100.0]

sd_incorporation = [0]

incorporators = pd.DataFrame({'taxon_name':taxon_name, 'library':library, 
                       'percent_incorporation':percent_incorporation, 'sd_incorporation':sd_incorporation})
incorporators = incorporators.reindex(columns = ['taxon_name', 'library', 'percent_incorporation', 'sd_incorporation'])
incorporators.to_csv("incorporators.txt", sep='\t', index=False)

!cat incorporators.txt<jupyter_output>taxon_name	library	percent_incorporation	sd_incorporation
Synechococcus_sp_WH7803	2	100.0	0<jupyter_text>Community composition fileCultures were isolates so they have a relative abundance of 100%<jupyter_code>library = [1, 2]

taxon_name = ['Synechococcus_sp_WH7803', 'Synechococcus_sp_WH7803'] 

rel_abund_perc = [100.0, 100.0]

rank = [1,2]

comm = pd.DataFrame({'library':library, 'taxon_name':taxon_name, 
                       'rel_abund_perc':rel_abund_perc, 'rank':rank})
comm = comm.reindex(columns = ['library', 'taxon_name', 'rel_abund_perc', 'rank'])
comm.to_csv("comm.txt", sep='\t', index=False)

!cat comm.txt<jupyter_output>library	taxon_name	rel_abund_perc	rank
1	Synechococcus_sp_WH7803	100.0	1
2	Synechococcus_sp_WH7803	100.0	2<jupyter_text>Generate fraction fileThese fractions are estimated from the paper figure<jupyter_code>%%R
fractions = read.table('/home/sam/data/SIPSim2_data/dev/validation/figure_estimates.txt', 
                          header=TRUE, sep="\t") %>%
    filter(Study == "Wawrik2009") %>%
    mutate(BD = round(BD, digits=3)) %>%
    select(library, BD, fraction) %>%
    group_by(library) %>%
    mutate(last_frac = max(fraction)) %>%
    ungroup()
for (i in 1:nrow(fractions)){
    fractions$BD_prev[i] = fractions$BD[i-1]
    fractions$BD_next[i] = fractions$BD[i+1]
    if (fractions$fraction[i] == 1){
        fractions$BD_prev[i] = (2*fractions$BD[i])-fractions$BD_next[i]
    }
    if (fractions$fraction[i] == fractions$last_frac[i]){
        fractions$BD_next[i] = (2*fractions$BD[i])-fractions$BD_prev[i]
    }
}
fractions = fractions %>%
    mutate(BD_min = (BD+BD_prev)/2, BD_max = (BD+BD_next)/2) %>%
    mutate(fraction_size = round(BD_max-BD_min, 4)) %>%
    select(library, fraction, BD_min, BD_max, fraction_size)
write.table(fractions, "fracs.txt", quote=F, sep="\t", row.names=F)
!cat fracs.txt<jupyter_output>library	fraction	BD_min	BD_max	fraction_size
1	1	1.687	1.693	0.006
1	2	1.693	1.6975	0.0045
1	3	1.6975	1.7005	0.003
1	4	1.7005	1.704	0.0035
1	5	1.704	1.708	0.004
1	6	1.708	1.712	0.004
1	7	1.712	1.7165	0.0045
1	8	1.7165	1.721	0.0045
1	9	1.721	1.725	0.004
1	10	1.725	1.7295	0.0045
1	11	1.7295	1.7345	0.005
1	12	1.7345	1.739	0.0045
1	13	1.739	1.7425	0.0035
1	14	1.7425	1.7455	0.003
2	1	1.687	1.693	0.006
2	2	1.693	1.6965	0.0035
2	3	1.6965	1.698	0.0015
2	4	1.698	1.7005	0.0025
2	5	1.7005	1.704	0.0035
2	6	1.704	1.7075	0.0035
2	7	1.7075	1.7105	0.003
2	8	1.7105	1.713	0.0025
2	9	1.713	1.7155	0.0025
2	10	1.7155	1.718	0.0025
2	11	1.718	1.7235	0.0055
2	12	1.7235	1.73	0.0065
2	13	1.73	1.7335	0.0035
2	14	1.7335	1.7365	0.003
2	15	1.7365	1.74	0.0035
2	16	1.74	1.7445	0.0045
2	17	1.7445	1.7495	0.005
2	18	1.7495	1.755	0.0055
2	19	1.755	1.761	0.006<jupyter_text>C) Make configuration file for simulationThis configuration file contains all the information needed for simulation. Gradient configuration is based on the methods section.<jupyter_code>config = ConfigParser.SafeConfigParser()

## Other parameters
config.add_section('Other')
config.set('Other', 'temp_directory', './tmp')
config.set('Other', 'threads', '1')
config.set('Other', 'logfile', 'SIPSim2_Wawrik.log')
config.set('Other', 'endpoint', 'fragment_list')

## Library parameters
config.add_section('Library')
config.set('Library', 'library_list', '1, 2')
config.set('Library', 'window_or_fraction', 'fraction')

## Fragment parameters
config.add_section('Fragment')
config.set('Fragment', 'path_to_SIPSim', '/home/sam/notebooks/SIPSim')
config.set('Fragment', 'genome_index_file', 'genome_index.txt')
config.set('Fragment', 'genomeDir', genomedir)
config.set('Fragment', 'frag_length_distribution', 'skewed-normal,9000,2500,-5')
config.set('Fragment', 'coverage_of_fragments', '100')
config.set('Fragment', 'fragment_file', 'frags.txt')
config.set('Fragment', 'temp_fragment_file', 'tmp.frags')
config.set('Fragment', 'number_of_iterations', '1')


## Gradient parameters
config.add_section('Gradient')
config.set('Gradient', 'temperature', '293.15')
config.set('Gradient', 'avg_density', '1.701')
config.set('Gradient', 'angular_velocity', '16860743')
config.set('Gradient', 'min_rotation_radius', '7.47')
config.set('Gradient', 'max_rotation_radius', '8.79')
config.set('Gradient', 'tube_angle', '0')
config.set('Gradient', 'tube_radius', '0.66')
config.set('Gradient', 'tube_height', '4.7')
config.set('Gradient', 'fraction_frag_in_DBL', '0.001')
config.set('Gradient', 'isotope', 'N')

## Model parameters
config.add_section('Model')
config.set('Model', 'min_bouyant_density', '1.67')
config.set('Model', 'max_bouyant_density', '1.775')
config.set('Model', 'bouyant_density_step', '0.0001')
config.set('Model', 'fraction_table_file', 'fracs.txt')

## Community parameters
config.add_section('Community')
config.set('Community', 'community_file', 'comm.txt')
config.set('Community', 'incorporator_file', 'incorporators.txt')

## Sequencing parameters
config.add_section('Sequencing')
config.set('Sequencing', 'number_of_genomes_per_iteration', '1')


# Writing our configuration file to 'example.cfg'
with open('simulation_parameters.cfg', 'wb') as configfile:
    config.write(configfile)<jupyter_output><empty_output><jupyter_text>D) Simulate fragments for each fractionUse the function SIPSim2_get_fragment_list_fractions.py to get the fragments for each genome in each fraction of every library<jupyter_code>!python /home/sam/notebooks/SIPSim_metagenome/bin/SIPSim_metagenome.py simulation_parameters.cfg<jupyter_output>Running SIPSim_metagenome
This program was writen by <NAME> (<EMAIL>)

This run was started on 22/05/19 at 13:10:17


You have chosen to get an list of simulated SIP metagenome fragments with adjusted buoyant densities.
You have selected to simulate metagenome sequences for each fraction.

Your community abundance file is: comm.txt

Your incorporator assignment file is: incorporators.txt

You are simulating with the isotope of N

It took 0.734 seconds to get these models.
It took 0.792 seconds to get these models.

Building fragments

It took 12.16 seconds to build the fragments

----------

Starting library 1

Starting library 1 fraction 1 BD:1.687-1.693
Writing fragments to file
It took 13.135 seconds to write fragment file.
It took 13.136 seconds to run library 1 fraction BD:1.687-1.693 iteration 1

It took 13.136 seconds to run the library 1 fraction 1 BD:1.687-1.693

Starting library 1 fraction 2 BD:1.693-1.6975
Writing fragments to file
It took 12.609 seconds t[...]<jupyter_text>Part 3: Original SIPSim simulation Compare SIPSim_metagenome to SIPSim<jupyter_code>%cd $workdir<jupyter_output>/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009<jupyter_text>Modify the comm.txt file so that you include both genomes<jupyter_code>%%bash
echo -e "library\ttaxon_name\trel_abund_perc\trank" > comm2.txt
echo -e "1\tSynechococcus_sp_WH7803\t100.0\t1" >> comm2.txt
echo -e "2\tSynechococcus_sp_WH7803\t100.0\t1" >> comm2.txt


cat comm2.txt<jupyter_output>library	taxon_name	rel_abund_perc	rank
1	Synechococcus_sp_WH7803	100.0	1
2	Synechococcus_sp_WH7803	100.0	1<jupyter_text>1) Get fragments<jupyter_code>!SIPSim fragments \
    genome_index.txt \
    --fp $genomedir/ \
    --fld skewed-normal,9000,2500,-5 \
    --flr None,None \
    --nf 100X \
    --tbl \
    --np 20 \
    2> Wawrik_SIPsim.log \
    > Wawrik_SIPsim.txt
!cat Wawrik_SIPsim.log<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
Processing: "Synechococcus_sp_WH7803"
  Genome name: Synechococcus_sp_WH7803
  Genome length (bp): 2366980
  Number of amplicons: None
  Number of fragments simulated: 33621<jupyter_text>2) Convert to 2d-KDE format<jupyter_code>!SIPSim fragment_KDE \
    Wawrik_SIPsim.txt \
    > Wawrik_SIPsim_kde.pkl
    
!ls -thlc Wawrik_SIPsim_kde.pkl<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
-rw-rw-r-- 1 sam sam 527K May 22 13:22 Wawrik_SIPsim_kde.pkl<jupyter_text>3) Add diffusion<jupyter_code>!SIPSim diffusion \
    Wawrik_SIPsim_kde.pkl \
    --np 20 \
    > Wawrik_SIPsim_kde_dif.pkl

!ls -thlc Wawrik_SIPsim_kde_dif.pkl<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
Index size: 90508
Processing: Synechococcus_sp_WH7803
-rw-rw-r-- 1 sam sam 3.9M May 22 13:22 Wawrik_SIPsim_kde_dif.pkl<jupyter_text>4) Add diffusive boundery layer<jupyter_code>!SIPSim DBL \
    Wawrik_SIPsim_kde_dif.pkl \
    --np 20 \
    -D 1.701 \
    -w 16860743 \
    --tube_height 4.7 \
    --r_min 7.47 \
    --r_max 8.79 \
    --vertical \
    > Wawrik_SIPsim_kde_dif_DBL.pkl


!ls -thlc Wawrik_SIPsim_kde_dif_DBL.pkl<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
/home/sam/anaconda3/envs/SIPSim2_py2/lib/python2.7/site-packages/numpy/lib/polynomial.py:680: RuntimeWarning: invalid value encountered in multiply
  y = y * x + p[i]
DBL_index file written: "DBL_index.txt"
Processing: Synechococcus_sp_WH7803
-rw-rw-r-- 1 sam sam 3.9M May 22 13:22 Wawrik_SIPsim_kde_dif_DBL.pkl<jupyter_text>5) Add isotope incorporation<jupyter_code>!cat incorporators.txt
%%bash
SIPSim incorp_config_example \
  --percTaxa 100 \
  --percIncorpMean 100 \
  --percIncorpSD 0 \
  --n_reps 1 \
  > incorp.config
  
cat incorp.config<jupyter_output>[1]
    # baseline: no incorporation
    treatment = control
    
    [[intraPopDist 1]]
        distribution = uniform
        
        [[[start]]]
            
            [[[[interPopDist 1]]]]
                distribution = uniform
                start = 0
                end = 0
        
        [[[end]]]
            
            [[[[interPopDist 1]]]]
                distribution = uniform
                start = 0
                end = 0
[2]
    # 'treatment' community: possible incorporation
    treatment = labeled
    max_perc_taxa_incorp = 100
    
    [[intraPopDist 1]]
        distribution = normal
        [[[mu]]]
            [[[[interPopDist 1]]]]
                start = 100
                distribution = uniform
                end = 100
        [[[sigma]]]
            [[[[interPopDist 1]]]]
                start = 0
                distribution = uniform
                end = 0<jupyter_text>Make file identifying incorporators for each library<jupyter_code>%%bash
echo -e "1\n2\tSynechococcus_sp_WH7803\n" > incorp_taxa.txt

cat incorp_taxa.txt<jupyter_output>1
2	Synechococcus_sp_WH7803<jupyter_text>Add isoptope incorporation<jupyter_code>!SIPSim isotope_incorp \
    --comm comm2.txt \
    --taxa incorp_taxa.txt \
    --isotope 15N \
    --np 20 \
    Wawrik_SIPsim_kde_dif_DBL.pkl \
    incorp.config \
    > Wawrik_SIPsim_kde_dif_DBL_inc.pkl

!ls -thlc Wawrik_SIPsim_kde_dif_DBL_inc.pkl 
!cat BD-shift_stats.txt<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
Loading KDE object...
Processing library: 1
Processing: Synechococcus_sp_WH7803
Processing library: 2
WARNING: config library 2 not found in KDEs.Using a different KDE object
Processing: Synechococcus_sp_WH7803
File written: BD-shift_stats.txt
-rw-rw-r-- 1 sam sam 7.7M May 22 13:22 Wawrik_SIPsim_kde_dif_DBL_inc.pkl
library	taxon	min	q25	mean	median	q75	max
1	Synechococcus_sp_WH7803	0.0	0.0	0.0	0.0	0.0	0.0
2	Synechococcus_sp_WH7803	0.0159999999999	0.016	0.016	0.016	0.016	0.016<jupyter_text>6) Making OTU table Get fraction table for each library<jupyter_code>!cat fracs.txt
!SIPSim OTU_table \
    --np 1 \
    --abs 1e5 \
    Wawrik_SIPsim_kde_dif_DBL_inc.pkl \
    comm2.txt \
    fracs.txt \
    > Wawrik_SIPsim_OTU.txt
    
!cat Wawrik_SIPsim_OTU.txt<jupyter_output>DEBUG:CACHEDIR=/home/sam/.cache/matplotlib
DEBUG:Using fontManager instance from /home/sam/.cache/matplotlib/fontList.json
DEBUG:backend agg version v2.2
Loading files...
Simulating OTUs...
Processing library: "1"
  Processing taxon: "Synechococcus_sp_WH7803"
   taxon abs-abundance:  100000
Processing library: "2"
  Processing taxon: "Synechococcus_sp_WH7803"
   taxon abs-abundance:  100000
library	taxon	fraction	BD_min	BD_mid	BD_max	count	rel_abund
1	Synechococcus_sp_WH7803	-inf-1.687	-inf	1.686	1.686	29	1
1	Synechococcus_sp_WH7803	1.687-1.693	1.687	1.69	1.693	25	1
1	Synechococcus_sp_WH7803	1.693-1.698	1.693	1.696	1.698	119	1
1	Synechococcus_sp_WH7803	1.698-1.700	1.698	1.699	1.7	211	1
1	Synechococcus_sp_WH7803	1.700-1.704	1.7	1.702	1.704	628	1
1	Synechococcus_sp_WH7803	1.704-1.708	1.704	1.706	1.708	1547	1
1	Synechococcus_sp_WH7803	1.708-1.712	1.708	1.71	1.712	3882	1
1	Synechococcus_sp_WH7803	1.712-1.716	1.712	1.714	1.716	17619	1
1	Synechococcus_sp_WH7803	1.716-1.721	1.716	1.719	1.721	[...]<jupyter_text>Part 3: AnalysisNow I'll see how well the simulation results correspond with the actual DNA distribution. Here we will also make easier to handle files for nice visualizationg in R. Get Wawrik dataAlso make a dataframe that indicates teh isotope for each reference in each library. This gives a bit more metadata that you can use when plotting<jupyter_code>%%R
real.df = read.table('/home/sam/data/SIPSim2_data/dev/validation/figure_estimates.txt', 
                          header=TRUE, sep="\t") %>%
    filter(Study == "Wawrik2009") %>%
    mutate(BD = round(BD, digits=3), data_type = "empirical") %>%
    select(taxon_name, Isotope, library, BD, Ratio_of_maximum, data_type, Study)

real.df
%%R
metadata = real.df %>%
    select(taxon_name, Isotope, library, Study) %>%
    unique
metadata<jupyter_output><empty_output><jupyter_text>metaSIPSim dataGet the metaSIPSim simulated data and estimate the concentration of DNA from each reference as a ratio to the maximum concentration.<jupyter_code>frag_table = pd.DataFrame(columns=['taxon_name', 'library', 'OriBD', 
                                   'percent_incorp', 'abundance',
                                   'fragment_start', 'fragment_length'])
for file in os.listdir(workdir):
    if file.endswith('fragments.txt.gz'):
        fraction = file.split('_')[3]
        frag_table_sub = pd.read_table(file, compression='gzip')
        frag_table_sub['fraction'] = fraction
        frag_table = frag_table.append(frag_table_sub, ignore_index=True)
        
frag_table.to_csv('allfragments.txt', sep='\t', index=False)
frag_table = None
%%R
frag.df = read.delim('allfragments.txt') %>%
    filter(abundance > 0)

head(frag.df)
%%R
taxa.df = frag.df %>%
    mutate(base_abundance = abundance*fragment_length) %>%
    group_by(taxon_name, fraction, library) %>%
    summarize(DNA_abundance = sum(base_abundance)) %>%
    as.data.frame
frag.df = NULL
taxa.df
%%R
fracs.df = read.delim('fracs.txt') %>%
    mutate(BD = (BD_min + BD_max)/2) %>%
    select(library, fraction, BD)
taxa.df = inner_join(taxa.df, fracs.df)
taxa.df
%%R
max.abds = taxa.df %>%
    group_by(library, taxon_name) %>%
    summarize(max_DNA = max(DNA_abundance)) %>%
    as.data.frame

taxa.df = inner_join(taxa.df, max.abds) %>%
    mutate(Ratio_of_maximum = DNA_abundance/max_DNA)
taxa.df
%%R
## Make dataframe that you can combine with original leuders data as well as SIPSim data. Write this to a file.
metaSIPSim.df = taxa.df %>%
    left_join(metadata, by=c("taxon_name", "library")) %>%
    mutate(data_type = "metaSIPSim") %>%
    select(taxon_name, Isotope, library, BD, Ratio_of_maximum, data_type, Study)

write.table(metaSIPSim.df, file="/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/sim_DNA_conc_metaSIPSim.txt", 
            quote=FALSE, sep="\t", col.names=TRUE, row.names=FALSE)

metaSIPSim.df<jupyter_output><empty_output><jupyter_text>Get SIPSim dataGet the SIPSim simulated data (OTU table) and estimate the DNA abundance for each reference as a ratio to the maximum abundance<jupyter_code>%%R
SIPSim.df = read.delim('Wawrik_SIPsim_OTU.txt') %>%
    filter(BD_min != "-Inf", BD_max != "Inf") %>%
    select(library, taxon, BD_mid, count) %>%
    rename(BD=BD_mid, taxon_name = taxon)

SIPSim.df
%%R
max_count = SIPSim.df %>%
    group_by(library, taxon_name) %>%
    summarize(max_count = max(count)) %>%
    as.data.frame

SIPSim.df = inner_join(SIPSim.df, max_count) %>%
    mutate(Ratio_of_maximum = count/max_count) %>%
    inner_join(metadata, by=c("taxon_name", "library")) %>%
    mutate(data_type = "SIPSim") %>%
    select(taxon_name, Isotope, library, BD, Ratio_of_maximum, data_type, Study)

write.table(SIPSim.df, file="/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/sim_DNA_conc_SIPSim.txt", 
            quote=FALSE, sep="\t", col.names=TRUE, row.names=FALSE)

SIPSim.df<jupyter_output>/home/sam/anaconda3/envs/SIPSim2_py2/lib/python2.7/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: Joining, by = c("library", "taxon_name")

  warnings.warn(x, RRuntimeWarning)<jupyter_text>Combine data into one dataframe<jupyter_code>%%R
full.df = rbind(metaSIPSim.df, SIPSim.df, real.df)
write.table(full.df, file="/home/sam/data/SIPSim2_data/dev/validation/Wawrik2009/DNA_conc_all_sets.txt", 
            quote=FALSE, sep="\t", col.names=TRUE, row.names=FALSE)

head(full.df)<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>%%R -h 350 -w 750

SIPSim2.plot = ggplot(data=metaSIPSim.df, aes(x=BD, y=Ratio_of_maximum, color=Isotope)) +
    geom_point() +
    geom_line() +
    xlim(1.688, 1.76)
SIPSim2.plot
%%R -h 350 -w 750

SIPSim.plot = ggplot(data=SIPSim.df, aes(x=BD, y=Ratio_of_maximum, color=Isotope)) +
    geom_point() +
    geom_line() +
    xlim(1.688, 1.76)
SIPSim.plot
%%R -h 350 -w 750
real.plot = ggplot(data=real.df, aes(x=BD, y=Ratio_of_maximum, color=Isotope)) +
    geom_point() +
    geom_line() +
    xlim(1.688, 1.76)

real.plot
%%R -h 350 -w 750
full.plot = ggplot(data=full.df, aes(x=BD, y=Ratio_of_maximum, color=Isotope, 
                                     shape=data_type, linetype=data_type)) +
    geom_point() +
    geom_line() +
    xlim(1.688, 1.76)

full.plot<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1265868533995934, 'fraction_numerical': 0.07973797153828778, 'mean_word_length': 4.115784608273631, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 143, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '23', 'dir': 'jupyter-structured-clean-dedup', 'id': '429012', 'license': 'MIT', 'n_tokens_mistral': 9920, 'n_tokens_neox': 8611, 'n_words': 1802, 'path': 'validation/Validation_Wawrik2009.ipynb', 'repo_name': 'seb369/MetaSIPSim', 'size': '22178'}