<jupyter_start><jupyter_text>Novel BGC candidates overview<jupyter_code>import numpy as np
import pandas as pd
from matplotlib import pyplot as plt 
from IPython.core.display import display, HTML
import seaborn as sns
import sys
sys.path.append('../bgc_detection/candidates')
from candidate_regions import candidate_regions
from biosynthetic_pfams import AS_BIO_PFAM_IDS
from ipywidgets import interact
from Bio import SeqIO
import matplotlib
font = {'family' : 'Arial', 'size'   : 12}
matplotlib.rc('font', **font)
import sys

sys.path.append('../bgc_detection/similarity')
import group_scatter_plot

from group_scatter_plot import group_scatter_plot, apply_group_color
pfams = pd.read_csv('../data/pfam/Pfam-A.clans.csv').set_index('pfam_id')
pfams.head(1)
DEEPBGC_PATH = '128lstm-100pfamdim-8pfamiter-posweighted-neg-10k-fpr2'
paths = [
    '{}/candidates.csv'.format(DEEPBGC_PATH),
    'antismash-minimal/candidates.csv',
    'antismash-clusterfinder/candidates.csv',
    #'antismash-clusterfinder-default/candidates.csv',
    #'antismash-export/candidates.csv',
    #'hmm_clusterfinder-tpr80-postproc/candidates.csv'
]
labels = [
    'DeepBGC',
    'antiSMASH',
    'ClusterFinder',
    #'antiSMASH+ClusterFinder',
    #'antiSMASH export'
    #'HMM Clusterfinder tpr=80%'
]
MODEL_COLORS = {
    'DeepBGC': 'red',
    'antiSMASH': 'orange',
    'ClusterFinder': 'dodgerblue'
}
colors = [MODEL_COLORS[l] for l in labels]

BIO_COLORS = {
        'Polyketide' : 'deeppink',
        'NRP' : 'lime',
        'RiPP' : 'blue',
        'Other' : 'darkcyan',
        'Saccharide' : 'orange',
        'Terpene': 'turquoise',
        'Alkaloid' : 'purple',
        'No class' : 'grey',
    'No confident class' : 'grey',
        'Hybrids' : 'yellow',
}
ACT_COLORS = {
    'Antibacterial': 'deeppink',
    'Antifungal': 'orange',
    'Cytotoxic': 'blue',
    'Inhibitor': 'darkcyan',
    'No class' : 'grey',
    'No confident class' : 'grey',
    'Hybrids' : 'yellow'
}
SPECIES = [
    'Mycobacterium','Bacillus','Pseudomonas','Streptomyces','Salmonella',
    'Escherichia','Paenibacillus','Amycolatopsis','Staphylococcus','Clostridium',
    'Xanthomonas','Klebsiella','Frankia','Nocardia'
]
SPECIES_COLORS = {s: matplotlib.cm.tab20.colors[i] for i, s in enumerate(SPECIES)}
candidates = []
for path, label in zip(paths, labels):
    cands = pd.read_csv('../data/bacteria/candidates/{}'.format(path))
    cands['model'] = label
    candidates.append(cands)
candidates = pd.concat(candidates)
candidates.groupby('model').size()
deepbgc = pd.read_csv('../data/bacteria/candidates/{}/candidates_merged.csv'.format(DEEPBGC_PATH))
deepbgc.head()
len(deepbgc)
all_classes = (deepbgc[['Alkaloid','NRP','Other','Polyketide','RiPP','Saccharide','Terpene','Antibacterial','Cytotoxic','Inhibitor','Antifungal']] >= 0.5).astype(np.int)
all_classes.head()
fig, ax = plt.subplots(1, 1, figsize=(7, 4))
sns.heatmap(all_classes.T.dot(all_classes), annot=True, fmt='.0f', cmap='Greys', ax=ax)
def plot_regions(cands, highlight_cand=None, start=None, end=None):
    if start:
        cands = cands[cands['nucl_end'] >= start]
    if end:
        cands = cands[cands['nucl_start'] <= end]
    fig = candidate_regions(cands, colors={m: c for m, c in zip(labels, colors)}, xlim=start or 0)
    if highlight_cand is not None:
        plt.axvspan(highlight_cand['nucl_start'], highlight_cand['nucl_end'], color='green', alpha=0.2)
    return fig
def plot_candidate_region(candidate_hash, start=None, end=None):
    cand = deepbgc[deepbgc['candidate_hash']==candidate_hash].iloc[0]
    cands = candidates[candidates['contig_id']==cand['contig_id']]
    return plot_regions(cands, highlight_cand=cand, start=start, end=end)<jupyter_output><empty_output><jupyter_text>DeepBGC Novel candidates<jupyter_code>mibig = pd.read_csv('../data/mibig/mibig_properties.csv').set_index('contig_id')
mibig.head()
novel = pd.read_csv('../data/bacteria/candidates/{}/novel.csv'.format(DEEPBGC_PATH))
print(len(novel))
novel.head()
ax = novel['nucl_length'].plot.box(showfliers=False)
ax = novel[['num_all_domains','num_bio_domains']].plot.box(showfliers=False)

novel_all_classes = (novel[['Alkaloid','NRP','Other','Polyketide','RiPP','Saccharide','Terpene','Antibacterial','Cytotoxic','Inhibitor','Antifungal']] >= 0.5).astype(np.int)
novel_all_classes.head()
fig, ax = plt.subplots(1, 1, figsize=(7, 4))
sns.heatmap(novel_all_classes.T.dot(novel_all_classes), annot=True, fmt='.0f', cmap='Greys', ax=ax)<jupyter_output><empty_output><jupyter_text>Biosynthetic class distribution<jupyter_code>def plot_class_counts(classes, colors, group_hybrids=False, title='', ax=None):
    classes = classes.fillna('No confident class')
    if group_hybrids:
        classes = classes.apply(lambda c: 'Hybrids' if ';' in c else c)
    class_counts = classes.value_counts()
    c = [colors.get(c, 'white') for c in class_counts.index]
    total = class_counts.sum()
    l = ['{} ({:.0%})'.format(c,v/total) for c,v  in class_counts.items()]
    ax = class_counts.plot.pie(wedgeprops=dict(width=0.5), ax=ax, colors=c, labels=l)
    ax.set_ylabel('')
    ax.set_xlabel(title)
def plot_species_counts(species, title='', ax=None):
    species = species.apply(lambda species: species.split('_')[0])
    top_species = species.value_counts()[:10].index
    top_species_counts = species.apply(lambda s: s if s in top_species else 'Others').value_counts()
    c = [SPECIES_COLORS.get(c, 'lightgrey') for c in top_species_counts.index]
    total = top_species_counts.sum()
    l = ['{} ({:.0%})'.format(c,v/total) for c,v  in top_species_counts.items()]
    ax = top_species_counts.plot.pie(wedgeprops=dict(width=0.5), colors=c, labels=l, ax=ax)
    ax.set_ylabel('')
    ax.set_xlabel(title)<jupyter_output><empty_output><jupyter_text>Activity class distribution<jupyter_code>plot_class_counts(deepbgc['classes1'], ACT_COLORS)
plot_class_counts(novel['classes1'], ACT_COLORS)
fig, axes = plt.subplots(2, 2, figsize=(13, 8))
plot_class_counts(deepbgc['classes'], BIO_COLORS, group_hybrids=True, ax=axes[0][0])
plot_class_counts(novel['classes'], BIO_COLORS, group_hybrids=True, ax=axes[0][1])

plot_species_counts(deepbgc['species'], title='All candidates', ax=axes[1][0])
plot_species_counts(novel['species'], title='Novel candidates', ax=axes[1][1])

fig.savefig('../data/figures/paper/supplementary/10a.pdf', bbox_inches='tight')
def plot_heatmap(heatmap_cands, columns=None, **kwargs):
    heatmap = pd.DataFrame()
    top_species = species.value_counts()[:10].index
    heatmap['species'] = heatmap_cands['species'].apply(lambda species: species.split('_')[0]).apply(lambda s: s if s in top_species else 'Others')
    heatmap['classes'] = heatmap_cands['classes'].fillna('No confident class')
    heatmap = heatmap.pivot_table(values='species', index='species', columns=['classes'], aggfunc=len)
    if columns is not None:
        heatmap = heatmap.transpose().loc[columns].fillna(0).transpose()
    heatmap = heatmap.fillna(0).astype(np.int)
    heatmap = heatmap.loc[top_species.tolist() + ['Others']]
    heatmap = heatmap.apply(lambda row: row / row.sum(), axis=1)
    sns.heatmap(heatmap, annot=True, fmt='.0%', cmap='Greys', **kwargs)
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
columns = deepbgc['classes'].sort_values().fillna('No confident class').unique()
plot_heatmap(deepbgc, columns=columns, ax=axes[0], cbar=False)
plot_heatmap(novel, columns=columns, ax=axes[1], cbar=False)
axes[1].set_ylabel('')
plt.tight_layout()
fig.savefig('../data/figures/paper/supplementary/10b.pdf', bbox_inches='tight')
ax = candidates.boxplot('nucl_length', by='model', showfliers=False)
ax.set_title('')
ax.figure.suptitle('')
def plot_coverage_distribution(coverage, figsize=(10, 2), cov_lim=None, cov_step=100, ax=None):
    dist = []
    thresholds = np.arange(0, 1, 0.001)
    for threshold in thresholds:
        num_over_threshold = coverage[coverage['coverage'] > threshold].groupby('model').size()
        dist.append(num_over_threshold)
    dist = pd.DataFrame(dist, index=thresholds)

    for i, name in enumerate(dist.columns):
        color = MODEL_COLORS.get(name, 'grey')
        is_postproc = 'postproc' in name
        linestyle = '--' if is_postproc else '-'
        ax.plot(dist[name], color=color, linestyle=linestyle, lw=1)
        
    if cov_lim:
        ax.set_yticks(list(range(0, cov_lim, cov_step)))
        ax.set_ylim([0, cov_lim])
    ax.set_xlim([-0.01,1.01])
    ax.set_xticklabels(['>{:.0f}%'.format(x*100) for x in ax.get_xticks()]);
    ax.set_ylabel('# BGCs')
    ax.set_xlabel('BGC Coverage')
    plt.tight_layout()
    return fig
cf_as_coverage = pd.read_csv('../data/bacteria/candidates/antismash-clusterfinder/coverage_antismash.csv')
cf_as_coverage['model'] = 'antiSMASH'
cf_db_coverage = pd.read_csv('../data/bacteria/candidates/antismash-clusterfinder/coverage_deepbgc.csv')
cf_db_coverage['model'] = 'DeepBGC'
cf_coverage = pd.concat([cf_as_coverage, cf_db_coverage])
cf_coverage.head()
db_as_coverage = pd.read_csv('../data/bacteria/candidates/{}/coverage_antismash.csv'.format(DEEPBGC_PATH))
db_as_coverage['model'] = 'antiSMASH'
db_cf_coverage = pd.read_csv('../data/bacteria/candidates/{}/coverage_clusterfinder.csv'.format(DEEPBGC_PATH))
db_cf_coverage['model'] = 'ClusterFinder'
db_coverage = pd.concat([db_as_coverage, db_cf_coverage])
db_coverage.head()
fig, axes = plt.subplots(1, 2, figsize=(7, 2))
plot_coverage_distribution(cf_coverage, ax=axes[0])
plot_coverage_distribution(db_coverage, ax=axes[1]);<jupyter_output><empty_output><jupyter_text>Selected candidates<jupyter_code>from Bio import SeqIO
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import FeatureLocation
import matplotlib.transforms as transforms
from dna_features_viewer import BiopythonTranslator, GraphicRecord, GraphicFeature

result_path = '../data/bacteria/candidates/{}/novel.gbk'.format(DEEPBGC_PATH)
"""
Credits: DNA features viewer
From: https://github.com/Edinburgh-Genome-Foundry/DnaFeaturesViewer
"""

import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.ticker import MaxNLocator

from dna_features_viewer.tools import (change_luminosity, get_text_box, compute_features_levels,
                    bokeh_feature_patch)
from dna_features_viewer.biotools import extract_translation
from dna_features_viewer.GraphicFeature import GraphicFeature

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import FeatureLocation, SeqFeature
from Bio.Alphabet import DNAAlphabet

import pandas as pd

class CustomGraphicRecord(GraphicRecord):

    def __init__(self, **kwargs):
        super().__init__(**kwargs)

    def plot_feature(self, ax, feature, level, linewidth=1.0):
        x1, x2 = feature.start, feature.end
        if feature.open_left:
            x1 -= 1
        if feature.open_right:
            x2 += 1
        if feature.strand == -1:
            x1, x2 = x2, x1
        x1, x2 = x1 - 0.5, x2 - 0.5
        is_undirected = feature.strand not in (-1, 1)
        head_is_cut = ((feature.strand == 1 and feature.open_right) or
                       (feature.strand == -1 and feature.open_left))
        head_length = (0.001 if (is_undirected or head_is_cut) else
                       max(0.3 * feature.thickness, 5))
        
        arrowstyle = mpatches.ArrowStyle.Simple(head_width=feature.thickness,
                                                tail_width=feature.thickness,
                                                head_length=head_length)
        y = self.feature_level_width * level
        patch = mpatches.FancyArrowPatch([x1, y], [x2, y],
                                         shrinkA=0.0, shrinkB=0.0,
                                         arrowstyle=arrowstyle,
                                         facecolor=feature.color, zorder=0,
                                         linewidth=feature.linewidth)
        ax.add_patch(patch)
        return patch

cand_classes = pd.read_csv('../data/bacteria/candidates/{}/classes/class-domain_set-rf/classes.csv'.format(DEEPBGC_PATH)).set_index('contig_id')['classes']
cand_classes.head()
pfam_classes = pd.read_csv('../data/clusterfinder/supplement/pfam_classes.tsv', sep='\t')
pfam_classes = pfam_classes.drop_duplicates(subset=['Accession'], keep='first')
pfam_classes = pfam_classes.set_index('Accession')['Type']
pfam_classes.head()
PFAM_COLORS = {
    'PKS':'darkred',
    'NRPS':'darkred',
    'Terpene':'darkred',
    'Oligosaccharide':'darkred',
    'Special':'darkred',
    'Ectoine':'darkred',
    'Siderophore':'darkred',
    'Butyrolactone':'darkred',
    'Beta-lactam':'darkred',
    'Cyclodipeptide':'darkred',
    'Aminoglycoside':'darkred',
    'Indolocarbazole':'darkred',
    'Nucleoside':'darkred',
    'Phosphoglycolipid':'darkred',
    'Redox':'blueviolet',
    'Group_Transfer':'orange',
    'Regulatory':'green',
    'Transport':'blue',
    'Ribpept':'darkred',
    'Other':'grey'
}
def group_genes(domains):
    first = domains.iloc[0]
    classes = pfam_classes.reindex(domains['pfam_id'].unique()).dropna().unique()
    if not len(classes):
        classes = ['Other']
    return pd.Series(dict(
        start=first['gene_start'], 
        end=first['gene_end'],
        strand=first['gene_strand'],
        pfam_ids=','.join(domains['pfam_id']),
        classes=','.join(classes),
        first_class=classes[0]
    ))
def plot_domain_genes(cand_domains, figure_width=10, **kwargs):
    cand_genes = cand_domains.groupby('protein_id').apply(group_genes)
    features = cand_genes.apply(lambda row: GraphicFeature(
        start=int(row['start']), 
        end=int(row['end']), 
        strand=int(row['strand']), 
        label='',
        color=PFAM_COLORS.get(row['first_class'], 'white')),
        axis=1
    )
    seq_len = cand_domains['gene_end'].max()
    graphic_record = CustomGraphicRecord(sequence_length=seq_len, features=features)
    graphic_record.feature_level_width = 0
    graphic_record.features
    ax, _ = graphic_record.plot(figure_width=figure_width, **kwargs)

    patches = []
    #patches.append(mpatches.Patch(facecolor='darkred', label='Class-specific', edgecolor='black', linewidth=1))
    for label, color in PFAM_COLORS.items():
        if label not in cand_genes['first_class'].values:
            continue
        patches.append(mpatches.Patch(facecolor=color, label=label, edgecolor='black', linewidth=1))
    
    ax.legend(handles=patches, loc='best', bbox_to_anchor=(1.1, 1), frameon=False)
    return ax.figure

def plot_genes(candidate, **kwargs):
    cand_domains = pd.read_csv('../data/bacteria/domains/{}.domains.csv'.format(candidate['contig_id']))
    cand_domains = cand_domains[cand_domains['gene_start'] >= candidate['nucl_start']]
    cand_domains = cand_domains[cand_domains['gene_end'] <= candidate['nucl_end']]
    cand_domains['gene_start'] -= candidate['nucl_start']
    cand_domains['gene_end'] -= candidate['nucl_start']
    return plot_domain_genes(cand_domains, **kwargs)
PALETTE = {
    'MIBiG Polyketide' : 'deeppink',
    'MIBiG NRP' : 'lime',
    'MIBiG RiPP' : 'blue',
    'MIBiG Other' : 'darkcyan',
    'MIBiG Saccharide' : 'orange',
    'MIBiG Terpene': 'turquoise',
    'MIBiG Alkaloid' : 'purple',
    'MIBiG ?' : None,
    'Predicted Polyketide' : 'deeppink:P',
    'Predicted NRP' : 'lime:P',
    'Predicted RiPP' : 'blue:P',
    'Predicted Other' : 'darkcyan:P',
    'Predicted Saccharide' : 'orange:P',
    'Predicted Terpene': 'turquoise:P',
    'Predicted Alkaloid' : 'purple:P'
}
len(novel[novel['num_all_domains'] >= 5])
novel_long = novel[novel['num_all_domains'] >= 5]
print(len(novel_long))
novel_long.head()<jupyter_output>227<jupyter_text>TODO move tsne to fpr2<jupyter_code>tsne = pd.read_csv('../data/bacteria/candidates/novel-{}/tsne/mibig_bgcs_all.pfam2vec-mean.tsne2.csv'.format(DEEPBGC_PATH)).set_index('contig_id')
#tsne = pd.read_csv('../data/bacteria/candidates/novel-{}/tsne/mibig_bgcs_all.domain_set.tsne2.csv'.format(DEEPBGC_PATH)).set_index('contig_id')
tsne.head()
novel_samples = tsne.loc[novel_long['candidate_hash'].unique()].copy()
prefix = 'Predicted '
novel_samples['marker'] = 'P'
novel_samples['alpha'] = 0.8
novel_samples['size'] = 30
novel_samples['group'] = cand_classes.loc[novel_samples.index].fillna('(no class)').apply(lambda c: prefix+c.replace(';', ';'+prefix))
novel_samples['color'] = novel_samples['group'].apply(lambda g: apply_group_color(g, PALETTE))
novel_samples.head()
mibig_samples = tsne.loc[mibig.index].dropna().copy()
prefix = 'MIBiG '
mibig_samples['marker'] = 'o'
mibig_samples['alpha'] = 0.15
mibig_samples['size'] = 20
mibig_samples['group'] = mibig.loc[mibig_samples.index]['classes'].apply(lambda c: prefix+c.replace(';', ';'+prefix) if isinstance(c, str) else c)
mibig_samples['color'] = mibig_samples['group'].apply(lambda g: apply_group_color(g, PALETTE))
mibig_samples.head()
first_cand = novel[novel['candidate_hash']=='ffd47679ea075ec67ed6bb995467cbcb'].iloc[0]
second_cand = novel[novel['candidate_hash']=='307e7e7ef2b7256aa60260c3f20ed974'].iloc[0]
third_cand = novel[novel['candidate_hash']=='ae833e7695e2a30dbf51c4db523a46b5'].iloc[0]
first_tsne = tsne.loc[first_cand['candidate_hash']]
second_tsne = tsne.loc[second_cand['candidate_hash']]
third_tsne = tsne.loc[third_cand['candidate_hash']]
samples = pd.concat([mibig_samples, novel_samples])
samples = samples.rename(columns={'a': 'pfam2vec mean: tsne1', 'b': 'pfam2vec mean: tsne2'})
fig = group_scatter_plot(samples=samples, groupby='group', figure_size=4.5, verbose=0, legend_hybrids=False, legend_loc='outside best', linewidth=0.8)
fig.axes[0].plot(first_tsne['a'], first_tsne['b'], marker='o', markersize=20, markeredgecolor='black', markerfacecolor=(0,0,0,0), lw=0.5)
#fig.axes[0].plot(third_tsne['a'], third_tsne['b'], marker='o', markersize=20, markeredgecolor='red', markerfacecolor=(0,0,0,0), lw=0.5)
#fig.axes[0].set_xlim([None, 115])
#fig.axes[0].set_ylim([-65, 60])
fig.axes[0].set_xticks([-50,-25,0,25,50])
fig.axes[0].set_yticks([-50,-25,0,25,50])
fig.axes[0].get_legend().remove()
fig.savefig('../data/figures/paper/main/5a.pdf', bbox_inches='tight')
fig = plot_genes(novel[novel['candidate_hash']=='ffd47679ea075ec67ed6bb995467cbcb'].iloc[0], figure_width=5)
#fig.savefig('../data/figures/paper/main/5b.pdf', bbox_inches='tight')
neighbors = candidates[candidates['contig_id'] == 'CP004886.1']
neighbors = neighbors[(neighbors['nucl_start'] > 270000) & (neighbors['nucl_end'] < 450000)]
neighbors = neighbors[['nucl_start','nucl_end','model','candidate_hash']].sort_values(by='nucl_start')
neighbors
314916-306057
390102-333380
fig = plot_candidate_region('ffd47679ea075ec67ed6bb995467cbcb', start=0, end=450000)
fig.set_figwidth(4)
fig.set_figheight(0.75)
#fig.savefig('../data/figures/paper/main/5c.pdf', bbox_inches='tight')
fig = plot_genes(novel[novel['candidate_hash']=='307e7e7ef2b7256aa60260c3f20ed974'].iloc[0])
fig = plot_candidate_region('307e7e7ef2b7256aa60260c3f20ed974', start=0, end=450000)
fig.set_figwidth(4)
fig.set_figheight(0.75)
candidate = novel[novel['candidate_hash']=='ae833e7695e2a30dbf51c4db523a46b5'].iloc[0]
plot_genes(candidate);
fig = plot_candidate_region('ae833e7695e2a30dbf51c4db523a46b5', start=3800000, end=4400000)
fig.set_figwidth(4)
fig.set_figheight(0.75)
# Not used, kept for reference
def plot_record_domains(record):
    translator = BiopythonTranslator()
    graphic_record = CustomGraphicRecord(sequence=record.seq, features=[translator.translate_feature(f) for f in record.features if f.type == 'CDS'])
    graphic_record.feature_level_width = 0
    for feature in graphic_record.features:
        feature.label = ''
        feature.thickness = 20
        feature.linewidth=1
    ax, _ = graphic_record.plot(figure_width=10, annotate_inline=True)
    graphic_record = CustomGraphicRecord(sequence=record.seq, features=[translator.translate_feature(f) for f in record.features if f.type == 'PFAM_domain'])
    graphic_record.feature_level_width = 2
    for feature in graphic_record.features:
        feature.label = ''
        feature.thickness = 13
        feature.color = 'orange'
        feature.linewidth=1
    ax, _ = graphic_record.plot(annotate_inline=True, ax=ax, level_offset=2)
    ax.set_ylim([-5, None])<jupyter_output><empty_output><jupyter_text>Similar to selected candidate<jupyter_code>from sklearn.metrics.pairwise import cosine_similarity, cosine_distances
domain_set = pd.read_pickle('../data/bacteria/candidates/{}/domain_set.pkl'.format(DEEPBGC_PATH))
domain_set.head()
first_domain_set = domain_set.loc[first_cand['candidate_hash']]
sum(first_domain_set)
first_similarity = pd.Series(cosine_similarity([first_domain_set], domain_set)[0], domain_set.index).sort_values(ascending=False)
first_similarity[:10]
cand
fig, axes = plt.subplots(11, 1, figsize=(5, 5.5))
i = 0
for candidate_hash, sim in first_similarity[:11].iteritems():
    matching = deepbgc[deepbgc['candidate_hash']==candidate_hash]
    cand = matching.iloc[0]
    novelty = db_coverage.set_index('candidate_id').loc[matching['candidate_id']].groupby('model')['coverage'].max()
    print('{}: Similarity {:.2f}%:'.format(cand['species'], sim*100))
    print(novelty)
    print('CF prediction: ', matching['cf_prediction'].max())
    plot_genes(cand, figure_width=7, ax=axes[i])
    axes[i].legend([]).remove()
    axes[i].set_xticks([])
    i += 1
fig.savefig('../data/figures/paper/supplementary/11a.pdf', bbox_inches='tight')
mibig_domain_set = pd.read_pickle('../data/training/positive/mibig_bgcs_all.domain_set.pkl')
mibig_domain_set.head()
common_features = np.intersect1d(domain_set.columns, mibig_domain_set.columns)
first_similarity_mibig = pd.Series(cosine_similarity([first_domain_set.loc[common_features]], mibig_domain_set[common_features])[0], mibig_domain_set.index).sort_values(ascending=False)
first_similarity_mibig[:10]
mibig_domains = pd.read_csv('../data/training/positive/mibig_bgcs_all.csv').groupby('contig_id')
for contig_id, sim in first_similarity_mibig[:5].iteritems():
    bgc = mibig_domains.get_group(contig_id)
    print('{}: Similarity {:.2f}%:'.format(contig_id, sim*100))
    fig = plot_domain_genes(bgc, figure_width=7)
    fig.set_figheight(0.4)
    fig.axes[0].legend([]).remove()
    plt.show()
deepbgc.set_index('candidate_hash').loc[first_similarity[:10].index]<jupyter_output><empty_output><jupyter_text>Overlap Venn diagram<jupyter_code>def get_single_contig_overlap(a_cands, b_cands):
    overlaps = []
    for c, cand in a_cands.iterrows():
        # For each suitable candidate from other model
        overlap = False
        for i, other in b_cands.iterrows():
            if other['nucl_start'] < cand['nucl_end'] and other['nucl_end'] > cand['nucl_start']:
                overlap = True
                break
        overlaps.append(overlap)
    return pd.Series(overlaps, index=a_cands.index)

def get_overlap(a_cands, b_cands):
    a_grouped = a_cands.groupby('contig_id')
    b_grouped = b_cands.groupby('contig_id')

    overlaps = []
    # Get overlap separately for all contigs
    for contig_id in a_grouped.groups:
        a_contig_cands = a_grouped.get_group(contig_id)
        b_contig_cands = b_grouped.get_group(contig_id) if contig_id in b_grouped.groups else pd.DataFrame()
        overlaps.append(get_single_contig_overlap(a_contig_cands, b_contig_cands))
    return pd.concat(overlaps)
def get_multiple_overlaps(candidates):
    overlaps = []
    for a_model, a_cands in candidates.groupby('model'):
        overlap = pd.DataFrame()
        for b_model, b_cands in candidates.groupby('model'):
            if a_model == b_model:
                continue
            overlap[b_model] = get_overlap(a_cands, b_cands)        
        overlap['model'] = a_model
        overlap[a_model] = True
        overlaps.append(overlap)
    return pd.concat(overlaps)
overlaps = get_multiple_overlaps(candidates)
overlaps.head()
x = deepbgc[deepbgc['coverage'] == 0].set_index('candidate_id')
len(x)
candidates[candidates['contig_id'] == 'CP002955.1'].groupby('model').size()
fig = plot_candidate_region('002d6748dad4145b1f09ca9dc93ead44', start=500000, end=2000000)
fig.set_figwidth(7)
fig.set_figheight(0.75)
overlap_bins = overlaps[['model']].copy()
num_overlap = overlaps.sum(axis=1)
overlap_bins['All overlap'] = num_overlap == 3
overlap_bins['Novel'] = num_overlap == 1
novel_or_all = overlap_bins['All overlap'] | overlap_bins['Novel']
overlap_bins['ClusterFinder overlap'] = overlaps['ClusterFinder'] & ~novel_or_all & (overlaps['model'] != 'ClusterFinder')
overlap_bins['DeepBGC overlap'] = overlaps['DeepBGC'] & ~novel_or_all & (overlaps['model'] != 'DeepBGC')
overlap_bins['antiSMASH overlap'] = overlaps['antiSMASH'] & ~novel_or_all & (overlaps['model'] != 'antiSMASH')
overlap_bins.head()
overlap_counts = overlap_bins.groupby('model').sum().astype(np.int)
overlap_counts
bin_colors = {
    'All overlap': '#444444',
    'Novel': '#00dd44',
    'ClusterFinder overlap': 'dodgerblue',
    'antiSMASH overlap': 'orange',
    'DeepBGC overlap': 'red',
    'ClusterFinder&DeepBGC overlap': 'purple'
}
counts = overlap_counts.loc['DeepBGC'][['Novel','ClusterFinder overlap', 'All overlap', 'antiSMASH overlap']]
col = [bin_colors[c] for c in counts.index]
total = counts.sum()
perc = ['{:.1f}%'.format(c/total*100) for c in counts]
ax = counts.plot.pie(colors=col, wedgeprops=dict(width=0.5), labels=perc, figsize=(3, 2), legend=False)
ax.set_ylabel('')
ax.figure.savefig('../data/figures/paper/main/5x.pdf', bbox_inches='tight')
counts = overlap_counts.loc['ClusterFinder'][['Novel','DeepBGC overlap', 'All overlap', 'antiSMASH overlap']]
col = [bin_colors[c] for c in counts.index]
total = counts.sum()
perc = ['{:.1f}%'.format(c/total*100) for c in counts]
ax = counts.plot.pie(colors=col, wedgeprops=dict(width=0.5), labels=perc, figsize=(3, 2), legend=False)
ax.set_ylabel('')
ax.figure.savefig('../data/figures/paper/main/5y.pdf', bbox_inches='tight')
fig, axes = plt.subplots(1, len(overlap_counts), figsize=(6, 5))
i = 0
for model, counts in overlap_counts.iterrows():
    col = [bin_colors[c] for c in overlap_counts.columns]
    total = counts.sum()
    perc = ['{:.0f}%'.format(c/total*100) for c in counts]
    ax = counts.plot.pie(colors=col, wedgeprops=dict(width=0.5), ax=axes[i], labels=perc)
    lgnd = ax.legend(overlap_counts.columns, loc='upper center', bbox_to_anchor=(0.5, -0.1))
    if i != 1:
        lgnd.remove()
    ax.set_ylabel('')
    ax.set_xlabel(model)
    i += 1
#plt.tight_layout()
fig.savefig('../data/figures/paper/main/5x.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15836627140974968, 'fraction_numerical': 0.0240918501788067, 'mean_word_length': 5.04184671366841, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 42, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '295479', 'license': 'MIT', 'n_tokens_mistral': 10402, 'n_tokens_neox': 9712, 'n_words': 1479, 'path': 'notebooks/NovelCandidates.ipynb', 'repo_name': 'prihoda/bgc-pipeline-1', 'size': '26565'}