<jupyter_start><jupyter_text>Goal* simulating amplicon fragments for genomes in non-singleton OTUs  Setting variables<jupyter_code>import os

workDir = '/var/seq_data/ncbi_db/genome/Jan2016/ampFragsGC/'
ampFragFile = '/var/seq_data/ncbi_db/genome/Jan2016/ampFrags_KDE.pkl'
otuFile = '/var/seq_data/ncbi_db/genome/Jan2016/rnammer_aln/otusn_map_nonSingle.txt'<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>import dill
import numpy as np
import pandas as pd
%load_ext rpy2.ipython
%load_ext pushnote
%%R
library(dplyr)
library(tidyr)
library(ggplot2)
if not os.path.isdir(workDir):
    os.makedirs(workDir)
    
%cd $workDir<jupyter_output>/var/seq_data/ncbi_db/genome/Jan2016/ampFragsGC<jupyter_text>gradient params<jupyter_code># max 13C shift
max_13C_shift_in_BD = 0.036
# min BD (that we care about)
min_GC = 13.5
min_BD = min_GC/100.0 * 0.098 + 1.66
# max BD (that we care about)
max_GC = 80
max_BD = max_GC / 100.0 * 0.098 + 1.66    # 80.0% G+C
max_BD = max_BD + max_13C_shift_in_BD
## BD range of values
BD_vals = np.arange(min_BD, max_BD, 0.001)<jupyter_output><empty_output><jupyter_text>Get GC distribution info<jupyter_code>infoFile = os.path.splitext(ampFragFile)[0] + '_info.txt'
infoFile = os.path.join(workDir, os.path.split(infoFile)[1])
!SIPSim KDE_info -s $ampFragFile > $infoFile

!wc -l $infoFile
!head -n 4 $infoFile
%%R -i infoFile

df.info = read.delim(infoFile, sep='\t') %>%
    mutate(genus_ID = gsub('_.+', '', taxon_ID),
           species_ID = gsub('^([^_]+_[^_]+).+', '\\1', taxon_ID))

df.info %>% head(n=3)<jupyter_output><empty_output><jupyter_text>Combining info table with OTU tabel<jupyter_code>%%R -i otuFile

df.OTU = read.delim(otuFile, sep='\t', header=FALSE) %>%
    mutate(genome_ID = gsub('\\.fna', '', V13)) %>%
    select(genome_ID, V2) %>%
    rename('OTU_ID' = V2)

df.info.j = inner_join(df.info, df.OTU, c('taxon_ID' = 'genome_ID'))
df.OTU = NULL

df.info.j %>% head(n=3)
%%R
df.info.j.f1 = df.info.j %>% 
    filter(KDE_ID == 1) %>%
    distinct(taxon_ID, OTU_ID) %>%
    group_by(OTU_ID) %>%
    mutate(n_taxa = n()) %>%
    ungroup() %>%
    filter(n_taxa > 1)

df.info.j.f1 %>% nrow %>% print
df.info.j.f1 %>% head(n=3) %>% as.data.frame
%%R -h 4000 -w 800

df.info.j.f1$taxon_ID = reorder(df.info.j.f1$taxon_ID, df.info.j.f1$genus_ID)
df.info.j.f1$OTU_ID = reorder(df.info.j.f1$OTU_ID, -df.info.j.f1$n_taxa)

ggplot(df.info.j.f1, aes(x=taxon_ID, y=median, 
                       ymin=percentile_25, ymax=percentile_75,
                       color=species_ID)) +
    geom_linerange() +
    geom_point(size=1) +
    facet_wrap(~ OTU_ID, scales='free_x', ncol=8) +
    theme_bw() +
    theme(
        axis.text.x = element_blank(),
        legend.position='none'
    )<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup