<jupyter_start><jupyter_text>Sanger - Sub-samplingThis Jupyter notebook reproduces the results from the B-ALL sub-sampling analysis (Supplementary Figure S5).<jupyter_code>%reload_ext autoreload
%autoreload 2

%matplotlib inline

import sys
sys.path.append('../src')

import numpy as np
import pandas as pd
import seaborn as sns

from matplotlib import pyplot as plt

sns.set_style('white')<jupyter_output><empty_output><jupyter_text>First, we read the insertions and CTGs that were identified at the different sequencing depths by IM-Fusion. These results have been pre-computed using an external Snakemake pipeline (see the documentation and the Makefile for more details).<jupyter_code>depths = [15, 30, 50, 70]

# Read insertions.
insertions = pd.concat((pd.read_csv('../data/processed/sanger/star-subsample/'
                                    '{}/insertions.txt'.format(depth), sep='\t')
                          .query('support >= 10')
                          .assign(depth=depth)
                        for depth in depths), axis=0)

# Read CTGs.
ctgs = pd.concat((pd.read_csv('../data/processed/sanger/star-subsample/'
                              '{}/ctgs.txt'.format(depth), sep='\t')
                    .assign(depth=depth)
                  for depth in depths), axis=0)<jupyter_output><empty_output><jupyter_text>Next, we plot the number of insertions and DE CTGs across the different sequencing depths to see how depth affects insertion/CTG detection:<jupyter_code>def plot_depth_ctg_count(ctgs, ax=None):
    if ax is None:
        _, ax = plt.subplots()
        
    counts = (ctgs.groupby('depth')['gene_name'].nunique()
                  .reset_index(name='count'))
    sns.barplot(data=counts, x='depth', y='count', ax=ax)
    
    ax.set_xlabel('Sequencing depth')
    ax.set_ylabel('Number of CTGs')

    ax.set_ylim(0, counts['count'].max())
    
    return ax


def plot_depth_insertion_count(insertions, ax=None, kind='bar'):
    if ax is None:
        _, ax = plt.subplots()
        
    counts = (insertions.groupby('depth')['id'].nunique()
                        .reset_index(name='count'))
    
    if kind == 'bar':
        sns.barplot(data=counts, x='depth', y='count', ax=ax)
    elif kind == 'reg':
        sns.regplot(data=counts, x='depth', y='count', ax=ax)
    else:
        raise ValueError('Unknown value for kind')
    
    ax.set_xlabel('Sequencing depth')
    ax.set_ylabel('Number of insertions')
    
    return ax


def plot_depth_overview(insertions, ctgs, axes=None):
    if axes is None:
        _, axes = plt.subplots(figsize=(9, 4), ncols=2, nrows=1)

    plot_depth_ctg_count(ctgs.query('de_pvalue < 0.05'), ax=axes[1])
    plot_depth_insertion_count(insertions, ax=axes[0], kind='reg')
    
    axes[0].set_title('Insertions')
    axes[1].set_title('DE CTGs')
    axes[1].set_ylabel('Number of DE CTGs')
    
    sns.despine()
    plt.tight_layout()
    
    return axes
    

fig, axes = plt.subplots(figsize=(9, 4), ncols=2)
plot_depth_overview(insertions, ctgs, axes=axes);


with sns.plotting_context('paper', font_scale=0.7):
    fig, axes = plt.subplots(figsize=(6, 3), ncols=2)
    plot_depth_overview(insertions, ctgs, axes=axes);

    fig.savefig('../reports/supplemental/figures/fig_s5ab_subsampling.pdf', bbox_inches='tight')
    plt.close(fig)<jupyter_output><empty_output><jupyter_text>This shows that with increasing depth we identify linearly increasing numbers of insertions, but that the number of DE CTGs does not change strongly. To check if these CTGs are consistent across depths, we also plot a venn diagram of the identified CTGs:<jupyter_code>from nbsupport.venn import venn

ctg_genes = {depth: set(grp.query('de_pvalue <= 0.05')['gene_name'])
              for depth, grp in ctgs.groupby('depth')}

fig = venn([ctg_genes[d] for d in depths],
           names=['{} million reads'.format(i) for i in depths],
           colors=sns.color_palette())
fig.suptitle('DE CTG overlap', fontsize=13, y=0.85)

with sns.plotting_context('paper', font_scale=0.7):
    fig = venn([ctg_genes[d] for d in depths],
               names=['{} million reads'.format(i) for i in depths],
               colors=sns.color_palette(), figsize=(6, 5))
    fig.suptitle('DE CTG overlap', fontsize=7, y=0.85)

    fig.savefig('../reports/supplemental/figures/fig_s5c_subsampling_venn.pdf', bbox_inches='tight')
    plt.close(fig)
    
    
ctg_genes<jupyter_output><empty_output><jupyter_text>This confirms that one additional CTG is identified with 50/70 million reads and that 5 CTGs are identified across all depths.To determine if we find similar support (in terms of number of samples) across the different depths, we also count the number of samples with insertions in each CTG across the sequencing depths:<jupyter_code>depth_overview = pd.pivot_table(ctgs.query('de_pvalue < 0.05'), index='gene_name',
                                columns='depth', values=['n_samples'])

depth_overview<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1239967897271268, 'fraction_numerical': 0.012239165329052969, 'mean_word_length': 3.3882042253521125, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '650178', 'license': 'MIT', 'n_tokens_mistral': 1733, 'n_tokens_neox': 1612, 'n_words': 403, 'path': 'notebooks/4. Sanger - Subsampling.ipynb', 'repo_name': 'jrderuiter/imfusion-analyses', 'size': '4984'}