5. Single-molecule data analysis 2

Contents

5. Single-molecule data analysis 2#

[1]:

from pathlib2 import Path
import papylio as pp
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

%matplotlib inline

Experiment import#

Note that the sequencing data is automatically imported.

[2]:

experiment_path = Path(r'C:\Users\user\Desktop\SPARXS example dataset')

[3]:

data_per_sequence_path = experiment_path / 'Analysis' / 'Datasets per sequence'
data_per_sequence_path

[3]:

WindowsPath('C:/Users/user/Desktop/SPARXS example dataset/Analysis/Datasets per sequence')

[4]:

exp = pp.Experiment(data_per_sequence_path)

Import files: 100%|█████████████████████████████████████████████████████████████| 7948/7948 [00:00<00:00, 46126.39it/s]


Initialize experiment:
C:\Users\user\Desktop\SPARXS example dataset\Analysis\Datasets per sequence

[5]:

exp.files[0:10]

[5]:

FileCollection([File(AAAAAGCG),
                File(AAAAAGTT),
                File(AAACAAAA),
                File(AAAGCATC),
                File(AAAGGAGT),
                File(AAATCGCT),
                File(AAATCTAT),
                File(AAATTCGG),
                File(AAATTGCG),
                File(AACAACGT)])

[6]:

file_HJ1 = exp.files.select('TTAGCCGA', 'name')[0]
file_HJ3 = exp.files.select('AATCGGCT', 'name')[0]
file_HJ7 = exp.files.select('GGCGCCGC', 'name')[0]
files_HJ137 = exp.files.select('TTAGCCGA|AATCGGCT|GGCGCCGC', 'name')

Evaluate traces#

[36]:

file.show_traces(plot_variables=['intensity_total', 'intensity', 'FRET', 'classification'],
                 ylims=[(0, 35000), (0, 35000), (0, 1), (-2.5,1.5)],
                 colours=[('k'), ('g', 'r'), ('b'), ('k')], selected=False, height=5)

%matplotlib inline

Full pipeline#

[37]:

def pipeline(file, red_intensity_cutoff, intensity_total_min, intensity_total_max):
    try:
        # Define selections
        selection_Cy5_active_start = file.intensity_red_before.sel(channel=1, drop=True) > red_intensity_cutoff
        selection_Cy5_active_start.name = 'selection_Cy5_active_start'

        selection_Cy5_active_end = file.intensity_red_after.sel(channel=1, drop=True) > red_intensity_cutoff
        selection_Cy5_active_end.name = 'selection_Cy5_active_end'

        intensity_total_rolling = file.intensity_total.rolling(frame=5, center=True).mean().dropna('frame')
        selection_intensity_total = (intensity_total_rolling < intensity_total_max).all('frame')
        selection_intensity_total.name = 'selection_intensity_total'

        # Save selections
        file.set_variable(selection_Cy5_active_start)
        file.set_variable(selection_Cy5_active_end)
        file.set_variable(selection_intensity_total)

        # Apply selections
        file.apply_selections(['selection_intensity_total', 'selection_Cy5_active_start', 'selection_Cy5_active_end'])

        # Define classifictions
        classification_donor_active = xr.DataArray(medfilt((file.intensity_total > intensity_total_min).astype(int), kernel_size=(1,11)).astype(bool),
                                                   name='classification_donor_active',
                                                   dims=('molecule','frame'))

        classification_single_dye = xr.DataArray(medfilt((file.intensity_total < intensity_total_max).astype(int), kernel_size=(1,11)).astype(bool),
                                                 name='classification_single_dye',
                                                 dims=('molecule','frame'))

        # Save classifications
        file.set_variable(classification_donor_active)
        file.set_variable(classification_single_dye)

        # Apply classifications
        file.apply_classifications(classification_donor_active=-1, classification_single_dye=-2)

        # Hidden markov modelling
        file.classify_hmm('FRET')

        # Apply additional selections produced by classify_hmm
        file.apply_selections(['selection_intensity_total', 'selection_Cy5_active_start', 'selection_Cy5_active_end',
                               'selection_complex_rates', 'selection_lower_rate_limit'])

        # Apply additional classification produced by classify_hmm
        file.apply_classifications(classification_donor_active=-1, classification_single_dye=-2, classification_hmm=[None, 0, 1])

    except:
        print('Error:', file)

files_HJ137.serial.map(pipeline)(red_intensity_cutoff, intensity_total_min, intensity_total_max)

  0%|                                                                                            | 0/3 [00:00<?, ?it/s]

Serial processing

100%|█████████████████████████████████████████████████████████████████████████████| 1006/1006 [00:08<00:00, 120.60it/s]
100%|███████████████████████████████████████████████████████████████████████████████| 753/753 [00:03<00:00, 210.48it/s]
100%|███████████████████████████████████████████████████████████████████████████████| 669/669 [00:06<00:00, 102.91it/s]
100%|████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:22<00:00,  7.41s/it]

Plot results#

[38]:

files_HJ137.plot_hmm_rates()

  0%|                                                                                            | 0/3 [00:00<?, ?it/s]

Serial processing

100%|████████████████████████████████████████████████████████████████████████████████████| 3/3 [00:00<00:00,  4.00it/s]

../_images/SPARXS_5_single_molecule_data_analysis_2_62_3.png

../_images/SPARXS_5_single_molecule_data_analysis_2_62_4.png

../_images/SPARXS_5_single_molecule_data_analysis_2_62_5.png