<jupyter_start><jupyter_text>TODO:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4750071/pdf/srep20948.pdf suggests about 75% on-target enrichment- 25% depletion of molecules not in bait library- 25% reduced efficiency of capture  Setup<jupyter_code>import pandas as pd
from pandas.api.types import CategoricalDtype, is_categorical_dtype
import numpy as np
from plotnine import *
import plotnine
import matplotlib as mp
import scipy.stats as ss
import warnings
warnings.filterwarnings('ignore')

def convert_params(mu, theta):
    """
    Convert mean/dispersion parameterization of a negative binomial to the ones scipy supports

    See https://en.wikipedia.org/wiki/Negative_binomial_distribution#Alternative_formulations
    """
    var = mu + theta * mu ** 2
    p = (var - mu) / var
    return theta, 1 - p<jupyter_output><empty_output><jupyter_text>Background and motivationMetagenomics is important because XXX.Sequencing sensitivity/depth is inherently proportional to input DNA species abundance.Therefore, detection of rare DNA species requires expensive, ultra-high depth sequencing. ObjectiveTo reduce cost and increase sensitivity, I propose using a simple based-based depletion strategy to remove high-abundance DNA species from a metagenomic NGS library. ApproachThis approach was designed to be agnostic of the input DNA content uses completed NGS libraries as input; thus, requiring minimal changes to molecular biology workflows.According to Castle et al. (10.1094/PBIOMES-09-17-0037-R), a 1:10 dilution followed by amplification of a soil fungal metagenome library resulted in approximately 20% loss of diversity. Their analysis showed that that low-abundance species in an undiluted control were the depleted constituents during the dilution step.DNA capture probes created from this biased library would reflect the DNA abundance distribution of the input metagenome excluding ca. 20% of low-abundance species. Therefore, upon depleting the input library with above capture probes, the resultant library will be depleted of high-abundance species.Sequencing of this depleted library will be more cost-efficient and sensitive because reads will be allocated to rare species.<jupyter_code># ingest published data
supp_df = pd.read_csv('s3://daniel.le-work/deplete_seq/Castle_etal_2018_qpcr.csv')
supp_df['dil_factor'] = supp_df['dil_factor'].astype(str)

# boxplots of observed OTUs
print(ggplot(supp_df)
      + theme_bw()
      + geom_boxplot(aes('dil_factor','observed OTUs', color = 'dilution'))
      + labs(x='dilution factor')
     )

# table of observed OTUs
ratio_df = (supp_df
            .drop('molecules (measured or estimated)', axis=1)
            .groupby(['dil_factor','dilution'])
            .median()
            .reset_index()
            .pivot(index = 'dil_factor', columns='dilution', values = 'observed OTUs'))
ratio_df['fraction_depleted'] = 1 - (ratio_df.Dilute / ratio_df.Undil)
display(ratio_df)<jupyter_output><empty_output><jupyter_text>Proposed methodAn adaptation of Roche NimbleGen SeqCap EZ protocol: 10.1002/0471142727.mb0721s112 Step 1: DiluteThe target NGS library (conc. TBD) is sampled and placed into dilution buffer, resulting in 1:>10 final dilution.  Step 2: Amplify and biotinylateAmplify diluted library with R1/R2 primers that contain 5'-biotin tag.   Step 3: CaptureHybridize capture probes to input solution: melt and slow anneal. Step 4: DepleteAdd magnetic steptavidin-coated beads to capture pool. Supernatant = final library.  Example empirical species distributionsData from E. Crawford's group: s3://alyden-bucket/AIH-Reseq/project-aihreseq_host-gene-counts.tar.gzs3://czb-seqbot/fastqs/190211_NB501961_0210_AHGWLVBGX9/rawdata/RR0026_002_SSa_S2_R1_001.fastq.gzThese data suggest species frequency follows a Poisson process.<jupyter_code># Dashed dataset
# s3://alyden-bucket/AIH-Reseq/project-aihreseq_host-gene-counts.tar.gz

fn1 = 'cut-filt-85-98-90_DASHL_AASLD-067_AIHReseq_DASHed_S212_L003_16157'
fn2 = 'cut-filt-85-98-90_AASLD-067_AIH_1302019_S276_L004_16164'
trial = 's1'

lib_df = pd.DataFrame()
for fn in [fn1, fn2]:
    if fn == fn1:
        group = 'dashed'
    else:
        group = 'not_dashed'
    tmp = pd.read_csv(f'/home/ubuntu/data/DL20190322_deplete_seq/{fn}', 
                           skiprows = [0,1,2,3],
                           sep = '\t',
                           header = None
                          )
    tmp.columns = ['species', 's1', 's2', 's3']
    lib_df[group] = tmp[trial].tolist()
lib_df['species'] = [f's_{x}' for x in range(len(lib_df))]

lib_df = lib_df.sample(200).sort_values('not_dashed', ascending = False)
plot_df = pd.melt(lib_df, id_vars='species')
plot_df['species'] = (plot_df['species']
                      .astype(str)
                      .astype(CategoricalDtype(lib_df.species.tolist(), ordered = True))
                     )

plotnine.options.figure_size = (3,3)
print(ggplot(plot_df)
      + theme_bw()
        + theme(axis_text_x = element_blank(),
          axis_ticks = element_blank(),
          panel_grid = element_blank()
         )
      + geom_point(aes('species', 'value', color = 'variable'))
#       + scale_y_log10()
      +labs(y='counts')
     )

# pig-pen DNA
# s3://czb-seqbot/fastqs/190211_NB501961_0210_AHGWLVBGX9/rawdata/RR0026_002_SSa_S2_R1_001.fastq.gz
## 27-mer table via KAT: k-mer frequency used as proxy for species frequency

kmer = pd.read_csv('/GB100_1/deplete_seq/test', skiprows=[x for x in range(6)],
                   header = None,
                   sep = ' '
                  )
kmer.columns = ['species', 'count']

plotnine.options.figure_size = (3,3)
print(ggplot(kmer.head(100))
      + theme_bw()
        + theme(axis_text_x = element_blank(),
          axis_ticks = element_blank(),
          panel_grid = element_blank()
         )
      + geom_point(aes('species', 'count'))
      + scale_y_log10()
      + labs(y = 'log10(counts)')
     )<jupyter_output><empty_output><jupyter_text>Simulated results  Simulated inputAssume negative binomial distribution of species abundance.  Model self vs non-self hybrid capture processAssume uniform probability of self vs non-self hybridization, proportional to the stoichiometry of self to non-self molecules.   Simulated outputAssume uniform reduction of respective sampling probabilities such that the proportion of depleted species is equal to that reported in above literature for a given dilution factor.Assume uniform capture efficient for a given input:probe stoichiometry.<jupyter_code># assume negative binomial expression level distribution
n, p = convert_params(30,0.5)
species = np.arange(ss.nbinom.ppf(0.001, n, p),ss.nbinom.ppf(0.999, n, p)) # Percent point function (inverse of cdf — percentiles)
species = species - np.min(species)
rv = ss.nbinom(n, p) # Random variates
prob = rv.pmf(species) # Probability mass function (function that describes probs for discrete variables)
prob = prob/sum(prob)
sample_size = 100000

# original distribution
raw_counts_df = pd.DataFrame(
                              pd.Series(
                                  np.random.choice(species, 
                                                   size = sample_size, 
                                                   p = prob)
                              )).rename(columns = {0:'species'})

# probs by species
theoretical_df = pd.DataFrame({'original_prob':prob,
                               'species':species
                              })
theoretical_df['species'] = theoretical_df['species'].astype(int)

# counts by species
theoretical_df = pd.merge(theoretical_df,(raw_counts_df['species']
                                          .value_counts()
                                          .reset_index()
                                          .rename(columns = {'species':'original_counts',
                                                             'index':'species'})
                                         ),
                          'left',
                          'species'
                         )
theoretical_df['original_counts'] = theoretical_df['original_counts'].replace(np.nan, 0)

plotnine.options.figure_size = (3,3)
print(ggplot()
      + theme_bw()
      + geom_histogram(raw_counts_df, aes('species', '..density..'), bins = 30)
      + geom_line(theoretical_df, aes('species', 'original_prob'), color = 'red')
      + labs(x = 'species', y = 'density')
      + ggtitle(f'Simulated library species distribution (sampling n = {sample_size})')
     )


# Model self vs non-self hybrid capture process
stoichiometry_df = pd.DataFrame({'self': [1]*19,
                                 'non_self':[x/10 for x in range(1,11,1)]+[x for x in range(2,11,1)]
                                })
stoichiometry_df['input_input'] = stoichiometry_df.self/(stoichiometry_df.self + stoichiometry_df.non_self)
stoichiometry_df['input_probe'] = 1-stoichiometry_df['input_input']
stoichiometry_df = pd.melt(stoichiometry_df, id_vars=['self', 'non_self'])

plotnine.options.figure_size = (3,3)
print(ggplot(stoichiometry_df)
      + theme_bw()
      + geom_line(aes('non_self','value', color = 'variable'))
      + geom_point(aes('non_self','value', color = 'variable'))
      + labs(x = '[probe]:[input] ratio', y = 'hybridization probability')
      + ylim((0,1))
     )

############# output #############
dilution_factor = 100
stoichiometry_factor = 3
nonspecific_fraction = 0.25

# depleted population based on expectation from literature (ie modeling dilution effect)
depletion_fraction = ratio_df[[x == str(dilution_factor) for x in ratio_df.index]]['fraction_depleted'].values[0]
theoretical_df = theoretical_df.sort_values('original_prob')
theoretical_df['keep_species'] = True
theoretical_df.loc[:int(len(theoretical_df) * (1-depletion_fraction)),'keep_species'] = False

# capture library
### uniform decrease in observing species in capture library to model drop out
uniform_loss = theoretical_df[theoretical_df.keep_species == False]['original_prob'].max()
theoretical_df['capature_prob'] = theoretical_df['original_prob'] - uniform_loss
### 0 min
theoretical_df['capature_prob'] = [x if x > 0 else 0 for x in theoretical_df['capature_prob']]
### normalize
theoretical_df['capature_prob'] = theoretical_df['capature_prob']/theoretical_df['capature_prob'].sum()

# depleted library
### depletion
theoretical_df['enriched_prob'] = theoretical_df['original_prob'] - theoretical_df['capature_prob']
### 0 min
theoretical_df['enriched_prob'] = [x if x > 0 else 0 for x in theoretical_df['enriched_prob']]
### normalize
theoretical_df['enriched_prob'] = theoretical_df['enriched_prob'] / theoretical_df['enriched_prob'].sum()
### uniform capture inefficiency 
uniform_efficiency = 1-stoichiometry_df[[x=='input_probe' and y==stoichiometry_factor for x,y in zip(stoichiometry_df.variable, stoichiometry_df.non_self)]]['value'].values[0]
theoretical_df['enriched_prob'] = theoretical_df['enriched_prob'] + (theoretical_df['original_prob']*uniform_efficiency)
## normalize for sampling
theoretical_df['enriched_prob'] = theoretical_df['enriched_prob'] / theoretical_df['enriched_prob'].sum()

# sampling
theoretical_df = theoretical_df.sort_values('species')
theoretical_df = pd.merge(theoretical_df,
                          pd.DataFrame(
                              pd.Series(
                                  np.random.choice(species, 
                                                   size = sample_size - int(nonspecific_fraction*sample_size),
                                                   p = theoretical_df['enriched_prob'].values)
                              ).value_counts()
                          )
                          .reset_index()
                          .rename(columns = {0:'enriched_counts',
                                             'index':'species'}),
                          'left',
                          'species'
                         )
theoretical_df['enriched_counts'] = theoretical_df['enriched_counts'].replace(np.nan, 0)


## capture nonspecificity (ref: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4750071/pdf/srep20948.pdf)
theoretical_df['enriched_counts'] = [int(y*nonspecific_fraction) + x for x,y in zip(theoretical_df['enriched_counts'],theoretical_df['original_counts'])]

# Plot showing enrichment vs original counts
plot_df = theoretical_df.sort_values('original_counts', ascending = False)
plot_df['species_cat'] = (plot_df['species']
                          .astype(str)
                          .astype(CategoricalDtype(plot_df['species'].astype(str).values.tolist(), ordered=True))
                         )
plot_df = pd.melt(plot_df, value_vars=['enriched_counts','original_counts'], id_vars=['species_cat'])

plotnine.options.figure_size = (3,3)
print(ggplot(plot_df)
      + theme_bw()
      + theme(axis_text_x = element_blank(),
              axis_ticks = element_blank(),
              panel_grid = element_blank()
             )
      + geom_point(aes('species_cat', 'value', color = 'variable'))
      + labs(x='species (ordered by abundance)', y='log10(frequency)')
      + scale_y_log10(breaks = [1,10,100,1000,10000])
      + ggtitle(f'Simulated enrichment distribution\n(sampling n = {sample_size}, dilution factor = {dilution_factor}, input:probe = 1:{stoichiometry_factor}, non-specific fraction = {nonspecific_fraction})')
     )

print(ggplot(plot_df)
      + theme_bw()
      + theme(axis_text_x = element_blank(),
              axis_ticks = element_blank(),
              panel_grid = element_blank()
             )
      + geom_point(aes('species_cat', 'value', color = 'variable'))
      + labs(x='species (ordered by abundance)', y='frequency')
#       + scale_y_log10()
      + ggtitle(f'Simulated enrichment distribution\n(sampling n = {sample_size}, dilution factor = {dilution_factor}, input:probe = 1:{stoichiometry_factor}, non-specific fraction = {nonspecific_fraction})')
     )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11236840220582969, 'fraction_numerical': 0.023777125259614697, 'mean_word_length': 3.146080760095012, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 20, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '151941', 'license': 'MIT', 'n_tokens_mistral': 4571, 'n_tokens_neox': 4075, 'n_words': 953, 'path': 'DL20190308_deplete_seq_simulation.ipynb', 'repo_name': 'czbiohub/biohub_scratch1', 'size': '13963'}