<jupyter_start><jupyter_text>Notebook 3 - Running BiG-SCAPE and Calculating Beta-DiversityThis notebook helps to run BiG-SCAPE for the annotated genomes and how to calculate beta-diversity for highlighting the most promising cyaniobacterial genomes<jupyter_code>import pandas as pd
import glob
import os
import subprocess
import re
# from Bio import SeqIO
import numpy as np
import networkx
from networkx.algorithms.components.connected import connected_components
from collections import defaultdict
import seaborn as sns<jupyter_output><empty_output><jupyter_text>**Downloading MIBiG BGCs** ```bashcd ./ncbi_antismash/nf_output/wget https://dl.secondarymetabolites.org/mibig/mibig_gbk_2.0.tar.gztar xvzf mibig_gbk_2.0.tar.gzrm mibig_gbk_2.0.tar.gz``` **Running BiG-SCAPE** ```bashpython /home/jovyan/jupyterdata/tiago/dorrestein_nuno/bigscape/BiG-SCAPE/bigscape.py \-i /home/jovyan/jupyterdata/tiago/cyanobiome/ncbi_antismash/nf_output/ \-o /home/jovyan/jupyterdata/tiago/cyanobiome/bigscape_outputs/ \-c 12 --include_singletons --pfam_dir /home/jovyan/jupyterdata/tiago/dorrestein_nuno/bigscape/``` **Gathering BiG-SCAPE Results**<jupyter_code>%%bash

cd ./bigscape_outputs/

rm ./bigscape_all_c030.txt
find ./network_files/2020-05-16_22-19-49_hybrids_glocal/ -type f -name *_c0.30.network | cat > ./network_files.txt
while read line; do cat $line >> ./pre_bigscape_all_c030.txt; done < ./network_files.txt
awk '/Clustername/&&c++>0 {next} 1' ./pre_bigscape_all_c030.txt >> ./bigscape_all_c030.txt
def get_bigscape_df(input_file):
    bigscape_df = pd.read_csv(input_file,sep='\t')
    bigscape_df.rename(columns=lambda x: re.sub(' ','_',x), inplace=True)
    bigscape_df = bigscape_df[bigscape_df.Clustername_1 != bigscape_df.Clustername_2]
    bigscape_df.reset_index(inplace=True,drop=True)
    return bigscape_df

bigscape_df = get_bigscape_df("./bigscape_outputs/bigscape_all_c030.txt")

bigscape_df[:5]<jupyter_output><empty_output><jupyter_text>**Replacing Distance by Similarity**<jupyter_code>sim_col = []

for i,r in bigscape_df.iterrows():
    similarity = 1 - float(bigscape_df['Raw_distance'].loc[i])
    sim_col.append(similarity)

bigscape_df['Raw_distance'] = sim_col
bigscape_df.rename(columns={'Raw_distance':'Similarity'}, inplace=True)

bigscape_df[:5]<jupyter_output><empty_output><jupyter_text>**Replacing Names for More Useful Labels**<jupyter_code>def parse_gbk_list(folder_list):
    gbk_list,new_name_list = [],[]
    for antismash_folder in folder_list:
        print("Parsing BGCs from folder %s"%antismash_folder)
        for root, dirs, files in os.walk(antismash_folder):
            count = 1
            for file in sorted(files):
                if file.endswith(".gbk"):
                    if 'region' in file:
                        if 'BGC' in file:
                            new_name = file.rstrip('.gbk')
                            new_name_list.append(new_name)
                            gbk_list.append(file.rstrip('.gbk'))
                        else:
                            strain_name = os.path.basename(os.path.normpath(root)).split('.')[0]
                            gbk_list.append(file.rstrip('.gbk'))
                            new_name = strain_name + '.region' + "{0:0=3d}".format(count)
                            count += 1
                            new_name_list.append(new_name)
    return gbk_list,new_name_list

folder_list = ["/Users/tiagoferreiraleao/Dropbox/tiago-NAS/cyanobiome/ncbi_antismash/nf_output/"]

gbk_list,new_name_list = parse_gbk_list(folder_list)
len(gbk_list),len(new_name_list)
name_dict = dict(zip(gbk_list, new_name_list))
new_col1, new_col2 = [],[]

for i,r in bigscape_df.iterrows():
    new_col1.append(name_dict[bigscape_df['Clustername_1'].loc[i]])
    new_col2.append(name_dict[bigscape_df['Clustername_2'].loc[i]])
    
bigscape_df['Clustername_1'] = new_col1
bigscape_df['Clustername_2'] = new_col2

bigscape_df[:5]<jupyter_output><empty_output><jupyter_text>**Obtaining Gene Cluster Families**<jupyter_code>def get_neighbors(target,dataframe,column1,column2):
    subset1 = dataframe[(dataframe[column1]==target)]
    subcat = subset1.append(dataframe[(dataframe[column2]==target)])
    temp_list = []
    for index,row in subcat.iterrows():
        temp_list.append(subcat[column1][index])
        temp_list.append(subcat[column2][index])
    temp_list = list(np.unique(temp_list))
    return temp_list

def to_edges(l):
    '''This function convert a list of group of nodes into edges
    Script adapted from https://stackoverflow.com/questions/4842613/merge-lists-that-share-common-elements/4842783
    '''
    it = iter(l)
    last = next(it)
    for current in it:
        yield last, current
        last = current

def to_graph(l):
    '''This function generate a graph using a list of group of nodes and the edges create with the function to_edges()
    Script adapted from https://stackoverflow.com/questions/4842613/merge-lists-that-share-common-elements/4842783
    '''
    G = networkx.Graph()
    for part in l:
        G.add_nodes_from(part)
        G.add_edges_from(to_edges(part))
    return G

def get_family_dict(components_list,dataframe,dictionary):
    count = 0
    for family in list(components_list):
        count += 1
        for fam_member in family:
            dictionary['GCF%s'%count].append(fam_member)
    return dictionary

def main_get_families(input_file,bigscape_df):
    targets_list = np.unique([bigscape_df.Clustername_1,bigscape_df.Clustername_2])
    neighbors_list = []
    for target in targets_list:
        neighbors_list.append(get_neighbors(target,bigscape_df,'Clustername_1','Clustername_2'))
    G = to_graph(neighbors_list)
    C = connected_components(G)
    gcf_dict = defaultdict(list)
    gcf_dict = get_family_dict(C,bigscape_df,gcf_dict)
    return gcf_dict
gcf_dict = main_get_families('./bigscape_outputs/bigscape_all_c030.txt',bigscape_df)
gcf_dict<jupyter_output><empty_output><jupyter_text>**Obtaining GCF Distribution in the Analized Genomes**<jupyter_code>cluster_list = []
strains = []

for item in gcf_dict.values():
    cluster_list = cluster_list + item

for item in cluster_list:
    if 'BGC' in item[:3]:
        strains.append('MIBiG')
    else:
        strains.append(item.split('.')[0])

col_names = list(np.unique(strains))

col_names.remove('')

col_names
family_rows = []
indexes = []

for key,value in gcf_dict.items():
    indexes.append(key)
    row = []
    for strain in col_names:
        if strain == "MIBiG":
            find_substring_in_list = [s for s in value if "BGC" in s]
            sublist = filter(lambda x: "BGC" in x, value)
            sublist = list(sublist)
        else:
            find_substring_in_list = [s for s in value if strain.lower() in s.lower()] # .lower() for case insensitive
            sublist = filter(lambda x: strain in x, value)
            sublist = list(sublist)
            for item in find_substring_in_list:
                if 'BGC' in item[:3]:
                    find_substring_in_list.remove(item)
                    sublist.remove(item)
        if find_substring_in_list:
            row.append(sublist)
        else:
            row.append(0)
    family_rows.append(row)
    
gcf_patterns = pd.DataFrame(family_rows,columns=col_names,index=indexes)
gcf_patterns
gcf_patterns.to_csv("./outputs/gcf_patterns-TFL201112.tsv",sep='\t')<jupyter_output><empty_output><jupyter_text>**Converting GCF Distribution into Binary for Plotting Heatmap**<jupyter_code>family_counts = []

for i,r in gcf_patterns.iterrows():
    new_row = []
    for item in r:
        if item == 0:
            new_row.append(0)
        else:
            new_row.append(1)
    family_counts.append(new_row)
    
gcf_binary = pd.DataFrame(family_counts,columns=col_names,index=indexes)

gcf_binary[:5]<jupyter_output><empty_output><jupyter_text>**Calculating Beta-Diversity and Creating Figure 2b**<jupyter_code>from scipy.spatial import distance
from scipy.cluster import hierarchy
import seaborn as sns
from skbio.diversity import beta_diversity
import matplotlib.pyplot as plt
sns.set()
#sns.set_context("talk")
%matplotlib inline

header = list(gcf_binary.columns.values)
# gcf_binary.index.names = ['GCFs']

ids = gcf_binary.columns.astype(str).tolist()
# data = gcf_binary.as_matrix()
data = gcf_binary[col_names].values
data = data.transpose()
data = data.astype(int)

bc_dm = beta_diversity('jaccard', data, ids)
bc_dm_df = bc_dm.to_data_frame()
sns.set_context("paper")

correlations = bc_dm_df.corr()
correlations_array = np.asarray(bc_dm_df.corr())

row_linkage = hierarchy.linkage(distance.pdist(correlations_array), method='average')

col_linkage = hierarchy.linkage(distance.pdist(correlations_array.T), method='average')

palette = sns.cubehelix_palette(10)

sns.clustermap(correlations, row_linkage=row_linkage, col_linkage=col_linkage,
               figsize=(13, 13),standard_scale=1,cmap=palette)

print("Jaccard Similarity (0.0-1.0):")

# plt.savefig("./big-scape/outputs/beta-diversity_%s.png"%(int(cutoff*10)),format="png")
cyanobiome_list = ['JAAHII01','JAAHIH01','JAAHGF01','JAAHHC01','JAAHHD01','JAAHFN01','JAAHGM01','JAAHGB01','JAAHFW01','JAAHFZ01','JAAHFP01','JAAHGH01','JAAHGC01','JAAHGA01','JAAHHJ01','JAAHGO01','JAAHGU01','JAAHGW01','JAAHGS01','JAAHHM01','JAAHHP01','JAAHHN01','JAAHGL01','JAAHFU01']
most_promising_color = {}

for key,value in diversity_scores.items():
    if key in cyanobiome_list:
        most_promising_color[key] = {"diversity":"cyanobiome", "subject":"Cyanobiome sample"}
        print(diversity_scores[key])
    else:
        if value[0] >= 0.95:
            most_promising_color[key] = {"diversity":"big_producer", "subject":"High diversity (avg. ≥ 95%)"}
        else:
            most_promising_color[key] = {"diversity":"small_producer", "subject":"Low diversity (avg. < 95%)"}
            
sample_md = pd.DataFrame.from_dict(most_promising_color, orient='index')

sample_md
from skbio.stats.ordination import pcoa

bc_pc = pcoa(bc_dm)
fig = bc_pc.plot(sample_md, 'subject', axis_labels=('PC 1', 'PC 2', 'PC 3'), cmap='jet', s=300)

fig.set_size_inches(20, 10)
len(sample_md[sample_md['diversity'] == 'big_producer'])<jupyter_output><empty_output><jupyter_text>**Obtaining List of High Diversity Samples for Dataset S1 (sheet 3)**<jupyter_code>metadata_df = pd.read_csv("./cyanobiome_metadata_df-TFL200507.tsv",sep='\t',index_col=0)

taxa_dict = dict(zip(metadata_df['GenomeID'],metadata_df['Taxa']))

taxa_dict
col1,col2,col3 = [],[],[]

high_sample_md = sample_md[sample_md['diversity'] == 'big_producer']

for item in high_sample_md.index:
    if item != 'MIBiG':
        bgc_count = len(glob.glob('./ncbi_antismash/nf_output/%s.fasta_output_results/*.region*.gbk'%item))
        col1.append(item)
        col2.append(taxa_dict[item])
        col3.append(bgc_count)
        frames = {"Genome ID":col1,"Taxa ID":col2,"# of BGCs":col3}
        high_bio_df = pd.DataFrame(data=frames)
    
high_bio_df
high_bio_df.to_csv("./outputs/high_bio_df-TFL201112.tsv",sep='\t')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12809618513831134, 'fraction_numerical': 0.02061110106671488, 'mean_word_length': 4.102859778597786, 'pattern_counts': {'":': 10, '<': 39, '<?xml version=': 0, '>': 44, 'https://': 3, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '465695', 'license': 'MIT', 'n_tokens_mistral': 4167, 'n_tokens_neox': 3983, 'n_words': 649, 'path': 'notebook3-running_bigscape.ipynb', 'repo_name': 'tiagolbiotech/cyanobiome', 'size': '11062'}