<jupyter_start><jupyter_text>Read sample data<jupyter_code>df = pd.read_csv(os.path.join(curr_pth, 'sample_data', 'bystander_sampledata.csv'))
df<jupyter_output><empty_output><jupyter_text>The models expect sequences (i.e. target sites) to be wrapped in a `pandas.DataFrame` with a column representing the `seq_id` of a sequence and the `input sequence itself` (i.e. `protospacer sequence`). In the above example, **`Inp_seq`** represents the input sequence column and **`seq_id`** represents the sequence id column.The sequences should be of length 20 (i.e. 20 bases) and represent the protospacer target site.<jupyter_code># create a directory where we dump the predictions of the models
csv_dir = create_directory(os.path.join(curr_pth, 'sample_data', 'predictions_haplo'))<jupyter_output><empty_output><jupyter_text>Specify device (i.e. CPU or GPU) to run the models on Specify device to run the model on. The models can run on `GPU` or `CPU`. We can instantiate a device by running `get_device(to_gpu,gpu_index)` function. - To run on GPU we pass `to_gpu = True` and specify which card to use if we have multiple cards `gpu_index=int` (i.e. in case we have multiple GPU cards we specify the index counting from 0). - If there is no GPU installed, the function will return a `CPU` device. We can get a detailed information on the GPU cards installed on the compute node by calling `report_available_cuda_devices` function.<jupyter_code>report_available_cuda_devices()
# instantiate a device using the only one available :P
device = get_device(True, 0)
device<jupyter_output><empty_output><jupyter_text>Create a BE-DICT (bystander) model  **Note**: We use`bystander` and `haplotype` interchangeably in this demo. In the sample data we are working with 👇, we have `Editor` column specifying what editor to use for the corresponding sequences. We also have `Outp_seq` column representing the observed *edited outcomes* of the input sequence with its corresponding probability (i.e. proportion of occurrence) of such outcome in `true_score` column.<jupyter_code>for gr_name, gr_df in df.groupby(by=['Editor']):
    print(gr_name)
    display(gr_df)<jupyter_output>ABE8e<jupyter_text>We can run the model for a specified `base editor` by following these steps:- First instantiate `SeqProcessConfig` object which takes the following arguments:    - protospacer sequence length (equal to `20` in this work)    - start and end position of the protospacer sequence such as `(1, 20)`     - start and end position of editing window such as `(1, 20)`    - numbering used to index the positions in the sequence above such as `1` based indexing (i.e. we start counting from 1). For example, we could use `0` based indexing or any indexing (including negative number) and in this case we provide the start and end of protospacer sequence and editing window based on the chosen indexing.<jupyter_code>seqconfig_dataproc = SeqProcessConfig(20, (1,20), (1,20), 1)
print(seqconfig_dataproc)<jupyter_output>+------------------------------------------------+-------+
|           Sequence processing Config           | Value |
+------------------------------------------------+-------+
|                sequence length                 |   20  |
|    sequence start index (0-based indexing)     |   0   |
|     sequence end index (0-based indexing)      |   19  |
| editable window start index (0-based indexing) |   0   |
|  editable window end index (0-based indexing)  |   19  |
|             offset start numbering             |   1   |
|              offset end numbering              |   20  |
+------------------------------------------------+-------+<jupyter_text>Note that if we know the target editing window is in a smaller range than the whole sequence, we can specify that in the editing window start and end positions above. For example, in [Song et al. data](https://doi.org/10.1038/s41587-020-0573-5), editing window is from positions 3-10, hence we define our `SeqProcessConfig` object as follows: ```pythonseqconfig_dataproc = SeqProcessConfig(20, (1,20), (3,10), 1)``` In fact sequences in the data frame above corresponding to `ABEmax` and `BE4max` are from Song et al. data and we will use this narrower edit range when running the models.<jupyter_code>seqconfig_dataproc = SeqProcessConfig(20, (1,20), (3,10), 1)
print(seqconfig_dataproc)<jupyter_output>+------------------------------------------------+-------+
|           Sequence processing Config           | Value |
+------------------------------------------------+-------+
|                sequence length                 |   20  |
|    sequence start index (0-based indexing)     |   0   |
|     sequence end index (0-based indexing)      |   19  |
| editable window start index (0-based indexing) |   2   |
|  editable window end index (0-based indexing)  |   9   |
|             offset start numbering             |   1   |
|              offset end numbering              |   20  |
+------------------------------------------------+-------+<jupyter_text>As you can see above ☝️, independent of the user provided indexing (see `offset start numbering`), the `SeqProcessConfig` will generate `0`-based internal indexing :) In other words, the supplied indexing is a `user friendly` way to refer to positions in the sequence (i.e. edit window and start of protospacer sequence) without doing the conversion on your own.    - Second we define `HaplotypeSeqProcessor` object which takes the following arguments:      - target editor name as in {`ABEmax`, `BE4max`, `ABE8e`, `Target-AID`}      - target canonical nucleotide and its expected conversion depending on the base editor used. For example,          - `(A,G)` for editors `ABEmax` and `ABE8e`         - `(C,T)` for editors `BE4max` and `Target-AID`.      - `SeqProcessConfig` instantiated above.<jupyter_code>teditor = 'ABEmax'
cnv_nucl = ('A', 'G')
seq_processor = HaplotypeSeqProcessor(teditor, cnv_nucl, seqconfig_dataproc)<jupyter_output>+--------------------------------------+--------+
|             Description              | Value  |
+--------------------------------------+--------+
|             Base editor              | ABEmax |
|          Target nucleotide           |   A    |
|        Conversion nucleotide         |   G    |
| Maximum number of targets considered |   12   |
+--------------------------------------+--------+
+------------------------------------------------+-------+
|           Sequence processing Config           | Value |
+------------------------------------------------+-------+
|                sequence length                 |   20  |
|    sequence start index (0-based indexing)     |   0   |
|     sequence end index (0-based indexing)      |   19  |
| editable window start index (0-based indexing) |   2   |
|  editable window end index (0-based indexing)  |   9   |
|             offset start numbering             |   1   |
|              offset end numbering              |   20  |
+---------[...]<jupyter_text>- Third we instantiate `SeqProcessConfig` similar to first step, to be used in generating data tensors      - protospacer sequence length (equal to `20` in this work)      - start and end position of the sequence such as `(1, 20)`       - start and end position of editing window such as `(1, 20)`      - numbering used to index the positions in the sequence above such as `1` based indexing (i.e. we start counting from 1).<jupyter_code>seqconfig_datatensgen = SeqProcessConfig(20, (1,20), (1 ,20), 1)<jupyter_output><empty_output><jupyter_text>The `SeqProcessConfig` used in the third step should mirror the object defined in first step except the `editing window` end position is always fixed to `20` (i.e. protospacer sequence length). We can change the `start` of `editing window` (i.e. using `3` instead of `1` given that we know the editing happens in the range `3-10`) or simply we specify equivalent start and end positions for both `protospacer sequence` and the `editing window` such as:```pythonseqconfig_datatensgen = SeqProcessConfig(20, (1,20), (1 ,20), 1) ``` Now we reach to the last step where we define the model `BEDICT_HaplotypeModel` and supply the objects defined above:   - `seq_processor` which is a `HaplotypeSeqProcessor` object defined in step 2.   - `seqconfig_datatensgen` which is a `SeqProcessConfig` object defined in step 3.   - `device` (i.e. CPU or GPU)<jupyter_code>bedict_haplo = BEDICT_HaplotypeModel(seq_processor, seqconfig_datatensgen, device)<jupyter_output><empty_output><jupyter_text>Run a BE-DICT (bystander) model  Given that we defined our model to run using `ABEmax` editor we will focus on sequences in our loaded data fame that have `Editor` column equal to `ABEmax`.<jupyter_code>sample_df = df.loc[df['Editor'] == 'ABEmax'].copy()
sample_df<jupyter_output><empty_output><jupyter_text>First, we prepare the data by running `prepare_data` method that takes the following args:   - sample data frame such as `sample_df`   - input sequence column name such as `Inp_seq`   - output sequence column name (optional) such as `Outp_seq`   - outcome column name (optional) representing probability of `output sequence`   - batch size representing the number of sequences in the batch  Main (default) mode The are different modes (i.e. options) to run the model. In the main (default) mode, we just provide the data frame that includes the input sequences and the model will generate the different outcome sequences and their corresponding probability.<jupyter_code>dloader = bedict_haplo.prepare_data(sample_df,
                                    ['seq_id','Inp_seq'],
                                    outpseq_col=None,
                                    outcome_col=None,
                                    renormalize=False,
                                    batch_size=500)<jupyter_output>225it [00:00, 18109.77it/s]            
2it [00:00, 146.11it/s]
250it [00:00, 11777.78it/s]            
100%|██████████| 2/2 [00:00<00:00, 211.34it/s]<jupyter_text>After we parse the data, we can run the model to get predictions. We have trained `5` different models for every `base editor` (denoted by runs) which we will use on this data and take their average predictions. The `trained` bystander/haplotype models are found under the `trained_models/bystander` directory.trained_models├── bystander│   ├── ABE8e│   │   └── train_val│   │       ├── run_0│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_1│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_2│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_3│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       └── run_4│   │           ├── config│   │           │   ├── exp_options.pkl│   │           │   └── mconfig.pkl│   │           └── model_statedict│   │               └── HaplotypeTransformer.pkl│   ├── ABEmax│   │   └── train_val│   │       ├── run_0│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_1│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_2│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_3│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       └── run_4│   │           ├── config│   │           │   ├── exp_options.pkl│   │           │   └── mconfig.pkl│   │           └── model_statedict│   │               └── HaplotypeTransformer.pkl│   ├── BE4max│   │   └── train_val│   │       ├── run_0│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_1│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_2│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       ├── run_3│   │       │   ├── config│   │       │   │   ├── exp_options.pkl│   │       │   │   └── mconfig.pkl│   │       │   └── model_statedict│   │       │       └── HaplotypeTransformer.pkl│   │       └── run_4│   │           ├── config│   │           │   ├── exp_options.pkl│   │           │   └── mconfig.pkl│   │           └── model_statedict│   │               └── HaplotypeTransformer.pkl│   └── Target-AID│       └── train_val│           ├── run_0│           │   ├── config│           │   │   ├── exp_options.pkl│           │   │   └── mconfig.pkl│           │   └── model_statedict│           │       └── HaplotypeTransformer.pkl│           ├── run_1│           │   ├── config│           │   │   ├── exp_options.pkl│           │   │   └── mconfig.pkl│           │   └── model_statedict│           │       └── HaplotypeTransformer.pkl│           ├── run_2│           │   ├── config│           │   │   ├── exp_options.pkl│           │   │   └── mconfig.pkl│           │   └── model_statedict│           │       └── HaplotypeTransformer.pkl│           ├── run_3│           │   ├── config│           │   │   ├── exp_options.pkl│           │   │   └── mconfig.pkl│           │   └── model_statedict│           │       └── HaplotypeTransformer.pkl│           └── run_4│               ├── config│               │   ├── exp_options.pkl│               │   └── mconfig.pkl│               └── model_statedict│                   └── HaplotypeTransformer.pkl|└── perbase    ├── ABE8e    │   └── train_val    │       ├── run_0    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_1    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_2    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_3    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       └── run_4    │           ├── config    │           │   ├── exp_options.pkl    │           │   └── mconfig.pkl    │           └── model_statedict    │               └── Transformer.pkl    ├── ABEmax    │   └── train_val    │       ├── run_0    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_1    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_2    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_3    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       └── run_4    │           ├── config    │           │   ├── exp_options.pkl    │           │   └── mconfig.pkl    │           └── model_statedict    │               └── Transformer.pkl    ├── BE4max    │   └── train_val    │       ├── run_0    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_1    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_2    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       ├── run_3    │       │   ├── config    │       │   │   ├── exp_options.pkl    │       │   │   └── mconfig.pkl    │       │   └── model_statedict    │       │       └── Transformer.pkl    │       └── run_4    │           ├── config    │           │   ├── exp_options.pkl    │           │   └── mconfig.pkl    │           └── model_statedict    │               └── Transformer.pkl    └── Target-AID        └── train_val            ├── run_0            │   ├── config            │   │   ├── exp_options.pkl            │   │   └── mconfig.pkl            │   └── model_statedict            │       └── Transformer.pkl            ├── run_1            │   ├── config            │   │   ├── exp_options.pkl            │   │   └── mconfig.pkl            │   └── model_statedict            │       └── Transformer.pkl            ├── run_2            │   ├── config            │   │   ├── exp_options.pkl            │   │   └── mconfig.pkl            │   └── model_statedict            │       └── Transformer.pkl            ├── run_3            │   ├── config            │   │   ├── exp_options.pkl            │   │   └── mconfig.pkl            │   └── model_statedict            │       └── Transformer.pkl            └── run_4                ├── config                │   ├── exp_options.pkl                │   └── mconfig.pkl                └── model_statedict                    └── Transformer.pkl     We run the model to get prediction using `predict_from_dloader` method that takes the following args:   - data loader such as `dloader` we generated above   - path to trained model such as `crispr/trained_models/ABEmax/train_val/run_0`   - outcome column name (optional) representing `true probability` of the output sequence. It is typically used if we need to keep the column in the generated prediction data frame to be later used for performance computation.  Moreover, after we generate `5` predictions we can average them using `compute_avg_predictions` method that takes a concatenated data frame from the predictions of each run.<jupyter_code>num_runs = 5 # number of models
pred_df_lst = [] # list to hold the prediction data frames for each model run
for run_num in range(num_runs):
    # specify model directory based on the chose editor
    model_dir = os.path.join(curr_pth, 'trained_models', 'bystander', teditor, 'train_val', f'run_{run_num}')
    print('run_num:', run_num)
    print('model_dir:', model_dir)

    pred_df = bedict_haplo.predict_from_dloader(dloader,
                                                model_dir, 
                                                outcome_col=None)
    pred_df['run_num'] = run_num
    pred_df_lst.append(pred_df)
    
# aggregate predictions in one data frame
pred_df_unif = pd.concat(pred_df_lst, axis=0, ignore_index=True)
# sanity check
check_na(pred_df_unif)

# compute average predictions
agg_df = bedict_haplo.compute_avg_predictions(pred_df_unif)
# sanity check
check_na(agg_df)

agg_df<jupyter_output><empty_output><jupyter_text>Visualize results<jupyter_code>from haplotype.data_preprocess import get_char
from IPython.core.display import HTML
from haplotype.data_preprocess import VizInpOutp_Haplotype<jupyter_output><empty_output><jupyter_text>To visualize the results, we use `visualize_haplotype` method that takes the following args:   - genereated predictions from the model (such as `agg_df`)   - list of sequence ids to visualize (we can pass all sequences such as `sample_df['seq_id']`)   - sequence id and input sequence column names in the generated `agg_df` (default is [`seq_id`,`Inp_seq`])   - output sequence column name in the generated `agg_df`(default is `Outp_seq`)   - prediction score column name in the generated `agg_df` (default is `pred_score`)   - prediction score threshold that filters out output sequences with probability below the specified threshold. If all output sequences are to be plotted then `predscore_thr` should be equal to `0.`.<jupyter_code>tseqids = sample_df['seq_id']
res_html = bedict_haplo.visualize_haplotype(agg_df, 
                                            tseqids, 
                                            ['seq_id','Inp_seq'], 
                                            'Outp_seq', 
                                            'pred_score', 
                                            predscore_thr=0.)
for seq_id in res_html:
    display(HTML(res_html[seq_id]))<jupyter_output>250it [00:00, 10263.15it/s]            
100%|██████████| 2/2 [00:00<00:00, 116.41it/s]<jupyter_text>`res_html` is a dictionary that has `input sequences` as `keys` and the corresponding `value` is the generated `html table`. We can also save the visualization of each sequence (or all sequences) in html format on disk using `HaplotypeVizFile` object.<jupyter_code>from haplotype.data_preprocess import HaplotypeVizFile<jupyter_output><empty_output><jupyter_text>`HaplotypeVizFile` object takes the path to `viz resources` needed to generate the html webpages. By default it is under `haplotype/viz_resources` folder.```pythonHaplotypeVizFile(os.path.join(curr_pth, 'haplotype', 'viz_resources'))``` To create and save the files on disk, we use `create` method that takes three arguments:   - generated `html table` from `visualize_haplotype` method   - path where to save the file   - name of the file (without extension)<jupyter_code>vf = HaplotypeVizFile(os.path.join(curr_pth, 'haplotype', 'viz_resources'))
for seq_id in res_html:
    display(HTML(res_html[seq_id]))
    vf.create(res_html[seq_id], csv_dir, f'{teditor}_{seq_id}_haplotype')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14184801684801684, 'fraction_numerical': 0.013030888030888031, 'mean_word_length': 1.9304448444330162, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 59, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '124419', 'license': 'MIT', 'n_tokens_mistral': 8263, 'n_tokens_neox': 7743, 'n_words': 2455, 'path': 'demo/BEDICT_haplotype_model_demoseqs.ipynb', 'repo_name': 'zxChouSean/crispr_bedict_reproduce', 'size': '22792'}