<jupyter_start><jupyter_text>Brunello Library Analysis The Brunello sgRNA library was downloaded from the publication by [Doench et al. Nature Biotechnology 2016](https://www.nature.com/articles/nbt.3437), in supplementary table 21. The start and end positions for the BED file from the provided "position of base after cut" because it is known for SpCas9 that the cut occurs 3 bp upstream of the PAM.<jupyter_code>import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

sns.set_context('paper')

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Load file.<jupyter_code>brunello = pd.read_excel('../dat/brunello/STable 21 Brunello.xlsx')
brunello.head()<jupyter_output><empty_output><jupyter_text>Convert to BED, including PAM.<jupyter_code>brunello_bed = pd.DataFrame()
brunello_bed['chrom'] = brunello['Genomic Sequence'].str.split('_').str[1].str.split('.').str[0].str.replace('0','')

starts = []
stops = []
# use this later to filter variants that occur in "N" position
n_positions = []

for ix, row in brunello.iterrows():
    # sense and antisense appear to be mixed up in spreadsheet
    if row['Strand'] == 'antisense':
        starts.append(row['Position of Base After Cut (1-based)'] - 17)
        stops.append(row['Position of Base After Cut (1-based)'] + 5)  # to incorporate PAM
        n_positions.append(row['Position of Base After Cut (1-based)'] + 3)
    else:
        starts.append(row['Position of Base After Cut (1-based)'] - 6)
        stops.append(row['Position of Base After Cut (1-based)'] + 16)
        n_positions.append(row['Position of Base After Cut (1-based)'] - 4)
                     
brunello_bed['start'] = starts
brunello_bed['stop'] = stops
brunello_bed['name'] = brunello.index.astype(str) + '_' + brunello['Target Gene Symbol']
brunello_bed.head()<jupyter_output><empty_output><jupyter_text>Save BED file of sgRNA protospacer and PAM coordinates.<jupyter_code>brunello_bed_chr = brunello_bed.copy()
brunello_bed_chr['chrom'] = 'chr' + brunello_bed['chrom'].astype(str)
brunello_bed.to_csv('../dat/brunello/brunello.bed', sep='\t', index=False, header=False)
brunello_bed_chr.to_csv('../dat/brunello/brunello_chr.bed', sep='\t', index=False, header=False)
brunello_bed = pd.read_csv('../dat/brunello/brunello.bed', sep='\t', header=None, names=['chrom','start','stop','name'])<jupyter_output><empty_output><jupyter_text>BCFtools command to determine which common variants from the 1000 Genomes Project occur in the sgRNAs.`bcftools view -R ../dat/brunello.bed --min-af 0.05 -G ~/path_to_1kg_dat/1kg_all_autosomal_variants.bcf -O v > ../dat/variants_in_guides_init.vcf` Filter these to make sure none of the variants we identified are from the "N" position of the "NGG" PAM. `n_positions` defined above stored the positions where the "N" in the "NGG" PAM is.<jupyter_code>unfilt_vars = pd.read_csv('../dat/brunello/variants_in_guides_init.vcf', sep='\t', header=None,
                         names=['chrom','position','rsid','ref','alt','score','passing','info'], comment='#')<jupyter_output><empty_output><jupyter_text>Determine the number of variants present before filtering out variants in the "N" position:<jupyter_code>len(unfilt_vars)<jupyter_output><empty_output><jupyter_text>Determine the number of variants that occur in the "N" position.<jupyter_code>len(unfilt_vars[unfilt_vars['position'].isin(n_positions)])<jupyter_output><empty_output><jupyter_text>Generate a filtered set of non-N position variants.<jupyter_code>filt_vars = unfilt_vars[~unfilt_vars['position'].isin(n_positions)].copy()<jupyter_output><empty_output><jupyter_text>Determine how many sgRNAs are impacted by the non-N variants.<jupyter_code>brunello_guides = brunello_bed.copy().sort_values(by='chrom')
brun_n_vars = []

for ix, row in brunello_guides.iterrows():
    chrom = row['chrom']
    start = row['start']
    stop = row['stop']
    guide_positions = list(range(start, stop))
    vars_in_guide = filt_vars.query('(chrom == @chrom) and (position >= @start) and (position <= @stop)')
    brun_n_vars.append(len(vars_in_guide))
brunello_guides['n_vars'] = brun_n_vars<jupyter_output><empty_output><jupyter_text>Identify the percentage of guides that have a common variant in the non-N position.<jupyter_code>100.0 * len(brunello_guides.query('n_vars > 0')) / len(brunello_guides)<jupyter_output><empty_output><jupyter_text>Save the results to a file.<jupyter_code>brunello_guides.to_csv('../dat/brunello/brunello_guides_results_1kgp.tsv', sep='\t')<jupyter_output><empty_output><jupyter_text>Visualize the number of variants per sgRNA.<jupyter_code>p = sns.countplot(x=brunello_guides['n_vars'], color='dodgerblue')
p.set_yscale('log')
plt.ylabel('Number of gRNAs')
plt.xlabel('Number of Common Variants (MAF >= 5%)')
plt.title('Common Variants from the\n 1000 Genomes Project\n in the Brunello gRNA library')<jupyter_output><empty_output><jupyter_text>WTC analysis of Brunello sgRNA libraryBCFtools command to figure out which variants in WTC are in Brunello gRNAs.`bcftools view -R ../dat/brunello_chr.bed ~/path_to_wtc_dat/wtc.bcf -O v > ../dat/wtc_variants_in_guides_init.vcf`<jupyter_code>unfilt_vars_wtc = pd.read_csv('../dat/brunello/wtc_variants_in_guides_init.vcf', sep='\t', header=None,
                         names=['chrom','position','random','ref','alt','rsid','score','passing','info','genotype'], comment='#')<jupyter_output><empty_output><jupyter_text>Determine the number of variants present before filtering out variants in the "N" position:<jupyter_code>len(unfilt_vars_wtc)<jupyter_output><empty_output><jupyter_text>Determine the number of variants that occur in the "N" position.<jupyter_code>len(unfilt_vars_wtc[unfilt_vars_wtc['position'].isin(n_positions)])<jupyter_output><empty_output><jupyter_text>Generate a filtered file with only non-N variants.<jupyter_code>filt_vars_wtc = unfilt_vars_wtc[~unfilt_vars_wtc['position'].isin(n_positions)].copy()<jupyter_output><empty_output><jupyter_text>Determine how many sgRNAs are impacted by the non-N variants.<jupyter_code>brunello_guides_wtc = brunello_bed_chr.copy().sort_values(by='chrom')

brun_n_vars_wtc = []

for ix, row in brunello_guides_wtc.iterrows():
    chrom = row['chrom']
    start = row['start']
    stop = row['stop']
    guide_positions = list(range(start, stop))
    vars_in_guide = filt_vars_wtc.query('(chrom == @chrom) and (position >= @start) and (position <= @stop)')
    brun_n_vars_wtc.append(len(vars_in_guide))
brunello_guides_wtc['n_vars'] = brun_n_vars_wtc<jupyter_output><empty_output><jupyter_text>Identify the percentage of guides that have a common variant in the non-N position.<jupyter_code>100.0 * len(brunello_guides_wtc.query('n_vars > 0')) / len(brunello_guides_wtc)<jupyter_output><empty_output><jupyter_text>Save the results to a file.<jupyter_code>brunello_guides_wtc.to_csv('../dat/brunello/brunello_guides_results_wtc.tsv', sep='\t')
p = sns.countplot(x=brunello_guides_wtc['n_vars'], color='dodgerblue')
p.set_yscale('log')
plt.ylabel('Number of gRNAs')
plt.xlabel('Number of Variants')
plt.title('Variants from WTC\n in the Brunello gRNA library')<jupyter_output><empty_output><jupyter_text>Put the plots together for 1000 Genomes and WTC and save to a file.<jupyter_code>brunello_guides['Dataset'] = 'Common variants\n (MAF >= 5%)\n 1000 Genomes'
brunello_guides_wtc['Dataset'] = 'WTC variants'
df_overall = pd.concat([brunello_guides, brunello_guides_wtc])
df_overall.to_csv('../dat/brunello/1kg_wtc_brunello_results.tsv', sep='\t', index=False)
p = sns.countplot(x=df_overall['n_vars'], hue=df_overall['Dataset'])
p.set_yscale('log')
plt.ylabel('Number of gRNAs')
plt.xlabel('Number of Variants')
plt.title('Variants from WTC and \nthe 1000 Genomes Project\n in the Brunello gRNA library')
p = sns.catplot('n_vars', kind='count', row='Dataset', hue='Dataset', 
                data=df_overall, height=3, aspect=2, sharex=False)
(p.set_axis_labels('Number of Variants', 'Number of gRNAs (log scale)')
 .set_titles('{row_name} {row_var}', verticalalignment='top'))
p.fig.get_axes()[0].set_yscale('log')
p.fig.suptitle('Variants in the Brunello gRNA library')
p.savefig('../figs/brunello_wtc_1kgp_countplots.pdf', dpi=300,
                    bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13909365558912387, 'fraction_numerical': 0.009305135951661632, 'mean_word_length': 6.310954063604241, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 89, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, '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': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '475846', 'license': 'MIT', 'n_tokens_mistral': 3165, 'n_tokens_neox': 2978, 'n_words': 598, 'path': 'manuscript_analyses/variants_in_published_libraries/src/Brunello Library Analysis.ipynb', 'repo_name': 'keoughkath/AlleleAnalyzer', 'size': '8275'}