<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Preparing the reads
[Loose et al] published their raw read files on ENA. This script uses four of these sets which contain reads of amplicons. These were processed using different "read until" scripts (or none at all), but that doesn't matter. What does matter is to get as much real reads as possible.
Step1:  Load metadata for 4 datasets
Step2:  The individual filenames will look like this
Step3:  Merging alignment data
[Loose et al] provide all the intermediate data files necessary to recreate their figures. Among these, there are some alignment files in SAM format.
Because it doesn't make sense to classify complement sequences/matches in the Read Until context, we only use the "Template" strands.
Step4:  Unfortunately filenames and sequence names tend to get a bit mangled when going from Fast5 to SAM, for various reasons. As of now, there is no particular convention for naming read files or naming the exported sequences. On the one hand I don't feel like it is a good idea to abuse filenames as character seperated database rows, on the other hand, using the unique read id from the Fast5 File isn't very human-friendly either.
To assign genomic coordinates to the reads, a regular expression extracts four numbers from the file name/query name  making each read unique and matchable.
Step5:  Visualizing the alignments
This is just a simple histogram showing where the "reference_start" values fall.
Step6:  Processing the amplicons
[Loose et al] pooled 11 amplicons. Each read has to be assigned retroactively to one of these, represented by number from 0 to 10.
Step7:  How many reads failed to be assigned?
Step8:  Purge these
Step9:  The number of viable reads is diminishing quickly. But this can't be helped.
How many reads longer than 500 bases are assigned to each amplicon?
Step10:  Unfortunately some amplicons are severely underrepresented, with one going as low as 635 reads.
This is a big problem for dividing the data into training and test sets, because blindly sampling from total pool may skew this balance even further. The algorithms will then bias against the least represented amplicons to gain a bit of extra accuracy, which is not what we want. With ten times as much data we could balance both the training and the test set. As it is, I chose to balance the test set only, to get a more realistic view of the performance. My assumption is that, over multiple repetitions of amplification / library preparation and sequencing runs, the amplicons should be roughly equally distributed.
To balance the test set, 200 reads from each amplicon are chosen. This makes for a very weak test set. But again, this can't be helped at this point.
<ASSISTANT_TASK:>
Python Code:
%load_ext autoreload
%autoreload 2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn
import porekit
import re
import pysam
import random
import feather
%matplotlib inline
Explanation: Preparing the reads
[Loose et al] published their raw read files on ENA. This script uses four of these sets which contain reads of amplicons. These were processed using different "read until" scripts (or none at all), but that doesn't matter. What does matter is to get as much real reads as possible.
End of explanation
directories = ["AmpliconOddEvenControl", "AmpliconOddReadUntil", "AmpliconEvenReadUntil", "Balanced"]
meta_frames = []
for d in directories:
    print(d)
    df = porekit.gather_metadata("/home/andi/nanopore/%s" % d, workers=4)
    print(df.shape[0])
    meta_frames.append(df)
meta = pd.concat (meta_frames)
for df in meta_frames:
    print (df.shape)
Explanation: Load metadata for 4 datasets
End of explanation
meta_frames[0].index.values[0]
meta_frames[1].index.values[0]
meta_frames[2].index.values[0]
meta_frames[3].index.values[0]
Explanation: The individual filenames will look like this:
End of explanation
def sam_to_dataframe(file_name):
    sam = pysam.AlignmentFile(file_name)
    records = []
    for i, segment in enumerate(sam):
        d = dict()
        for k in ["query_name", "reference_start", "reference_end", "mapping_quality", ]:
            d[k] = getattr(segment, k)
        records.append(d)
    alignments = pd.DataFrame.from_records(records)
    return alignments
base  = "/home/andi/nanopore/RUFigs/data"
bams = ["/fig3/RU_dudu/RU_dudu_Template.bam",
        "/fig3/RU_udud/RU_udud_Template.bam",
        "/fig3/NO_RU/NO_RU_Template.bam",
        "/fig4/200/200_Template.bam",
        ]
alignments = pd.concat([sam_to_dataframe(base+file_name) for file_name in bams])
Explanation: Merging alignment data
[Loose et al] provide all the intermediate data files necessary to recreate their figures. Among these, there are some alignment files in SAM format.
Because it doesn't make sense to classify complement sequences/matches in the Read Until context, we only use the "Template" strands.
End of explanation
regexp = re.compile(r'_(?P<a>\d+)_(?P<b>\d+)_ch(?P<c>\d+)_file(?P<d>\d+)')
def extract(s):
    try:
        return "_".join(regexp.search(s).groups())
    except:
        return ""
alignments["alignment_key"] = alignments.query_name.map(extract)
meta["alignment_key"] = meta.index.map(extract)
alignments["alignment_key"].map(lambda s: s.split("_")[0]).unique()
meta["run_number"] = meta["alignment_key"].map(lambda s: s.split("_")[0])
meta2 = meta.reset_index().merge(alignments).set_index("filename")
meta2.shape
meta = meta2
Explanation: Unfortunately filenames and sequence names tend to get a bit mangled when going from Fast5 to SAM, for various reasons. As of now, there is no particular convention for naming read files or naming the exported sequences. On the one hand I don't feel like it is a good idea to abuse filenames as character seperated database rows, on the other hand, using the unique read id from the Fast5 File isn't very human-friendly either.
To assign genomic coordinates to the reads, a regular expression extracts four numbers from the file name/query name  making each read unique and matchable.
End of explanation
f, ax = plt.subplots()
f.set_figwidth(13)
ax.hist(meta.reference_start, bins=110);
Explanation: Visualizing the alignments
This is just a simple histogram showing where the "reference_start" values fall.
End of explanation
amplicons = [(52,1980),
             (2065,3965),
             (4070,5989),
             (6059,7981),
             (8012,9947),
             (10008,11963),
             (12006,13941),
             (14011,15945),
             (16076,17987),
             (18022,19972),
             (20053,21979),
            ]
def amplicon_from_position(pos):
    for i,c in enumerate(amplicons):
        a,b = c
        if a<=pos<=b:
            return i
meta["amplicon"] = meta.reference_start.map(amplicon_from_position)
Explanation: Processing the amplicons
[Loose et al] pooled 11 amplicons. Each read has to be assigned retroactively to one of these, represented by number from 0 to 10.
End of explanation
meta.amplicon.isnull().sum()
Explanation: How many reads failed to be assigned?
End of explanation
meta = meta[np.isnan(meta.amplicon)==False]
meta.shape
Explanation: Purge these:
End of explanation
meta.query("template_length>500").groupby("amplicon").format.count()
Explanation: The number of viable reads is diminishing quickly. But this can't be helped.
How many reads longer than 500 bases are assigned to each amplicon?
End of explanation
sufficient = meta.query("template_length>=500")
all_files = sufficient.index.values
test_files = []
for i in range(11):
    sub = sufficient[sufficient.amplicon==i]
    test_files += list(np.random.choice(sub.index.values, 200))
    
training_files = list(set(sufficient.index.values) - set(test_files))
len(training_files), len(test_files)
test_data = sufficient.ix[np.array(test_files)]
feather.write_dataframe(test_data, "amplicon_test_metadata.feather")
training_data = sufficient.ix[np.array(training_files)]
feather.write_dataframe(training_data, "amplicon_training_metadata.feather")
Explanation: Unfortunately some amplicons are severely underrepresented, with one going as low as 635 reads.
This is a big problem for dividing the data into training and test sets, because blindly sampling from total pool may skew this balance even further. The algorithms will then bias against the least represented amplicons to gain a bit of extra accuracy, which is not what we want. With ten times as much data we could balance both the training and the test set. As it is, I chose to balance the test set only, to get a more realistic view of the performance. My assumption is that, over multiple repetitions of amplification / library preparation and sequencing runs, the amplicons should be roughly equally distributed.
To balance the test set, 200 reads from each amplicon are chosen. This makes for a very weak test set. But again, this can't be helped at this point.
End of explanation
<END_TASK>