<jupyter_start><jupyter_text>Tutorial sgRNA analysis  Python Imports<jupyter_code>1
from antenna import antenna_count_reads
import pandas as pd
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Main Analysis Specify the bam file to examine. This file must have been pre-processed with antenna_tag_reads and must be indexed<jupyter_code>counted_bam_path = ''<jupyter_output><empty_output><jupyter_text>Examine the reads in the file and collect TRS alignement scores and relevant metadata in a pandas array<jupyter_code>sgRNA_scores = antenna_count_reads.load_sgRNA_scores(counted_bam_path)<jupyter_output><empty_output><jupyter_text>Examine the distribution of scores in all possible configurations. Specifically:    - Forward / Reverse: orientation of read wrt reference- R1/R2: Whether this read is read 1 or read 2- 5'/3': If the TRS was identified in the 3' or 5' flanking region of the read- TRS/RC TRS/R TRS/C TRS: Orientation of sequences that was found<jupyter_code>fig1 = antenna_count_reads.plot_read_score_dist(sgRNA_scores,n_bins=40)<jupyter_output><empty_output><jupyter_text>Ideas for automating cutoff selection- Remove 0s- Convert to z-scores- Report fraction of reads that exceed 2 sds- Report mean score of all reads above  Merge information from read pairs<jupyter_code>merged_read_information = antenna_count_reads.merge_read_information(sgRNA_scores)
merged_read_information.plot.scatter(x='p5_o_score_r1',y='p5_o_score_r2')
import numpy as np
merged_read_information.plot.scatter(x='p5_o_score_r1',y='p5_o_score_r2')<jupyter_output><empty_output><jupyter_text>Notes:- 2d histogram matplot hist2d- add alpha + jitter- contour on log scale -- pre-bin this- report proportion of reads in each quadrant- https://numpy.org/doc/stable/reference/generated/numpy.histogram2d.html- https://plotly.com/python/2d-histogram-contour/<jupyter_code>cutoff = 30
# Add histogram (over/under axis)
def annotate_merged_read_information(merge_sg_read_info, cutoff, sgRNA_bam_tag_name="TO"):
    """Count sgRNAs stratifying by orientation"""

    columns_to_count = [
        "p5_o_score_r1",
        "p3_o_score_r1",
        "p5_o_score_r2",
        "p3_o_score_r2",
    ]
    
    TRS_found_column = []

    for read_pair in merge_sg_read_info.itertuples():
        count_read_pair = False

        for cc in columns_to_count:
            if getattr(read_pair, cc) > cutoff:
                count_read_pair = True
                continue  # This should only jump out of inner
        
        TRS_found_column.append(count_read_pair)
        
    
    merge_sg_read_info['TRS_found'] = TRS_found_column
    
    return merge_sg_read_info
        

annot_read_info = annotate_merged_read_information(merged_read_information, 50)
fig, ax = plt.subplot()
plt.hist(annot_read_info[annot_read_info['TRS_found']]['reference_start_r1'],log=True,bins=1000);
plt.hist(annot_read_info[annot_read_info['TRS_found'] == False]['reference_start_r1'],log=True,bins=1000);
merged_read_information
trs_intervals = antenna_count_reads.count_sgRNA(merged_read_information, '../../data/orf_start.bed', cutoff)
intervals_counts = antenna_count_reads.summarize_trs_intervals(trs_intervals)
intervals_counts
# Make jupyter notebook and run 
# nbconvert to generate report<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup