<jupyter_start><jupyter_text>Estimating the number of 16S rRNA genes per ng of DNA* Used to estimate the number of 16S genes in a CsCl gradient   * Used ~5 ug of DNA per gradient   Calculations* Number of genomes in X DNA  * MW of DNA  * Size distribution of bacterial genomes* Number of 16S gene copies in X DNA  * Distribution of 16S gene copy number per genome  Setting variables<jupyter_code>import os

workDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/SSU_genes_per_ng_DNA/'
rnammerDir = os.path.join(workDir + 'rnammer')
genomeDir = '/home/nick/notebook/SIPSim/dev/bac_genome1210/genomes/'<jupyter_output><empty_output><jupyter_text>Init<jupyter_code>import glob
import pyfasta
import numpy as np
import pandas as pd
from collections import Counter
import matplotlib.pyplot as plt
import scipy.stats as ss
from fitter import Fitter
from functools import partial

%matplotlib inline  
%load_ext rpy2.ipython
%%R
library(dplyr) 
library(tidyr) 
library(ggplot2) 
if not os.path.isdir(workDir):
    os.makedirs(workDir)
    
if not os.path.isdir(rnammerDir):
    os.makedirs(rnammerDir)<jupyter_output><empty_output><jupyter_text>Size distribution of bacterial genomes<jupyter_code>p = os.path.join(genomeDir, '*.fasta')
genomeFiles = glob.glob(p)

print 'Number of genome files: {}'.format(len(genomeFiles))<jupyter_output>Number of genome files: 1210<jupyter_text>Distribution of 16S gene copies per genome<jupyter_code>total_seq_len = lambda x: sum([len(y) for y in x.values()])

def total_genome_lens(genome_files):
    genome_lens = {}
    for fasta in genome_files:
        name = os.path.split(fasta)[-1]
        name = os.path.splitext(name)[0]

        pyf = pyfasta.Fasta(fasta)
        genome_lens[name] = [total_seq_len(pyf)]
    return genome_lens

genome_lens = total_genome_lens(genomeFiles)
df_genome_len = pd.DataFrame(genome_lens).transpose()
df_genome_len
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(df.ix[:,0], bins=20)<jupyter_output><empty_output><jupyter_text>Fitting distribution<jupyter_code>fo = Fitter(df_genome_len.ix[:,0])
fo.fit()
fo.summary()
genome_len_best_fit = fo.fitted_param['rayleigh']
genome_len_best_fit
# test of distribution
x = ss.rayleigh.rvs(*genome_len_best_fit, size=10000)
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(x, bins=50)
fig.show()<jupyter_output><empty_output><jupyter_text>Distribution of 16S gene copies per genome  rnammer run<jupyter_code>%%bash -s "$genomeDir" "$rnammerDir"

find $1 -name "*fasta" | \
    perl -pe 's/.+\/|\.fasta//g' | \
    xargs -n 1 -I % -P 30 bash -c \
    "rnammer -S bac -m ssu -gff $2/%_rrn.gff -f $2/%_rrn.fna -xml $2/%_rrn.xml < $1/%.fasta"
## Summarizing the results

!cd $rnammerDir; \
    egrep -v "^#" *.gff | \
    grep "16s_rRNA" | \
    perl -pe 's/:/\t/' > ssu_summary.txt

inFile = os.path.join(rnammerDir, 'ssu_summary.txt')

inFH = open(inFile, 'rb')
df_ssu = pd.read_csv(inFH, sep='\t', header=None)
df_ssu.head()
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(df_ssu.ix[:,6], bins=50)
fig.show()
# filtering by gene length of >= 1000 bp
df_ssu_f = df_ssu.loc[df[6] >= 1000]
df_ssu_f.head()
# counting number of 16S genes per genome
ssu_count = Counter(df_ssu_f[1])
ssu_max = max(ssu_count.values())

# plotting distribution
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(ssu_count.values(), bins=ssu_max)
fig.show()<jupyter_output><empty_output><jupyter_text>Fitting distribution<jupyter_code>fo = Fitter(ssu_count.values())
fo.fit()
fo.summary()
ssu_ray_fit = fo.fitted_param['rayleigh']
ssu_ray_fit
# test of distribution
x = ss.rayleigh.rvs(*ssu_ray_fit, size=10000)
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(x, bins=50)
fig.show()
ssu_beta_fit = fo.fitted_param['beta']
ssu_beta_fit
# test of distribution
x = ss.beta.rvs(*ssu_beta_fit, size=10000)
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(x, bins=50)
fig.show()<jupyter_output><empty_output><jupyter_text>Notes* Using rayleigh distribution  Monte Carlo estimation of 16S gene copies per ng of DNA* M.W. of dsDNA = ( nucleotides x 607.4) + 157.9<jupyter_code># example of calculations
gradient_DNA_conc = 1e-9    # g of DNA
avogadro = 6.022e23         # molecules/mole

genome_len = 4000000
mw_genome = genome_len * 607.4 + 157.9

n_genomes = gradient_DNA_conc / mw_genome * avogadro

ssu_copy_per_genome = 4

n_genomes * ssu_copy_per_genome
def SSU_copies_in_ng_DNA(DNA_conc, genome_len, ssu_copy_per_genome):
    
    DNA_conc__g = DNA_conc * 1e-9     # ng --> g of DNA
    avogadros = 6.022e23              # molecules/mole

    mw_genome = genome_len * 607.4 + 157.9
    n_genomes = DNA_conc__g / mw_genome * avogadros
    ssu_copies = n_genomes * ssu_copy_per_genome
    
    return ssu_copies

# run
SSU_copies_in_ng_DNA(1, 4000000, 4)
def SSU_copies_MC(DNA_conc, genome_len_dist, ssu_copy_dist, n=100000):
    n_copy_dist = []
    for i in range(n):
        genome_len = genome_len_dist(size=1)[0]
        ssu_copy_per_genome = ssu_copy_dist(size=1)[0]
        n_copies = SSU_copies_in_ng_DNA(DNA_conc, genome_len, ssu_copy_per_genome)
        n_copy_dist.append(n_copies)                     
    return n_copy_dist
# distribution functions
genome_len_dist = partial(ss.rayleigh.rvs, *genome_len_best_fit)
ssu_copy_dist = partial(ss.rayleigh.rvs, *ssu_ray_fit)

# monte carlo estimation of ssu copies in a gradient
gradient_dna_conc__ng = 5000
n_copy_dist = SSU_copies_MC(gradient_dna_conc__ng, genome_len_dist, ssu_copy_dist, n=10000)
fig = plt.figure()
ax = plt.subplot(111)
ax.hist(n_copy_dist, bins=50)
fig.show()
median_copy = int(np.median(n_copy_dist))
std_copy = int(np.std(n_copy_dist))

print 'Number of SSU copies in {} ng of DNA: {} +/- {}'.format(gradient_dna_conc__ng, median_copy, std_copy)
def median_confidence_interval(data, confidence=0.95):
    a = 1.0*np.array(data)
    n = len(a)
    m, se = np.median(a), ss.sem(a)
    h = se * ss.t._ppf((1+confidence)/2., n-1)
    return m, m-h, m+h
mci = median_confidence_interval(n_copy_dist)
mci = map(int, mci)

# lci,hci = ss.norm.interval(0.05, loc=np.mean(n_copy_dist), scale=np.std(n_copy_dist))
# copy_median = np.median(n_copy_dist)
# mci = [copy_median, copy_median - lci, copy_median + hci]

print 'Number of SSU copies in {} ng of DNA: {:,d} (low:{:,d}, high:{:,d})'.format(gradient_dna_conc__ng, *mci)<jupyter_output>Number of SSU copies in 5000 ng of DNA: 5,097,287,217 (low:4,810,572,731, high:5,384,001,703)
{'fraction_non_alphanumeric': 0.12734668335419275, 'fraction_numerical': 0.03520025031289111, 'mean_word_length': 4.708035714285714, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 2}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '574095', 'license': 'MIT', 'n_tokens_mistral': 2825, 'n_tokens_neox': 2554, 'n_words': 539, 'path': 'ipynb/bac_genome/SSU_genes_per_ng_DNA.ipynb', 'repo_name': 'arischwartz/test', 'size': '6392'}