<jupyter_start><jupyter_text>genome size<jupyter_code>for genome_name, faifn, gtffn, chroms in GENOMES:
    t = pd.read_csv(faifn, delimiter='\t', header=None,
                    names=['name','len','offst','linebases','linewidth'])
    total_genome_len = t.loc[t['name'].isin(chroms),'len'].sum()
    print(genome_name, faifn)
    print(chroms)
#     display(t)
    print('total',total_genome_len)
# TRANSFN = "../datafiles/transcript_list_{}".format(GENOME)
# tlist = pd.read_csv(TRANSFN, header=None, names=['transcript_id'])
# t = tlist['transcript_id'].str.rsplit('-', n=1, expand=True)
# t.columns = ['gene', 'splice_id']
# tlist = pd.concat((tlist,t), axis=1)
for genome_name, faifn, gtffn, chroms in GENOMES:
    print(genome_name, gtffn)
    d = pd.read_csv(gtffn, sep='\t', comment='#', header=None,
                names=['seqid',
                       'source',
                       'type',
                       'start',
                       'end',
                       'score',
                       'strand',
                       'phase',
                       'attributes'],
                    dtype={'seqid':str})

    # total number of genes
    t = d.loc[d['type']=='gene' ,:].copy(deep=True)
    t['gene_id'] = t['attributes'].apply(
                        lambda x: dict([_.strip().split() for _ in
                        x.split(';') if _])['gene_id'].strip('"'))
    print("total number of genes in gtf:", t['gene_id'].unique().shape[0])

    # list of CDS
    d = d.loc[d['type']=='CDS' ,:]
    d['gene_id'] = d['attributes'].apply(
                        lambda x: dict([_.strip().split() for _ in
                        x.split(';') if _])['gene_id'].strip('"'))
    d['transcript_id'] = d['attributes'].apply(
                        lambda x: dict([_.strip().split() for _ in
                        x.split(';') if _])['transcript_id'].strip('"'))
    # filter to only CRHOMS
    d = d.loc[d['seqid'].isin(chroms) ,:]
    # cdslist = d

    print('Total coding size (sum of CDS lens)',(d['end']-d['start']+1).sum())

    # assert len(set(cdslist['transcript_id'])-set(tlist['transcript_id'])) == 0 # transcripts should be superset of CDS
    # tcdslist = set(tlist['transcript_id']) & set(cdslist['transcript_id'])
    # tcdsdf = tlist[tlist['transcript_id'].isin(tcdslist)]
    # print('# trascripts with CDS:', len(tcdslist))

    # num_coding_genes = tcdsdf['gene'].unique().shape[0]
    # print('# genes with a transcript with a CDS:', num_coding_genes)<jupyter_output>AgamP4.11 ../datafiles/Anopheles-gambiae-PEST_BASEFEATURES_AgamP4.11.gtf
total number of genes in gtf: 13822
Total coding size (sum of CDS lens) 25728530
Aaegypti_L5.1 ../datafiles/Aedes-aegypti-LVP_AGWG_BASEFEATURES_AaegL5.1.gtf
total number of genes in gtf: 19712
Total coding size (sum of CDS lens) 59121227<jupyter_text>Results summaries from json files<jupyter_code>SETNAMES = [
    'VGL-gam',
    'Ag1000g-gam',
    'VGL-col',
    'Ag1000g-col',
    'VGL-Aaeg',
]
dat_any = OrderedDict()
dat_onecall = OrderedDict()
for setname in SETNAMES:
    try:
        with open(setname+'_results.json','r') as fh:
            dat_any[setname] = OrderedDict(json.load(fh))
    except FileNotFoundError as e:
        print(e)
        pass
    with open(setname+'_ignore-one-call_results.json','r') as fh:
        dat_onecall[setname] = OrderedDict(json.load(fh))
#         display(dat[setname])
s = ''
for label,dat in [['any variation',dat_any], 
                  ['ignore single call variants',dat_onecall],
                 ]:
    s += '<h3>'+label+'</h3>\n'
    s += '<table border="1">\n'

    s += '<tr>\n'
    s += '<th></th>\n'
    for k,d in dat.items():
        s += '<th>{:s}</th>\n'.format(k)
    s += '</tr>\n'

    s += '<tr>\n'
    s += '<th>'+'N samples'+'</th>\n'
    for k,d in dat.items():
        s += '<td>{:d}</td>\n'.format(dat[k]['number of samples'])
    s += '</tr>\n'

    s += '<tr>\n'
    s += '<th>'+'good unique targets'+'</th>\n'
    for k,d in dat.items():
        s += '<td>{:d}<br>{:.2f}% T<br>{:.2f}% P</td>\n'.format(
            dat[k]['good unique targets'],
            dat[k]['good unique targets % of total unique'],
            dat[k]['good unique targets % of potential'])
    s += '</tr>\n'

    s += '<tr>\n'
    s += '<th>'+'coding transcripts with good targets'+'</th>\n'
    for k,d in dat.items():
        s += '<td>{:d}<br>{:.2f}%</td>\n'.format(
            dat[k]['coding transcripts w/ good targets'],
            dat[k]['coding transcripts w/ good targets %'])
    s += '</tr>\n'

    s += '<tr>\n'
    s += '<th>'+'coding genes with good targets'+'</th>\n'
    for k,d in dat.items():
        s += '<td>{:d}<br>{:.2f}%</td>\n'.format(
            dat[k]['coding genes w/ good targets'],
            dat[k]['coding genes w/ good targets %'])
    s += '</tr>\n'

    s += '</table>'

display(HTML(s))
# print(s)<jupyter_output><empty_output><jupyter_text>filtering by probability of target being reference sequence<jupyter_code>def good_counts(tout, num_total_unique_targets, num_coding_genes, min_refAF=None, min_p_ref=None):
    """return counts of good targets and genes derived from tout
    min_refAF filter targets based on the minimum reference allele frequency for for all variants in target
    thereby disregarding low frequency variants
    ...or...
    min_p_ref filter targets based on probability of reference sequence for whole target
    p_ref = the product of the reference freq for all variants in a target 
         ~= probability of reference sequence for whole target
    min_p_ref of 1 means any variants in a target will invalidate it
    """
    
    results_dict = OrderedDict()
    results_dict['total unique targets'] = num_total_unique_targets
    results_dict['coding genes'] = num_coding_genes

    t = tout[['gene','transcript_id','chrom','pos','nv','refAF','p_ref','Genomic location']]
    if min_p_ref is not None:
        good_targets = t[t['p_ref'] >= min_p_ref]
    if min_refAF is not None:
        good_targets = t[t['refAF'].apply(lambda x: min(x,default=1)) >= min_refAF]

    # by target (not put in results_dict)
    print('num potential targets (including non-unique):', t.shape[0])
    print('num good targets (including non-unique):', good_targets.shape[0])
    print("% good targets of potential (including non-unique): {:d}/{:d} = {:0.3f}%".format(
        good_targets.shape[0], t.shape[0], 100*good_targets.shape[0]/t.shape[0]))

    # by unique target
    num_potential_unique_targets = t['Genomic location'].unique().shape[0]
    results_dict['potential unique targets'] = num_potential_unique_targets

    num_good_unique_targets = good_targets['Genomic location'].unique().shape[0]
    results_dict['good unique targets'] = num_good_unique_targets
    results_dict['good unique targets % of total unique'] = (100*
        results_dict['good unique targets']/results_dict['total unique targets'])
    results_dict['good unique targets % of potential'] = (100*
        results_dict['good unique targets']/results_dict['potential unique targets'])

    # by gene
    num_coding_genes = results_dict['coding genes']
    num_potential_genes = t['gene'].unique().shape[0]
    num_good_genes = good_targets['gene'].unique().shape[0]
    results_dict['coding genes w/ potential targets'] = num_potential_genes
    results_dict['coding genes w/ potential targets %'] = (100*
        num_potential_genes/num_coding_genes)
    results_dict['coding genes w/ good targets'] = num_good_genes
    results_dict['coding genes w/ good targets %'] = (100*
        num_good_genes/num_coding_genes)

    return results_dict
SETNAMES = [
    'VGL-gam',
    'Ag1000g-gam',
    'VGL-col',
    'Ag1000g-col',
    'VGL-Aaeg',
]

minRefAF_thresholds = [1, 0.995, 0.990, 0.985, 0.980, 0.5, 0]

saved_results_any = OrderedDict()
res_by_minRefAF = OrderedDict()

for setname in SETNAMES:#[1:2]:
    print(setname)

    # get the main target to variation results... takes a while to load
    tout = pd.read_msgpack(setname+'_tout.msgpack')
    
    # also need a few values from the 'any variation' results
    with open(setname+'_results.json','r') as fh:
        saved_results_any[setname] = OrderedDict(json.load(fh))
    
    for minRefAF_threshold in minRefAF_thresholds:
        if setname not in res_by_minRefAF:
            res_by_minRefAF[setname] = OrderedDict()
        res_by_minRefAF[setname][minRefAF_threshold] = good_counts(tout, 
                                       saved_results_any[setname]['unique location target sites'],
                                       saved_results_any[setname]['coding genes'],
                                       min_refAF=minRefAF_threshold)
for k,v in res_by_minRefAF.items():
    display(k)
    for k2,v2 in v.items():
        display(k2)
        display(v2)
# load the ignore-one-call results
res_onecall = OrderedDict()
for setname in SETNAMES:
    with open(setname+'_ignore-one-call_results.json','r') as fh:
        res_onecall[setname] = OrderedDict(json.load(fh))
res_onecall
fig,axs = mpl.pyplot.subplots(2,1, squeeze=False, figsize=(6,6),
                             gridspec_kw={'height_ratios': [3, 1]})

barwidth = 0.25

ax = axs[0,0]
y = [res_by_minRefAF[s][0.99]['coding genes w/ good targets %'] for s in SETNAMES]
x = np.arange(0, len(y))
ax.bar(x+0*barwidth, y, width=barwidth*.9, label='allowing <1% variants')

y = [res_onecall[s]['coding genes w/ good targets %'] for s in SETNAMES]
ax.bar(x+1*barwidth, y, width=barwidth*.9, label='allowing 1-call variants')

y = [res_by_minRefAF[s][1]['coding genes w/ good targets %'] for s in SETNAMES]
ax.bar(x+2*barwidth, y, width=barwidth*.9, label='no detected variants')

ax.set_xticks(np.arange(0, len(y))+barwidth)
ax.set_xticklabels(SETNAMES)
ax.set_ylabel('% of genes with good targets')
ax.set_ylim((0,100))

ax_leg = axs[1,0]
lh,ll = ax.get_legend_handles_labels()
ax_leg.legend(lh,ll, loc='center')
ax_leg.set_axis_off()

SETNAMES, LABELS = zip(*[
    ['VGL-gam', 'VGL\n'+r'$\mathit{An. gambiae}$'],
    ['Ag1000g-gam', 'Ag1000g\n'+r'$\mathit{An. gambiae}$'],
    ['VGL-col', 'VGL\n'+r'$\mathit{An. coluzzii}$'],
    ['Ag1000g-col', 'Ag1000g\n'+r'$\mathit{An. coluzzii}$'],
    ['VGL-Aaeg', 'VGL\n'+r'$\mathit{Ae. aegypti}$'],
])
fig,axs = mpl.pyplot.subplots(2,1, squeeze=False, figsize=(6,6),
                             gridspec_kw={'height_ratios': [100, 1]})

barwidth = 1

ax = axs[0,0]
y1 = np.array([res_by_minRefAF[s][0.99]['coding genes w/ good targets %'] for s in SETNAMES])
y2 = np.array([res_onecall[s]['coding genes w/ good targets %'] for s in SETNAMES])
y3 = np.array([res_by_minRefAF[s][1]['coding genes w/ good targets %'] for s in SETNAMES])
x = np.arange(0, len(y1))

ax.bar(x, y1-y2, bottom=y2, width=barwidth*.8, label='allowing <1% variants', edgecolor='k', linewidth=0.5)
ax.bar(x, y2-y3, bottom=y3, width=barwidth*.8, label='allowing 1-call variants', edgecolor='k', linewidth=0.5)
ax.bar(x, y3, bottom=0, width=barwidth*.8, label='no detected variants', edgecolor='k', linewidth=0.5)

ax.set_xticks(np.arange(0, len(y)))#+barwidth)
ax.set_xticklabels([LABELS[i]+ #s.replace('-','\n')+
                    '\n$N='+str(res_onecall[s]['number of samples'])+'$' for i,s in enumerate(SETNAMES)])
# ax.set_xticklabels([s.replace('-','\n')+
#                     '\n$N='+str(res_onecall[s]['number of samples'])+'$' for s in SETNAMES])
ax.set_ylabel('% of genes with good targets')
ax.set_ylim((0,100))

ax_leg = axs[1,0]
lh,ll = ax.get_legend_handles_labels()
ax_leg.legend(lh,ll, loc='center')
ax_leg.set_axis_off()
fig.tight_layout()

fig.savefig('percent_good_genes_overview.svg')
fig.savefig('percent_good_genes_overview.png')
SETNAMES, LABELS = zip(*[
    ['VGL-gam', 'VGL\n'+r'$\mathit{An. gambiae}$'],
    ['Ag1000g-gam', 'Ag1000g\n'+r'$\mathit{An. gambiae}$'],
    ['VGL-col', 'VGL\n'+r'$\mathit{An. coluzzii}$'],
    ['Ag1000g-col', 'Ag1000g\n'+r'$\mathit{An. coluzzii}$'],
    ['VGL-Aaeg', 'VGL\n'+r'$\mathit{Ae. aegypti}$'],
])
SETNAMES, LABELS = zip(*[
    ['VGL-gam', 'VGL\n'+r'$\mathit{gambiae}$'],
    ['Ag1000g-gam', 'Ag1000g\n'+r'$\mathit{gambiae}$'],
    ['VGL-col', 'VGL\n'+r'$\mathit{coluzzii}$'],
    ['Ag1000g-col', 'Ag1000g\n'+r'$\mathit{coluzzii}$'],
    ['VGL-Aaeg', 'VGL\n'+r'$\mathit{aegypti}$'],
])

barwidth = .425
b2off = .3

fig,axs = mpl.pyplot.subplots(2,1, squeeze=False, figsize=(4,4),
                             gridspec_kw={'height_ratios': [100, 1]})

ax = axs[0,0]
y1 = np.array([res_by_minRefAF[s][0.99]['coding genes w/ good targets %'] for s in SETNAMES])
y2 = np.array([res_by_minRefAF[s][1]['coding genes w/ good targets %'] for s in SETNAMES])
x = np.arange(0, len(y1))

# ax.bar(x, y1-y2, bottom=y2, width=barwidth*.8, label='allowing <1% variants', edgecolor='k', linewidth=0.5)
ax.bar(x, y1, bottom=0, width=barwidth, label='allowing DRAs with <1% frequency', edgecolor='k', linewidth=0.5)
ax.bar(x+b2off, y2, bottom=0, width=barwidth, label='no detected DRAs', edgecolor='k', linewidth=0.5)

ax.set_xticks(np.arange(0, len(y))+b2off/2)
ax.set_xticklabels([LABELS[i]+ #s.replace('-','\n')+
                    '\n$N='+str(res_onecall[s]['number of samples'])+'$' for i,s in enumerate(SETNAMES)])
ax.tick_params(axis='both', labelsize=8)
# ax.set_xticklabels([s.replace('-','\n')+
#                     '\n$N='+str(res_onecall[s]['number of samples'])+'$' for s in SETNAMES])
ax.set_ylabel('% of genes with good targets')
ax.set_ylim((0,100))

ax_leg = axs[1,0]
lh,ll = ax.get_legend_handles_labels()
ax_leg.legend(lh,ll, loc='center')
ax_leg.set_axis_off()
fig.tight_layout()

fig.savefig('percent_good_genes_overview.svg')
fig.savefig('percent_good_genes_overview.png')
fig,axs = mpl.pyplot.subplots(2,1, squeeze=False, figsize=(6,6),
                             gridspec_kw={'height_ratios': [100, 1]})

barwidth = 1

ax = axs[0,0]
x = np.arange(0, len(y))

y1 = np.array([res_by_minRefAF[s][0.99]['good unique targets % of potential'] for s in SETNAMES])
y2 = np.array([res_onecall[s]['good unique targets % of potential'] for s in SETNAMES])
y3 = np.array([res_by_minRefAF[s][1]['good unique targets % of potential'] for s in SETNAMES])
# y1 = np.array([res_by_minRefAF[s][0.99]['coding genes w/ good targets %'] for s in SETNAMES])
# y2 = np.array([res_onecall[s]['coding genes w/ good targets %'] for s in SETNAMES])
# y3 = np.array([res_by_minRefAF[s][1]['coding genes w/ good targets %'] for s in SETNAMES])

ax.bar(x, y1-y2, bottom=y2, width=barwidth*.8, label='allowing <1% variants', edgecolor='k', linewidth=0.5)
ax.bar(x, y2-y3, bottom=y3, width=barwidth*.8, label='allowing 1-call variants', edgecolor='k', linewidth=0.5)
ax.bar(x, y3, bottom=0, width=barwidth*.8, label='no detected variants', edgecolor='k', linewidth=0.5)

ax.set_xticks(np.arange(0, len(x)))#+barwidth)
ax.set_xticklabels([s.replace('-','\n')+'\n$N='+str(res_onecall[s]['number of samples'])+'$' for s in SETNAMES])
ax.set_ylabel('% of good potential targets')
ax.set_ylim((0,100))

ax_leg = axs[1,0]
lh,ll = ax.get_legend_handles_labels()
ax_leg.legend(lh,ll, loc='center')
ax_leg.set_axis_off()
fig.tight_layout()

# fig.savefig('percent_good_genes_overview.svg')
n_samples = [res_onecall[s]['number of samples'] for s in SETNAMES]
gg = np.array([res_by_minRefAF[s][0.99]['coding genes w/ good targets %'] for s in SETNAMES])
# y2 = np.array([res_onecall[s]['coding genes w/ good targets %'] for s in SETNAMES])
# y3 = np.array([res_by_minRefAF[s][1]['coding genes w/ good targets %'] for s in SETNAMES])

model = sm.OLS(gg, sm.add_constant(n_samples)).fit()
display(model.summary())

fig, ax = mpl.pyplot.subplots(1,1)
ax.plot(n_samples, gg, '.')<jupyter_output>/home/travc/miniconda3/envs/allel/lib/python3.6/site-packages/statsmodels/stats/stattools.py:71: ValueWarning: omni_normtest is not valid with less than 8 observations; 5 samples were given.
  "samples were given." % int(n), ValueWarning)<jupyter_text>targets per gene stats<jupyter_code># setname = 'VGL-gam'
setname = 'Ag1000g-gam'
tout = pd.read_msgpack(setname+'_tout.msgpack')
display(pd.Series(tout['gene'].value_counts()).describe())
setname = 'VGL-Aaeg'
tout = pd.read_msgpack(setname+'_tout.msgpack')
display(pd.Series(tout['gene'].value_counts()).describe())
1-(1-0.03)**47<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18129806213472777, 'fraction_numerical': 0.022946785604429407, 'mean_word_length': 3.690132717830352, 'pattern_counts': {'":': 0, '<': 58, '<?xml version=': 0, '>': 55, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '503997', 'license': 'MIT', 'n_tokens_mistral': 6408, 'n_tokens_neox': 6024, 'n_words': 1224, 'path': 'variants/results_summary.ipynb', 'repo_name': 'travc/Cas9_target_site_survey', 'size': '16255'}