<jupyter_start><jupyter_text>This tutorial takes you through the basics of analysing Mitty data with some help from cytoolz and pandas<jupyter_code>%load_ext autoreload
%autoreload 2
import time
import matplotlib.pyplot as plt
import cytoolz.curried as cyt
from bokeh.plotting import figure, show, output_file
import mitty.analysis.bamtoolz as bamtoolz
import mitty.analysis.bamfilters as mab
import mitty.analysis.plots as mapl
# import logging
# FORMAT = "[%(filename)s:%(lineno)s - %(funcName)20s() ] %(message)s"
# logging.basicConfig(format=FORMAT, level=logging.DEBUG)
# fname = '../../../mitty-demo-data/filter-demo-data/HG00119-bwa.bam'
# scar_fname = '../../../mitty-demo-data/generating-reads/HG00119-truncated-reads-corrupt-lq.txt'

# 180255 read pairs
fname = '../../../mitty-demo-data/alignment-accuracy/HG00119-bwa.bam'
scar_fname = '../../../mitty-demo-data/generating-reads/HG00119-reads-corrupt-lq.txt'<jupyter_output><empty_output><jupyter_text>Ex1: Simple processing pipeline---Simple example showing how to create a processing pipeline and run it using cytoolz.pipe The components of the pipeline are1. Parse read qname2. Compute d_err3. Categorize reads into some categories4. Count 'em<jupyter_code>max_d = 200
n1 = mab.parse_read_qnames(scar_fname)  # This gives us a curried function
n2 = mab.compute_derr(max_d=max_d)
f_dict = {
    'd = 0': lambda mate: mate['d_err'] == 0,
    '0 < d <= 50': lambda mate: 0 < mate['d_err'] <= 50,
    '50 < d': lambda mate: 50 < mate['d_err'] <= max_d,
    'WC': lambda mate: mate['d_err'] == max_d + 1,
    'UM': lambda mate: mate['d_err'] == max_d + 2
}
n3 = mab.categorize_reads(f_dict)
n4 = mab.count_reads

pipeline = [n1, n2, n3, n4]
%%time 
counts = cyt.pipe(
    bamtoolz.read_bam_st(bam_fname=fname), *pipeline)
print(counts, sum(counts.values(), 0))<jupyter_output>Counter({'d = 0': 359585, 'nocat': 304, '0 < d <= 50': 276, 'WC': 190, '50 < d': 155}) 360510
CPU times: user 25 s, sys: 454 ms, total: 25.5 s
Wall time: 27.2 s<jupyter_text>Ex2: Single reads and `scatter`---We now take the original pipeline and pass it to scatter to perform the same operation, but in parallel<jupyter_code>%%time 
counts = sum(
    (c for c in bamtoolz.scatter(
        pipeline, bam_fname=fname, ncpus=4)
    ), mab.Counter())
print(counts, sum(counts.values(), 0))
max_d = 200
n1 = mab.parse_read_qnames(scar_fname)  # This gives us a curried function
n2 = mab.compute_derr(max_d=max_d)
f_dict = {
    'd = 0': lambda mate: mate['d_err'] == 0,
    '0 < d <= 50': lambda mate: 0 < mate['d_err'] <= 50,
    '50 < d': lambda mate: 50 < mate['d_err'] <= max_d,
    'WC': lambda mate: mate['d_err'] == max_d + 1,
    'UM': lambda mate: mate['d_err'] == max_d + 2
}
n3 = mab.categorize_reads(f_dict)
n4 = mab.count_reads

pipeline = [n1, n2, n3, n4]
%%time 
counts = sum((c for c in bamtoolz.scatter(pipeline, bam_fname=fname, paired=True, ncpus=2)), mab.Counter())
print(counts, sum(counts.values(), 0))<jupyter_output>Counter({'d = 0': 74249, 'WC': 981, 'UM': 304, '50 < d': 130, 'nocat': 116, '0 < d <= 50': 18}) 75798
CPU times: user 36.9 ms, sys: 18 ms, total: 54.8 ms
Wall time: 1.85 s<jupyter_text>Ex3: Alignment histogram----------------------<jupyter_code>max_d = 200
n1 = mab.parse_read_qnames(scar_fname)  # This gives us a curried function
n2 = mab.compute_derr(max_d=max_d)
p1 = mab.default_histogram_parameters()  # The full 8D histogram
p1
p2 = mab.default_histogram_parameters()  # We'll do a detailed dive into MQ and d_err here
p2['mq_bin_edges'] = list(range(62))
p2['v_bin_edges'] = None
p2['t_bin_edges'] = None
p2['xt_bin_edges'] = None
p2['v_bin_edges'] = None
n3 = mab.histograminator(histogram_def=[p1, p2])
pipeline = [n1, n2, n3]
%%time 
h8 = cyt.pipe(
    bamtoolz.read_bam_paired_st(bam_fname=fname), 
    *pipeline)
h8[0].attrs
h8[1].sum()
p2 = mab.collapse(h8[0], xd1=None, xd2=None)
mapl.plot_hist2D(p2)
plt.show()
p2 = mab.collapse(h8[1], mq1=None)
mapl.plot_hist1D(p2)
plt.show()
p2 = mab.collapse(h8[1], mq2=None)
mapl.plot_hist1D(p2)
plt.show()
p2 = mab.collapse(h8[1], mq1=None, mq2=None)
mapl.plot_hist2D(p2)
plt.show()
max_d = 200
n1 = mab.parse_read_qnames(scar_fname)  # This gives us a curried function
n2 = mab.compute_derr(max_d=max_d)
p1 = mab.initialize_pah(name='FullHist')
n3 = mab.histogramize(pah=p1)

pipeline = [n1, n2, n3]
%%time
h8 = None
for h in bamtoolz.scatter(
    pipeline, 
    bam_fname=fname, paired=True, ncpus=2):
    print(h.sum())
    if h8 is None:
        h8 = h
    else:
        h8 += h   
print(h8.sum())
p2 = mab.collapse(h8, xd1=None, xd2=None)
mapl.plot_hist(p2)
plt.show()
h8.sum()
h10.sum()<jupyter_output><empty_output><jupyter_text>Ex2---More involved example showing, in sequence:1. Reading from a BAM2. Pairing up the reads3. Parse qnames (this is a simulated data set)4. Compte d_err for the reads5. Categorize reads based on d_err6. Count reads in each category7. Filter to keep non-reference reads only. Keep a pair only if both reads are non-ref8. Re-Categorize reads based on d_err9. Re-Count reads in each categoryAt the end the category counts are comparatively plotted.<jupyter_code>r1 = mab.read_bam(bam_fname=fname)
r2 = mab.make_pairs
r3 = mab.parse_read_qnames(scar_fname)
r4 = mab.compute_derr(max_d=200)

f_dict = {
    'd = 0': lambda mate: mate['d_err'] == 0,
    '0 < d <= 50': lambda mate: 0 < mate['d_err'] <= 50,
    '50 < d': lambda mate: 50 < mate['d_err'] < 200,
    'WC': lambda mate: 200 < mate['d_err'],
    'UM': lambda mate: mate['read'].is_unmapped
}
r5 = mab.categorize_reads(f_dict)
all_counts = {}
r6 = mab.count_reads(all_counts) 
r7 = mab.filter_reads(mab.non_ref(), all)
r8 = mab.categorize_reads(f_dict)
nr_counts = {}
r9 = mab.count_reads(nr_counts)

for r in cyt.pipe(r1, r2, r3, r4, r5, r6, r7, r8, r9):
    pass<jupyter_output><empty_output><jupyter_text>The new concept here is the use of a dictionary of filters to supply to the categorization function. The result is stored in the `all_counts` and `nr_counts` dictionaries which need to be preallocated and passed to the counting function which modifes them.<jupyter_code>mapl.plot_read_counts(ax=plt.subplot(1, 1, 1), 
                      counts_l=[all_counts, nr_counts],
                      labels=['All', 'non-ref'],
                      keys=['d = 0', '0 < d <= 50', '50 < d', 'WC', 'UM'], 
                      colors=None)
plt.show()<jupyter_output><empty_output><jupyter_text>Ex3---Alignment metrics plotting example1. Read BAMs2. Compute PairedAlignmentHistogram3. Plot slices of the alignment<jupyter_code>max_d = 200
r1 = mab.read_bam(bam_fname=fname)
r2 = mab.make_pairs
r3 = mab.parse_read_qnames(scar_fname)
r4 = mab.compute_derr(max_d=200)

p1 = mab.initialize_pah(name='FullHist')
mab.histogramize(pah=p1)(cyt.pipe(r1, r2, cyt.take(100000), r3, r4))
#mab.save_pah(p1, 'test_hist.pklz')
p2 = mab.collapse(p1, xd1=None, xd2=None)
mapl.plot_hist(p2)
plt.show()<jupyter_output><empty_output><jupyter_text>The plot suggests that the error in aligning mate1 is largely unrelated to the error in aligning mate2 except for some cases where both mates are off by about the same amount  in the same direction (top, right and bottom left corners) or in opposite directions.<jupyter_code>p2 = mab.collapse(p1, xd1=None, xd2=None)
mapl.plot_hist(p2)
plt.show()
p2 = mab.collapse(p1, mq1=None, mq2=None)
mapl.plot_hist(p2)
plt.show()
p2.attrs
p2 = mab.collapse(p1, xd1=None, mq1=None)
mapl.plot_hist(p2)
plt.show()
p1.coords['v1']
p2_r = mab.collapse(p1, xt=None, v1=(5, 6), v2=(5, 6))
p2_snp = mab.collapse(p1, xt=None, v1=(2, 3), v2=(5, 6))

plt.semilogy(p2_r/p2_r.max(), label='Ref')
plt.semilogy(p2_snp/p2_snp.max(), label='SNP(1)')

plt.legend(loc='best')
plt.xticks(range(p2_r.shape[0]), p2_r.coords['xt'].data, rotation='vertical')
plt.show()
p2_1 = mab.collapse(p1, xt=None, mq1=None)
p2_2 = mab.collapse(p1, xt=None, mq2=None)

fig, ax = plt.subplots(1, 2, figsize=(13, 5))
plt.subplots_adjust(wspace=0.5)

mapl.plot_hist(p2_1.T, ax[0])
mapl.plot_hist(p2_2.T, ax[1])

plt.show()
p2 = mab.collapse(p1, t=None, mq1=None)
mapl.plot_hist(p2)
plt.show()
p2 = mab.collapse(p1, mq1=None)
plt.semilogy(p2)
plt.xticks(range(p2.coords['mq1'].values.size), p2.coords['mq1'].values, rotation='vertical')
plt.show()
p2 = mab.collapse(p1, xd1=None, xd2=None, v1=(1, 2), v2=(5, 6))
mab.plot_hist(p2)
plt.show()
r1 = mab.read_bam(bam_fname=fname)

for r in cyt.pipe(r1, cyt.take(20)):
    tostring(r[0]['read'])
read = r[0]['read']
read.tostring(mab.pysam.AlignmentFile(fname))
mab.pysam.AlignmentFile?

mab.plot_hist(p2)
plt.show()
p2.coords[p2.dims[0]].values
ax = plt.subplot(1,1,1)
mapl.plot_mean_MQ_vs_derr(ax=ax, dmv_mat=dmv_mat, fmt='yo', ms=5)
plt.show()
p1 = mab.initialize_pah()
len(p1.coords)
ax1 = plt.subplot(1,1,1)
mapl.plot_perr_vs_MQ(ax=ax1, dmv_mat=dmv_mat, yscale='log')
ax2 = ax1.twinx()
mapl.plot_read_count_vs_MQ(ax=ax2, dmv_mat=dmv_mat)
ax2.set_ylabel('Read count', color='r')
ax2.tick_params('y', colors='r')
ax1.legend(loc='lower right', fontsize=9)
plt.show()
for n, v_bin_label in enumerate(
    ['Ref', 'SNP', 'DEL <= 10']):
    ax = plt.subplot(1, 3, n + 1)
    mapl.hengli_plot(ax=ax, dmv_mat=dmv_mat, v_bin_label=v_bin_label)
plt.show()<jupyter_output><empty_output><jupyter_text>Ex4: Dataframes--------------<jupyter_code>r1 = mab.read_bam(bam_fname=fname, sidecar_fname=scar_fname)
r2 = mab.compute_derr(max_d=200)
r3 = mab.to_df(tags=['NM'])
df = cyt.pipe(r1, r2, cyt.take(20), r3)
df
r1 = mab.read_bam(bam_fname=fname, sidecar_fname=scar_fname)
r2 = mab.compute_derr(max_d=200)
r3 = mab.make_pairs
r4 = mab.to_df(tags=['NM'])
df = cyt.pipe(r1, r2, r3, cyt.take(20), r4)
df
import io
import pysam
fp = pysam.AlignmentFile(fname)
fout_str = io.StringIO()
pysam.AlignmentFile(fout_str, 'r')
read.tostring(mab.pysam.AlignmentFile(fname)).split('\t')
def fromstring(s, ref_dict):
  """Inverse of pysam.AlignedSegment.tostring(): given a string, create an aligned segment

  :param s:
  :param ref_dict: ref_dict = {r: n for n, r in enumerate(fp.references)}
  :return:
  """
  def _split(_s):
    qname, flag, rname, pos, \
    mapping_quality, cigarstring, \
    rnext, pnext, template_length, seq, qual, *_tg = _s.split('\t')

    flag = int(flag)
    rname = ref_dict[rname]
    pos = int(pos)
    mapping_quality = int(mapping_quality)
    rnext = ref_dict[rnext]
    pnext = int(pnext)
    template_length = int(template_length)

    return qname, flag, rname, pos, \
      mapping_quality, cigarstring, \
      rnext, pnext, template_length, seq, qual, _tg

  # So close, pysam.tostring, so close
  def _tags(_t):
    _tl = _t.split(':')
    if _tl[1] == 'i':
      _tl[2] = int(_tl[2])
    elif _tl[1] == 'f':
      _tl[2] = float(_tl[2])

    return _tl[0], _tl[2], _tl[1]

  r = pysam.AlignedSegment()
  r.qname, r.flag, r.rname, r.pos, \
  r.mapping_quality, r.cigarstring, \
  r.rnext, r.pnext, r.template_length, r.seq, r.qual, tags = _split(s)
  r.set_tags([_tags(t) for t in tags])
  return r
read2 = fromstring(read.tostring(mab.pysam.AlignmentFile(fname)), ref_dict)
read2.tostring(mab.pysam.AlignmentFile(fname)).split()
read.tostring(mab.pysam.AlignmentFile(fname)).split()
fp = bamtoolz.pysam.AlignmentFile(fname)
fp.unmapped
import logging
logging.debug("hello world")
mab.save_histogram(p2, 'test_me.pklz')
p3 = mab.load_histogram('test_me.pklz')
p3.dims
p3.dims
p3<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16051011433597184, 'fraction_numerical': 0.044239226033421286, 'mean_word_length': 4.648782911077993, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '619854', 'license': 'Apache-2.0', 'n_tokens_mistral': 5114, 'n_tokens_neox': 4841, 'n_words': 1048, 'path': 'docs/filter-based-analysis-tutorial/filter-based-analysis-tutorial.ipynb', 'repo_name': 'ozemsbg/Mitty', 'size': '11370'}