<jupyter_start><jupyter_text>GC-IGR Graphing Worksheet Purpose:This jupyter notebook will aid in choosing a genome, download the relevant files from NCBI, extract intergenic regions, match IGRs to known Rfam annotations, and then use support vector machine classifiers to select IGRs for further analysis.    Necessary Imports and Configuration<jupyter_code>%cd '/home/jovyan/work'

import sys
import os
import pandas as pd
import plotly as py
import plotly.graph_objs as go
import numpy as np
import shutil
import tarfile
import pickle
from sklearn import svm
from sqlalchemy import or_
from ipywidgets import interactive
from Bio import SeqIO, SeqRecord

from src.data.make_dataset import extract_igrs, annotate_igrs, download_genome
from src.visualization.visualize import graph_genome, graph_layout
from src.data.rfam_db import rfam_session, Genome

py.io.orca.config.use_xvfb = True
pd.set_option('display.max_columns', 60)<jupyter_output>/home/jovyan/work<jupyter_text>Step 1: Review Bacterial and Archaeal Genomes with Rfam annotations and select UniprotID (upid)<jupyter_code># Create a connection to local or remote Rfam Database
session = rfam_session()

# Get list of bacterial and archaeal genomes and save them in Genome DF
genome_query = session.query(Genome).filter(or_(Genome.kingdom=='archaea', Genome.kingdom=='bacteria'))

# If necessary filter for completely assembled genomes.
genome_query = genome_query.filter(Genome.assembly_level == 'complete-genome')
genome_list = genome_query.all()
session.close()
genome_df = pd.read_sql_query(genome_query.statement, genome_query.session.bind)

# Display the the genomes numbered 0-9 from the above criteria. 
genome_df.iloc[0:10]<jupyter_output><empty_output><jupyter_text>Step 2: Enter the upid of the genome of interest and graph genome.Look up the Uniprot ID (upid) for a genome of interest here: http://rfam.xfam.org/search?q=entry_type:%22Genome%22 and enter it below for the upid:<jupyter_code>upid = 'UP000001174'

session = rfam_session()
genome =  session.query(Genome).get(upid)
session.close()
download_genome(genome)
igr_df = extract_igrs(genome, igr_length_cutoff=1)
annotated_df = annotate_igrs(genome, igr_df)
scatter_plots = graph_genome(annotated_df)
layout = graph_layout(genome)
fig = go.FigureWidget(data=scatter_plots, layout=layout)
fig<jupyter_output><empty_output><jupyter_text>Step 3: SVM Selection of Genomic Regions Enriched for ncRNAsIf necessary, modify the SVM selection hyperparameters (class_weight_mod, gamma_exp, and c_exp) using the sliders below.<jupyter_code># Define function for interactive modification of selection parameters. 

def prepare_selection(annotated_df):
    y = (annotated_df['category'] != 'No Known RNA') & (annotated_df['category'] != 'sRNA') 
    total_igrs = len(y)
    total_knowns = y.sum()
    total_unknowns = total_igrs - total_knowns 
    
    return (y, total_igrs, total_knowns, total_unknowns)

def build_interactive_fn (annotated_df):
    
    y, total_igrs, total_knowns, total_unknowns = prepare_selection(annotated_df)    

    def interactive_fn(class_weight_mod=0.5,  c_exp=2, gamma_exp=-2,):
        class_weight = {False: total_knowns / total_igrs, True: (total_unknowns / total_igrs * class_weight_mod)}
        svm_clf = svm.SVC(C=10**c_exp, class_weight=class_weight, gamma=10**(gamma_exp), random_state=0)
        svm_clf.fit(annotated_df.loc[:, ["gc", "log_length"]], y)
        selection = pd.Series(svm_clf.predict(annotated_df.loc[:, ["gc", "log_length"]]))
        scatter_plots = graph_genome(annotated_df, selection=selection)
        fig = go.FigureWidget(data=scatter_plots, layout=layout)
        display(fig)
    
    return interactive_fn


interactive_fn = build_interactive_fn(annotated_df)
interactive_plot = interactive(interactive_fn, class_weight_mod=(0.05, 2.0, 0.05), gamma_exp=(-5, 5, 0.25), c_exp=(-5,5,0.25))
interactive_plot<jupyter_output><empty_output><jupyter_text>Step 4: Finalize Selection and Build Blast Script/Data TarfileAfter finalizing the selection in the interactive graph above, execute the following two blocks of code to extract the selected intergenic regions, and prepare a tarfile with the collection of data and scripts necessary for blast analysis.<jupyter_code># Extract the values from the interactive plot
class_weight_mod = interactive_plot.kwargs["class_weight_mod"]
c_exp = interactive_plot.kwargs["c_exp"]
gamma_exp = interactive_plot.kwargs["gamma_exp"]

output_folder="data/interim/{}/selection_{}_{}_{}".format(genome.assembly_acc, class_weight_mod, c_exp, gamma_exp)
if not os.path.exists(output_folder + '/data_files'):
    os.makedirs(output_folder + '/data_files')

# Re-create the selection
y, total_igrs, total_knowns, total_unknowns = prepare_selection(annotated_df)    
class_weight = {False: total_knowns / total_igrs, True: (total_unknowns / total_igrs * class_weight_mod)}
svm_clf = svm.SVC(C=10**c_exp, class_weight=class_weight, gamma=10**gamma_exp, probability=True, random_state=0)
svm_clf.fit(annotated_df.loc[:, ["gc", "log_length"]], y)

# Save the selection classifier to a pickle                                                          
svm_pickle = pickle.dumps(svm_clf)
with open("{}/data_files/svmclf.pickle".format(output_folder,genome.assembly_acc, class_weight_mod, c_exp, gamma_exp), 'wb') as svm_pickle_file:
    svm_pickle_file.write(svm_pickle)
selection = pd.Series(svm_clf.predict(annotated_df.loc[:, ["gc", "log_length"]]))
                                                          
# Save a graph of the genome.
scatter_plots = graph_genome(annotated_df, selection=selection)
layout = graph_layout(genome)
fig = go.FigureWidget(data=scatter_plots, layout=layout)
fig.write_image("{}/data_files/{}_plot.svg".format(output_folder,genome.scientific_name.replace(' ','_')))
py.io.write_json(fig, "{}/data_files/{}_json.plotly".format(output_folder,genome.scientific_name.replace(' ','_')))
                                                          
selected_unknowns = selection & (annotated_df['category'] == 'No Known RNA')

# Save a fasta file with all the selected IGRs
selected_igr_list = [SeqRecord.SeqRecord(row.sequence, id=("{}/{}-{}".format(row.accession, row.start +1, row.end))) 
                     for row in annotated_df.loc[selected_unknowns, ["accession","start","end","sequence"]].itertuples()]

if not os.path.exists(output_folder + '/igr_fastas'):
    os.makedirs(output_folder + '/igr_fastas')

for igr in selected_igr_list:
    outputfilename = "{}/igr_fastas/{}.fasta".format(output_folder, igr.id.replace('/','_'))
    SeqIO.write(igr, outputfilename, "fasta")
    
annotated_df.to_csv("{}/data_files/annotated_igrs.csv".format(output_folder), index=False)
if not os.path.exists(output_folder + '/scripts'):
    os.makedirs(output_folder + '/scripts')

shutil.copy('src/shell/blast_source_template.sh', '{}/scripts/blast_source.sh'.format(output_folder))
shutil.copy('src/shell/blast_run_template.sh', '{}/blast_run.sh'.format(output_folder))

if not os.path.exists(output_folder + '/blast_xml'):
    os.makedirs(output_folder + '/blast_xml')
if not os.path.exists(output_folder + '/output'):
    os.makedirs(output_folder + '/output')

with open("{}/scripts/blast_jobfile.sh".format(output_folder), 'w') as jobfile:
    for igr in selected_igr_list:
        fasta_filename = "igr_fastas/{}.fasta".format(igr.id.replace('/','_'))
        xml_filename =  "blast_xml/{}.xml".format(igr.id.replace('/','_'))
        jobfile.write("source scripts/blast_source.sh; $BLASTCMD {} > {}\n".format(fasta_filename, xml_filename))
        
with tarfile.open("data/export/{}_{}_selection_{}_{}_{}_blastdata.tar.gz".format('_'.join(genome.scientific_name.split(' ')[0:2]), genome.assembly_acc, class_weight_mod, c_exp, gamma_exp), "w:gz") as tar:
    tar.add(output_folder, arcname="{}_{}_selection_{}_{}_{}_blastdata".format('_'.join(genome.scientific_name.split(' ')[0:2]), genome.assembly_acc, class_weight_mod, c_exp, gamma_exp))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12207922029239035, 'fraction_numerical': 0.007747094839435212, 'mean_word_length': 5.2875098193244305, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 6}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '637465', 'license': 'MIT', 'n_tokens_mistral': 2755, 'n_tokens_neox': 2607, 'n_words': 566, 'path': 'notebooks/1-Genome-IGR-Selection.ipynb', 'repo_name': 'snmalk/dimpl', 'size': '8003'}