<jupyter_start><jupyter_text>SARS-CoV-2 section of manuscript<jupyter_code>import os
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import numpy as np
import scipy.stats as stats
import statsmodels
import glob

# Suppress FutureWarnings
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)

# SMEAGOL imports
import smeagol.io, smeagol.models, smeagol.scan, smeagol.enrich, smeagol.visualize, smeagol.variant<jupyter_output><empty_output><jupyter_text>Read genome, PWMs, sites, and enrichment<jupyter_code>genome_file='../DATA/Genomes/Genomes/IV/Monopartite/Coronaviridae/GCF_009858895.2/GCF_009858895.2_ASM985889v3_genomic.fna.gz'
genome = smeagol.io.read_fasta(genome_file)
pwm_file='../DATA/PWMs/attract_rbpdb_encode_filtered_human_pwms.h5'
pwms = pd.read_hdf(pwm_file, key="data")
print(len(pwms))
print(len(set(pwms.Gene_name)))
# Representative PWMs
reps = list(pd.read_csv('../DATA/PWMs/attract_rbpdb_encode_representative_matrices.txt', header=None)[0])
sites = pd.read_csv('../RESULTS/sites_attract_rbpdb_encode_filtered_human_pwms.tsv', sep='\t')
sites = sites.merge(pwms.loc[:, ['Matrix_id', 'Gene_name', 'source']])
all_enr = pd.read_csv('../RESULTS/genome_enrichment_results.tsv', sep='\t')
pre_len = len(all_enr)
all_enr = all_enr.merge(pwms.loc[:, ['Matrix_id', 'Gene_name', 'source']])
assert len(all_enr) == pre_len

all_enr['sig05'] = all_enr.fdr<.05
all_enr['sig1'] = all_enr.fdr<.1
all_enr['dir'] = (all_enr.num > all_enr.avg).map({True: 'up', False: 'dn'})
all_enr['fc'] = all_enr.num/all_enr.avg
all_enr['log2fc'] = np.log2(all_enr.fc)

enr = all_enr[all_enr.assembly=='GCF_009858895.2'].reset_index(drop=True)
cov_enr = all_enr[(all_enr.family == 'Coronaviridae') & \
                  (all_enr.assembly!='GCF_009858895.2')].reset_index(drop=True)
uenr = pd.read_csv('../RESULTS/utr_enrichment_results.tsv', sep='\t')
uenr = uenr[uenr.assembly == 'GCF_009858895.2']
uenr['sig05'] = uenr.fdr<.05
uenr['sig1'] = uenr.fdr<.1
uenr['dir'] = (uenr.num > uenr.avg).map({True: 'up', False: 'dn'})
uenr['fc'] = uenr.num/uenr.avg
uenr['log2fc'] = np.log2(uenr.fc)<jupyter_output>/Users/avantika/Desktop/viruslandscape/vir2/lib/python3.7/site-packages/pandas/core/arraylike.py:358: RuntimeWarning: divide by zero encountered in log2
  result = getattr(ufunc, method)(*inputs, **kwargs)<jupyter_text>Binding sites (threshold 0.8)<jupyter_code>sites.sense.value_counts()
# How many RBPs are predicted to bind?
pos_bound_rbps = set(sites.Gene_name[sites.sense=='+'])
neg_bound_rbps = set(sites.Gene_name[sites.sense=='-'])
all_bound_rbps = set(sites.Gene_name)
print("RBPs predicted to bind to the + strand: {}\n".format(len(pos_bound_rbps)), 
      "RBPs predicted to bind to the - strand: {}\n".format(len(neg_bound_rbps)), 
      "RBPs predicted to bind to either strand: {}\n".format(len(all_bound_rbps)))
# Which RBPs have the most sites on the + strand?
sites.loc[sites.sense=='+', ['Matrix_id', 'Gene_name']].value_counts().head()
# How many sites and RBPs are predicted for the UTRs?
usites = pd.read_csv('../RESULTS/gpIV_utr_sites.tsv', sep='\t')
usites = usites[usites.genome_file == 'GCF_009858895.2_ASM985889v3_genomic.fna.gz']
print("Number of sites on the 5'UTR: {}\n".format(len(usites[usites.region=='five_prime_UTR'])), 
      "RBPs predicted to bind to the 5'UTR: {}\n".format(len(set(usites.Gene_name[usites.region=='five_prime_UTR']))),
      "Number of sites on the 3'UTR: {}\n".format(len(usites[usites.region=='three_prime_UTR'])), 
      "RBPs predicted to bind to the 3'UTR: {}\n".format(len(set(usites.Gene_name[usites.region=='three_prime_UTR']))))<jupyter_output>Number of sites on the 5'UTR: 187
 RBPs predicted to bind to the 5'UTR: 63
 Number of sites on the 3'UTR: 423
 RBPs predicted to bind to the 3'UTR: 70<jupyter_text>Enrichment<jupyter_code>print("Significant motifs: {}\n".format(len(set(enr.Matrix_id[enr.sig05]))),
      "Significant RBPs: {}\n".format(len(set(enr.Gene_name[enr.sig05]))))
# How many enriched and depleted RBPs on each strand?
enr.loc[enr.sig05, ['Gene_name', 'sense', 'dir']].drop_duplicates().iloc[:, 1:].value_counts()
print("Significant motifs in UTRs: {}\n".format(sum(uenr.sig05)),
      "Significant RBPs in UTRs: {}\n".format(len(set(uenr.RBP[uenr.sig05]))))<jupyter_output>Significant motifs in UTRs: 9
 Significant RBPs in UTRs: 9<jupyter_text>Define enriched sets<jupyter_code>pos_up = set(enr.Gene_name[(enr.sense == '+') & (enr.sig05) & (enr.dir == 'up')])
neg_up = set(enr.Gene_name[(enr.sense == '-') & (enr.sig05) & (enr.dir == 'up')])
pos_dn = set(enr.Gene_name[(enr.sense == '+') & (enr.sig05) & (enr.dir == 'dn')])
neg_dn = set(enr.Gene_name[(enr.sense == '-') & (enr.sig05) & (enr.dir == 'dn')])

f_up=set(uenr.RBP[(uenr.region=='five_prime_UTR') & (uenr.dir == 'up') & (uenr.sig05)])
t_up=set(uenr.RBP[(uenr.region=='three_prime_UTR') & (uenr.dir == 'up') & (uenr.sig05)])
f_dn=set(uenr.RBP[(uenr.region=='five_prime_UTR') & (uenr.dir == 'dn') & (uenr.sig05)])
t_dn=set(uenr.RBP[(uenr.region=='three_prime_UTR') & (uenr.dir == 'dn') & (uenr.sig05)])

pos_up1 = set(enr.Gene_name[(enr.sense == '+') & (enr.fdr<0.1) & (enr.dir == 'up')])
neg_up1 = set(enr.Gene_name[(enr.sense == '-') & (enr.fdr<0.1) & (enr.dir == 'up')])
pos_dn1 = set(enr.Gene_name[(enr.sense == '+') & (enr.fdr<0.1) & (enr.dir == 'dn')])
neg_dn1 = set(enr.Gene_name[(enr.sense == '-') & (enr.fdr<0.1) & (enr.dir == 'dn')])

f_up1=set(uenr.RBP[(uenr.region=='five_prime_UTR') & (uenr.dir == 'up') & (uenr.fdr<.1)])
t_up1=set(uenr.RBP[(uenr.region=='three_prime_UTR') & (uenr.dir == 'up') & (uenr.fdr<.1)])
f_dn1=set(uenr.RBP[(uenr.region=='five_prime_UTR') & (uenr.dir == 'dn') & (uenr.fdr<.1)])
t_dn1=set(uenr.RBP[(uenr.region=='three_prime_UTR') & (uenr.dir == 'dn') & (uenr.fdr<.1)])

all_pos_up = pos_up.union(f_up).union(t_up)
all_up = all_pos_up.union(neg_up)
all_dn = pos_dn.union(f_dn).union(t_dn).union(neg_dn)

all_pos_up1 = pos_up1.union(f_up1).union(t_up1)
all_up1 = all_pos_up1.union(neg_up1)
all_dn1 = pos_dn1.union(f_dn1).union(t_dn1).union(neg_dn1)

print("Enriched:")
print("In + : {} ".format(len(pos_up)), 
      "In -: {} ".format(len(neg_up)), 
      "In 5'UTR: {} ".format(len(f_up)), 
      "In 3'UTR: {} ".format(len(t_up)), 
      "In all + : {} ".format(len(all_pos_up)), 
      "In all: {} \n".format(len(all_up)))
print("Depleted:")
print("In + : {} ".format(len(pos_dn)), 
      "In -: {} ".format(len(neg_dn)), 
      "In 5'UTR: {} ".format(len(f_dn)), 
      "In 3'UTR: {}".format(len(t_dn)), 
      "In all: {}\n".format(len(all_dn)))

print("Enriched at FDR .1:")
print("In + : {} ".format(len(pos_up1)), 
      "In -: {} ".format(len(neg_up1)), 
      "In 5'UTR: {} ".format(len(f_up1)), 
      "In 3'UTR: {} ".format(len(t_up1)), 
      "In all + : {} ".format(len(all_pos_up1)), 
      "In all: {}\n".format(len(all_up1)))
print("Depleted at FDR .1:")
print("In + : {} ".format(len(pos_dn1)), 
      "In -: {} ".format(len(neg_dn1)), 
      "In 5'UTR: {} ".format(len(f_dn1)), 
      "In 3'UTR: {}".format(len(t_dn1)), 
      "In all: {}\n".format(len(all_dn1)))
# Differences between + and -
print(set(pos_up).difference(neg_up))
print(set(neg_up).difference(pos_up))
print(set(pos_dn).difference(neg_dn))
print(set(neg_dn).difference(pos_dn))
# Differences between UTR and full
print(f_up.difference(pos_up))
print(t_up.difference(pos_up))
print(f_up.intersection(pos_dn))
print(t_up.intersection(pos_dn))<jupyter_output>{'ZRANB2'}
{'PPIE', 'PABPN1L', 'KHDRBS3'}
set()
set()<jupyter_text>Intersect with known binders<jupyter_code>schmidt_b_core = set(pd.read_csv('../../sars-cov-2/binding_data/schmidt_core.txt', header=None)[0])
schmidt_b_exp = set(pd.read_csv('../../sars-cov-2/binding_data/schmidt_expanded.txt', header=None)[0])
schmidt_b = schmidt_b_core.union(schmidt_b_exp)

kamel_b = pd.read_excel('~/Desktop/viruslandscape/sars-cov-2/binding_data/kamel.xlsx')
kamel_b = set(kamel_b['Gene.name'][(kamel_b['Sig.vRNA.Interactors.FDR10'] == 'Yes') & (~kamel_b['Gene.name'].str.startswith('sp|'))])

lee_b = set(pd.read_csv('../../sars-cov-2/binding_data/lee.txt', header=None)[0])

all_validated = schmidt_b.union(kamel_b).union(lee_b)


print("Validated binders from Schmidt study : {}\n".format(len(schmidt_b)),
      "Validated binders from Schmidt study with PWMs: {}\n".format(len(set(pwms.Gene_name).intersection(schmidt_b))),
      "Validated binders from Kamel study : {}\n".format(len(kamel_b)),
      "Validated binders from Kamel study with PWMs: {}\n".format(len(set(pwms.Gene_name).intersection(kamel_b))),
      "Validated binders from Lee study : {}\n".format(len(lee_b)),
      "Validated binders from Lee study with PWMs: {}\n".format(len(set(pwms.Gene_name).intersection(lee_b))),
      "Validated binders from all studies : {}\n".format(len(all_validated)),
      "Validated binders from all studies with PWMs: {}\n".format(len(set(pwms.Gene_name).intersection(all_validated))),
      "Validated binders from all studies with sites: {}".format(len(set(sites.Gene_name).intersection(all_validated))))

print(len(all_up), len(all_dn))
print("Validated binders enriched at fdr .05 : {}\n".format(len(all_up.intersection(all_validated))),
      "Validated binders depleted at fdr .05 : {}\n".format(len(all_dn.intersection(all_validated))),
      "Validated binders enriched at fdr .1 : {}\n".format(len(all_up1.intersection(all_validated))),
      "Validated binders depleted at fdr .1 : {}\n".format(len(all_dn1.intersection(all_validated))))
all_dn.intersection(all_validated)<jupyter_output><empty_output><jupyter_text>Compare to CRISPR screens<jupyter_code>crispr_df = pd.read_csv('../RESULTS/SARS_CoV_2_validated_RBP_effectors.tsv', sep='\t')
crispr_df.head(3)
crispr_df['proviral'] = (crispr_df.Lee_et_al_proviral) | (crispr_df.Flynn_et_al_proviral) \
| (crispr_df.Wei_et_al_proviral)
crispr_df['antiviral'] = (crispr_df.Lee_et_al_antiviral) | (crispr_df.Flynn_et_al_antiviral) \
| (crispr_df.Wei_et_al_antiviral)
all_pv = set(crispr_df.Gene[crispr_df['proviral']])
all_av = set(crispr_df.Gene[crispr_df['antiviral']])
all_crispr = set(crispr_df.Gene)
print(

"Antiviral in any study : {}\n".format(len(all_av)),
"Antiviral in any study with PWMs: {}\n".format(sum((crispr_df.antiviral) & (crispr_df.in_RBP_db))),
"Antiviral in any study with enrichment: {}\n".format(all_up.intersection(all_av)),
"Antiviral in any study with depletion: {}\n".format(all_dn.intersection(all_av)),
"Antiviral in any study with depletion at .1: {}\n".format(all_dn1.intersection(all_av)),

"\n",

"Proviral in any study : {}\n".format(len(all_pv)),
"Proviral in any study with PWMs: {}\n".format(sum((crispr_df.proviral) & (crispr_df.in_RBP_db))),
"Proviral in any study with enrichment: {}\n".format(all_up.intersection(all_pv)),
"Proviral in any study with enrichment at .1: {}\n".format(all_up1.intersection(all_pv)),
"Proviral in any study with depletion: {}\n".format(all_dn.intersection(all_pv)),

"\n",

"Found in any study : {}\n".format(len(all_crispr)),
"Found in any study with PWMs: {}\n".format(len(crispr_df[crispr_df.in_RBP_db]))
)
# Antiviral and depleted
enr.loc[(enr.Gene_name.isin(all_av)) & (enr.sig1) & (enr.dir=='dn'), 
        ['RBP', 'Matrix_id', 'fc']].sort_values('fc')
smeagol.visualize.plot_pwm(pwms, 'M150_0.6', 5)
smeagol.visualize.plot_pwm(pwms, 'encode_FUBP3_2', 5)
smeagol.visualize.plot_pwm(pwms, 's101', 5)
smeagol.visualize.plot_pwm(pwms, 'encode_PCBP2_2', 5)
smeagol.visualize.plot_pwm(pwms, 'encode_HNRNPA2B1_1', 5)
smeagol.visualize.plot_pwm(pwms, 'encode_TARDBP_1', 5)
smeagol.visualize.plot_pwm(pwms, 'M013_0.6', 5)
smeagol.visualize.plot_pwm(pwms, '75', 5)<jupyter_output><empty_output><jupyter_text>Novel findings<jupyter_code>enr['binder'] = enr.RBP.isin(all_validated)
uenr['binder'] = uenr.RBP.isin(all_validated)
enr['in_crispr'] = enr.RBP.isin(all_crispr)
uenr['in_crispr'] = uenr.RBP.isin(all_crispr)
enr['proviral'] = enr.RBP.isin(all_pv)
enr['antiviral'] = enr.RBP.isin(all_av)
uenr['proviral'] = uenr.RBP.isin(all_pv)
uenr['antiviral'] = uenr.RBP.isin(all_av)
# RBPs which are clearly enriched, not found in screens
enr_novel_hits = enr[(enr.sense=='+') & (enr.sig05) & (enr.dir=='up') & (enr.fc >= 2) & \
                    (~enr.in_crispr) & (~enr.RBP.isin(all_dn))]
enr_novel_hits[['RBP', 'Matrix_id', 'num', 'avg', 'fc', 'binder']]
uenr_novel_hits = uenr[(uenr.sig05) & (uenr.dir=='up') & (uenr.fc >= 2) & \
                 (~uenr.in_crispr) & (~uenr.RBP.isin(all_dn1))]
uenr_novel_hits[['RBP', 'Matrix_id', 'region', 'num', 'avg', 'fc', 'binder']]
#Clearly Depleted, not found in screens
dep_novel_hits = enr[(enr.sense=='+') & (enr.sig05) & (enr.dir=='dn') & (enr.fc <= 0.5) &
               (~enr.in_crispr) & (~enr.RBP.isin(all_up1))]
dep_novel_hits[['RBP', 'Matrix_id', 'num', 'avg', 'fc', 'binder']]
set(dep_novel_hits.RBP)<jupyter_output><empty_output><jupyter_text>Compare to other coronaviruses<jupyter_code>cov_pos_up = set(cov_enr.RBP[(cov_enr.sig05) & (cov_enr.sense=='+') & (cov_enr.dir=='up')])
cov_pos_dn = set(cov_enr.RBP[(cov_enr.sig05) & (cov_enr.sense=='+') & (cov_enr.dir=='dn')])
cov_neg_up = set(cov_enr.RBP[(cov_enr.sig05) & (cov_enr.sense=='-') & (cov_enr.dir=='up')])
cov_neg_dn = set(cov_enr.RBP[(cov_enr.sig05) & (cov_enr.sense=='-') & (cov_enr.dir=='dn')])
print(pos_up.difference(cov_pos_up))
print(pos_dn.difference(cov_pos_dn))
print(neg_up.difference(cov_neg_up))
print(neg_dn.difference(cov_neg_dn))
set(['YBX1', 'SFPQ', 'ZCRB1']).intersection(all_validated)
enr.loc[(enr.Gene_name=='YBX1') & (enr.sig05), 
        ['RBP', 'Matrix_id', 'num', 'avg', 'fc', 'binder']]
smeagol.visualize.plot_pwm(pwms, 'M111_0.6', 6)<jupyter_output><empty_output><jupyter_text>Analyze variant effects<jupyter_code># Select PWMs based on the following criteria:
## (Clearly, strongly depleted novel hits OR antiviral with predicted sites on the + genome)
## AND fewer than 10 binding sites on the + sense genome

strong_dep_mats = set(dep_novel_hits.Matrix_id[dep_novel_hits.fc <= .5])
av_with_sites_mats = set(enr.Matrix_id[(enr.antiviral) & (enr.sense=='+') & (~enr.proviral)])
few_sites_mats = set(enr.Matrix_id[(enr.sense=='+') & (enr.num<10) & (enr.num>0)])

select_mats = (strong_dep_mats.union(av_with_sites_mats)).intersection(few_sites_mats)
enr.loc[(enr.Matrix_id.isin(select_mats)) & (enr.sense=='+'),
   ['RBP', 'Matrix_id', 'num', 'avg', 'fc', 'binder', 'antiviral', 'in_crispr']]
select_mats.remove('M024_0.6')
select_mats
len(select_mats)
set(pwms.Gene_name[pwms.Matrix_id.isin(select_mats)]).intersection(all_av)
select_sites = sites[(sites.Matrix_id.isin(select_mats)) & (sites.sense=='+')]
select_sites.reset_index(drop=True, inplace=True)
len(select_sites)
select_pwms = pwms[pwms.Matrix_id.isin(select_mats)]
select_pwms.reset_index(drop=True, inplace=True)
# Load mutations from GESS
all_var = pd.read_csv('../../sars-cov-2/SNV_till20210507_uniq.tsv', sep='\t')
all_var.rename(columns={'alteration': 'alt', 'position':'pos'}, inplace=True)
all_var['name'] = 'NC_045512.2'
len(all_var)
all_var.head(3)
effects = smeagol.variant.variant_effect_on_sites(select_sites, all_var, genome, select_pwms)
effects.head(3)
escape = effects[effects.variant_score<.5].copy().reset_index(drop=True)
print(len(escape))
escape.head(3)
escape.drop(columns=['max_score', 'number of samples carrying the SNV', 'AA', 'id', 'width'], 
            inplace=True)
escape.rename(columns={'Gene_name': 'RBP'}, inplace=True)
escape.head()
escape.to_csv('../RESULTS/sarscov2_escape_variants.tsv', sep='\t', index=False)
escape_snvs = escape[['SNV', 'location', 'annotations']].drop_duplicates()
len(escape_snvs)
escape_snvs.location.value_counts()
escape_snvs.annotations.value_counts()
escape.RBP.value_counts()<jupyter_output><empty_output><jupyter_text>Focus on position 11078<jupyter_code>select_sites = select_sites[(select_sites.start <= 11078) & (select_sites.end >= 11078)].reset_index(drop=True)
select_var = all_var.loc[all_var.pos==11078, ['SNV', 'pos', 'reference', 'alt', 'name']]
select_var = pd.concat([select_var, 
   pd.DataFrame.from_dict({'SNV':'T11078A','pos':11078,'reference':'T','alt':'A', 'name':'NC_045512.2'}, orient='index').T]).reset_index(drop=True)
e = smeagol.variant.variant_effect_on_sites(select_sites, select_var, genome, select_pwms)
e['diff'] = e.variant_score - e.score
e = e[e.Matrix_id.isin(reps)]
e = e[['SNV', 'alt', 'Matrix_id', 'Gene_name', 'score', 'diff']].drop_duplicates().sort_values('alt').reset_index(drop=True)
e = e.groupby(['SNV', 'alt', 'Gene_name']).max('diff').reset_index()
e
colors = {'C': 'b', 'G': 'r', 'A':'y'}
fig, ax = plt.subplots()
width = 0.35
x = np.arange(3) 
rects1 = ax.bar(x - width/2, e.loc[e.alt=='C', 'diff'], width, label='T11078C')
rects2 = ax.bar(x + width/2, e.loc[e.alt=='G', 'diff'], width, label='T11078G')
ax.set_xticks(x)
ax.set_xticklabels(['CPEB4 M149_0.6', 'ELAVL1 M232.08', 'RALY M150_0.6'])
plt.xticks(rotation = 90)
ax.legend(loc=2)
ax.set_ylabel('Change in PWM score')
fig.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>Visualize binding site distribution<jupyter_code>simple_pwms = pd.read_csv('../RESULTS/simple_pwms.tsv', sep='\t')
enr.loc[(enr.sig05) & (enr.sense=='+') & (enr.dir=='dn') & (enr.RBP.isin(all_validated)), 
        ['Matrix_id', 'RBP', 'num']]
enr.loc[(enr.sig05) & (enr.sense=='+') & (enr.dir=='dn') & (enr.RBP.isin(all_validated)), 
        ['Matrix_id', 'RBP']]
urich_dn = enr.loc[(enr.sig05) & (enr.sense=='+') & (enr.dir=='dn') & \
        (enr.RBP.isin(all_validated)) & \
        (enr.Matrix_id.isin(simple_pwms.Matrix_id[simple_pwms.nucleotide=='T'])), 
        ['Matrix_id', 'RBP', 'num', 'avg', 'log2fc', 'fdr']]
urich_dn
enr.loc[(enr.sense=='+') & (enr.log2fc<(-2)) & (enr.fdr<0.05),
   ['Matrix_id', 'RBP', 'num', 'avg', 'log2fc', 'fdr']]
urich_dn[urich_dn.Matrix_id.isin(reps)]
smeagol.visualize.plot_pwm(pwms, 'M232_0.6')
smeagol.visualize.plot_pwm(pwms, 'M149_0.6')
smeagol.visualize.plot_pwm(pwms, 'M150_0.6')
smeagol.visualize.plot_pwm(pwms, 'encode_TIA1_1')
smeagol.visualize.plot_pwm(pwms, 'encode_HNRNPC_1')
sns.set_theme()
plt.figure(figsize=(8,6))
plt.hist(enr_to_plot.log2fc, range=[-3, 4], bins=50, alpha = 0.5, 
         label='all')
plt.hist(enr_to_plot.log2fc[enr_to_plot.sig05], range=[-3, 4], bins=50, 
         alpha = 1.0, label='significant')
plt.xlabel("log2 foldchange", size=14)
plt.ylabel("Count", size=14)
plt.legend(loc='upper left')
plt.savefig("../RESULTS/SARSCoV2_foldchanges_hist.pdf")
plt.show()
mat='encode_TIA1_1'
r = smeagol.enrich.enrich_in_sliding_windows(sites[sites.sense=="+"], genome, mat, 250)

r['window'] = r.index
r['sig'] = r.padj<0.05
r['window_name'] = r.apply(lambda row: row.id + ':' + str(row.start) + '-' + str(row.end), axis=1)
smeagol.visualize.sliding_window_enrichment_plot(sliding_window_df=r, x_var='window', 
                        y_var='count', xticklabels='window_name', title=mat)
mat='M150_0.6'
r = smeagol.enrich.enrich_in_sliding_windows(sites[sites.sense=="+"], genome, mat, 500)

r['window'] = r.index
r['sig'] = r.padj<0.05
r['window_name'] = r.apply(lambda row: row.id + ':' + str(row.start) + '-' + str(row.end), axis=1)
smeagol.visualize.sliding_window_enrichment_plot(sliding_window_df=r, x_var='window', 
                        y_var='count', xticklabels='window_name', 
                        title=mat + ' (RALY) counts per genome window',
                        file_path='../RESULTS/SARSCoV2_RALY_window_plot.pdf')
mat='M149_0.6'
r = smeagol.enrich.enrich_in_sliding_windows(sites[sites.sense=="+"], genome, mat, 250)

r['window'] = r.index
r['sig'] = r.padj<0.05
r['window_name'] = r.apply(lambda row: row.id + ':' + str(row.start) + '-' + str(row.end), axis=1)
smeagol.visualize.sliding_window_enrichment_plot(sliding_window_df=r, x_var='window', 
                        y_var='count', xticklabels='window_name', title=mat)
mat='encode_HNRNPC_1'
r = smeagol.enrich.enrich_in_sliding_windows(sites[sites.sense=="+"], genome, mat, 250)

r['window'] = r.index
r['sig'] = r.padj<0.05
r['window_name'] = r.apply(lambda row: row.id + ':' + str(row.start) + '-' + str(row.end), axis=1)
smeagol.visualize.sliding_window_enrichment_plot(sliding_window_df=r, x_var='window', 
                        y_var='count', xticklabels='window_name', title=mat)
mat='M232_0.6'
r = smeagol.enrich.enrich_in_sliding_windows(sites[sites.sense=="+"], genome, mat, 250)

r['window'] = r.index
r['sig'] = r.padj<0.05
r['window_name'] = r.apply(lambda row: row.id + ':' + str(row.start) + '-' + str(row.end), axis=1)
smeagol.visualize.sliding_window_enrichment_plot(sliding_window_df=r, x_var='window', 
                        y_var='count', xticklabels='window_name', title=mat)
genome[0].seq[11070:11086].__str__().replace('T', 'U')
sites[(sites.Gene_name=='RALY') & (sites.sense=='+')]
sites[(sites.Matrix_id=='M232_0.6') & (sites.sense=='+')]
sites[(sites.Gene_name=='CPEB4') & (sites.sense=='+')]
gff = pd.read_csv('../DATA/Genomes/Genomes/IV/Monopartite/Coronaviridae/GCF_009858895.2/GCF_009858895.2_ASM985889v3_genomic.gff.gz',
                 compression='gzip', sep='\t', skiprows=7, header=None)
gff<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.20121807280109533, 'fraction_numerical': 0.026061092488551058, 'mean_word_length': 5.830699774266366, 'pattern_counts': {'":': 0, '<': 69, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '295535', 'license': 'MIT', 'n_tokens_mistral': 9666, 'n_tokens_neox': 8920, 'n_words': 1254, 'path': 'notebooks/SARS-CoV-2_example.ipynb', 'repo_name': 'gruber-sciencelab/VirusHostInteractionAtlas', 'size': '21181'}