<jupyter_start><jupyter_text>IntroductionIn this Task we will imagine that we are in a lab studying Clostridioides difficile. We are tasked with creating a set of gRNAs which will target an inserted plasmid (Addgene Plasmid 54622) expressing an enhanced GFP (EGFP). We will find a set of protospacers that are specific to the EGFP coding region and excluding those with a potental off-target effect within the Thermofilum genome.<jupyter_code>import sys
sys.path.append('../')
import crseek
from crseek import utils
from crseek import estimators
from crseek import annotators
from Bio import SeqIO
from Bio import Seq
import gzip
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sbn

sbn.set(style = 'white', font_scale=1.5)

import warnings
warnings.simplefilter('ignore', category = Seq.BiopythonWarning)

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load SequencesStandard Biopython tools can be used to load the Clostridioides difficile 630, complete genome (RefSeq: NC_009089.1) and the AddGene plasmid.<jupyter_code>with open('data/Clostridioides_difficile_630.gb') as handle:
    genome = list(SeqIO.parse(handle, 'genbank'))[0]
    
with open('data/addgene-plasmid-54622-sequence-158642.gbk') as handle:
    plasmid = list(SeqIO.parse(handle, 'genbank'))[0]

for rec in plasmid.features:
    if rec.qualifiers.get('product', [''])[0].startswith('enhanced GFP'):
        print('Found')
        egfp_feature = rec
        break
        
egfp_record = egfp_feature.extract(plasmid)<jupyter_output>Found<jupyter_text>Extract possible targetsUsing the EGFP sequence the `crisprtree` can extract all potential protospacers by finding all sequences upstream of a PAM sequence on each strand.<jupyter_code>possible_targets = utils.extract_possible_targets(egfp_record, pams=('NGG',))   
print('%i possible targets' % len(possible_targets))<jupyter_output>117 possible targets<jupyter_text>Exclude Off-targetsIn order to avoid effecting the bacterium with the gRNA it is important to check for potential off-target effects in the C. diff genome. The `CFDEstimator` can be used to assess this likelihood taking into account the position and nucleotide specific penalties for mismatches.<jupyter_code>#Build the predictor
estimator = estimators.CFDEstimator.build_pipeline()<jupyter_output><empty_output><jupyter_text>`cas-offinder` is used to rapidly find all positions in the C. diff genome that are <=5 mismatches for each of the potential gRNAs.<jupyter_code>possible_binding = utils.cas_offinder(possible_targets, 5, locus = [genome])

# Score each hit across the genome
possible_binding['Score'] = estimator.predict_proba(possible_binding.values)<jupyter_output><empty_output><jupyter_text>Higher numbers imply a larger likelihood of binding. Since we are looking for gRNAs with no off-target effects we can summarize a gRNA by its worst-case off-target hit.<jupyter_code>results = possible_binding.groupby('spacer')['Score'].agg('max')
results.sort_values(inplace=True)
results.head()<jupyter_output><empty_output><jupyter_text>Here are 5 protospacers that are identical to the EGFP plasmid gene but have no predicted effect across the C. diff genome.  VisualizationThe `crisprtree` module is fully back-compatable with the BioPython module. As such, building visualizations becomes trivial.<jupyter_code>from Bio.Seq import reverse_complement
# Label the top-5 best protospacers on the plasmid backbone.

for key, val in results.head().to_dict().items():
    annotators.annotate_grna_binding(key, plasmid, 
                                     estimator, 
                                     exhaustive = False,
                                     extra_qualifiers = {'Off Target Score': val})
from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import GenomeDiagram
from Bio import SeqIO

gd_diagram = GenomeDiagram.Diagram("AddGene Plasmid")
gd_track_for_features = gd_diagram.new_track(1, name="Annotated Features")
gd_feature_set = gd_track_for_features.new_set()

for feature in plasmid.features:
    if feature.type == 'CDS':
        # Label genes in Blue
        color = colors.blue
    elif 'spacer' in feature.qualifiers:
        # Label gRNAs in Green
        color = colors.green
    else:
        # Everything else in Light Blue
        color = colors.lightblue
    gd_feature_set.add_feature(feature, color=color,
                               label=True, label_size=12)



gd_diagram.draw(format="circular",  pagesize=(20*cm,20*cm),
                fragments=4, start=0, end=len(plasmid), circle_core=0.7)

gd_diagram.write("plasmid_circular.png", "png")

fig, ax = plt.subplots(1,1, figsize= (10, 10))

img = plt.imread("plasmid_circular.png")
ax.imshow(img)
ax.set_xticks([])
ax.set_yticks([])
sbn.despine(ax=ax, left=True, bottom=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08522843679573858, 'fraction_numerical': 0.01188281089940586, 'mean_word_length': 4.598623853211009, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '466937', 'license': 'MIT', 'n_tokens_mistral': 1575, 'n_tokens_neox': 1464, 'n_words': 454, 'path': 'paper/task_01.ipynb', 'repo_name': 'DamLabResources/crisprtree', 'size': '4881'}