Spaces:
Runtime error
Runtime error
a96123155
commited on
Commit
·
c2ce74d
1
Parent(s):
dedea73
app
Browse files
.DS_Store
CHANGED
|
Binary files a/.DS_Store and b/.DS_Store differ
|
|
|
app.py
ADDED
|
@@ -0,0 +1,550 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
import streamlit as st
|
| 2 |
+
from io import StringIO
|
| 3 |
+
from Bio import SeqIO
|
| 4 |
+
|
| 5 |
+
st.title("IRES-LM prediction and mutation")
|
| 6 |
+
|
| 7 |
+
# Input sequence
|
| 8 |
+
st.subheader("Input sequence")
|
| 9 |
+
|
| 10 |
+
seq = st.text_area("FASTA format only", value=">vir_CVB3_ires_00505.1\nTTAAAACAGCCTGTGGGTTGATCCCACCCACAGGCCCATTGGGCGCTAGCACTCTGGTATCACGGTACCTTTGTGCGCCTGTTTTATACCCCCTCCCCCAACTGTAACTTAGAAGTAACACACACCGATCAACAGTCAGCGTGGCACACCAGCCACGTTTTGATCAAGCACTTCTGTTACCCCGGACTGAGTATCAATAGACTGCTCACGCGGTTGAAGGAGAAAGCGTTCGTTATCCGGCCAACTACTTCGAAAAACCTAGTAACACCGTGGAAGTTGCAGAGTGTTTCGCTCAGCACTACCCCAGTGTAGATCAGGTCGATGAGTCACCGCATTCCCCACGGGCGACCGTGGCGGTGGCTGCGTTGGCGGCCTGCCCATGGGGAAACCCATGGGACGCTCTAATACAGACATGGTGCGAAGAGTCTATTGAGCTAGTTGGTAGTCCTCCGGCCCCTGAATGCGGCTAATCCTAACTGCGGAGCACACACCCTCAAGCCAGAGGGCAGTGTGTCGTAACGGGCAACTCTGCAGCGGAACCGACTACTTTGGGTGTCCGTGTTTCATTTTATTCCTATACTGGCTGCTTATGGTGACAATTGAGAGATCGTTACCATATAGCTATTGGATTGGCCATCCGGTGACTAATAGAGCTATTATATATCCCTTTGTTGGGTTTATACCACTTAGCTTGAAAGAGGTTAAAACATTACAATTCATTGTTAAGTTGAATACAGCAAA")
|
| 11 |
+
st.subheader("Upload sequence file")
|
| 12 |
+
uploaded = st.file_uploader("Sequence file in FASTA format")
|
| 13 |
+
|
| 14 |
+
# augments
|
| 15 |
+
global output_filename, start_nt_position, end_nt_position, mut_by_prob, transform_type, mlm_tok_num, n_mut, n_designs_ep, n_sampling_designs_ep, n_mlm_recovery_sampling, mutate2stronger
|
| 16 |
+
output_filename = st.text_input("output a .csv file", value='IRES_LM_prediction_mutation')
|
| 17 |
+
start_nt_position = st.number_input("The start position of the mutation of this sequence, the first position is defined as 0", value=0)
|
| 18 |
+
end_nt_position = st.number_input("The last position of the mutation of this sequence, the last position is defined as length(sequence)-1 or -1", value=-1)
|
| 19 |
+
mut_by_prob = st.checkbox("Mutated by predicted Probability or Transformed Probability of the sequence", value=True)
|
| 20 |
+
transform_type = st.selectbox("Type of probability transformation",
|
| 21 |
+
['', 'sigmoid', 'logit', 'power_law', 'tanh'],
|
| 22 |
+
index=2)
|
| 23 |
+
mlm_tok_num = st.number_input("Number of masked tokens for each sequence per epoch", value=1)
|
| 24 |
+
n_mut = st.number_input("Maximum number of mutations for each sequence", value=3)
|
| 25 |
+
n_designs_ep = st.number_input("Number of mutations per epoch", value=10)
|
| 26 |
+
n_sampling_designs_ep = st.number_input("Number of sampling mutations from n_designs_ep per epoch", value=5)
|
| 27 |
+
n_mlm_recovery_sampling = st.number_input("Number of MLM recovery samplings (with AGCT recovery)", value=1)
|
| 28 |
+
mutate2stronger = st.checkbox("Mutate to stronger IRES variant, otherwise mutate to weaker IRES", value=True)
|
| 29 |
+
|
| 30 |
+
if not mut_by_prob and transform_type != '':
|
| 31 |
+
st.write("--transform_type must be '' when --mut_by_prob is False")
|
| 32 |
+
transform_type = ''
|
| 33 |
+
|
| 34 |
+
|
| 35 |
+
# Import necessary libraries
|
| 36 |
+
# import matplotlib
|
| 37 |
+
# import matplotlib.pyplot as plt
|
| 38 |
+
import numpy as np
|
| 39 |
+
import os
|
| 40 |
+
import pandas as pd
|
| 41 |
+
# import pathlib
|
| 42 |
+
import random
|
| 43 |
+
# import scanpy as sc
|
| 44 |
+
# import seaborn as sns
|
| 45 |
+
import torch
|
| 46 |
+
import torch.nn as nn
|
| 47 |
+
import torch.nn.functional as F
|
| 48 |
+
# from argparse import Namespace
|
| 49 |
+
from collections import Counter, OrderedDict
|
| 50 |
+
from copy import deepcopy
|
| 51 |
+
from esm import Alphabet, FastaBatchedDataset, ProteinBertModel, pretrained, MSATransformer
|
| 52 |
+
from esm.data import *
|
| 53 |
+
from esm.model.esm2 import ESM2
|
| 54 |
+
# from sklearn import preprocessing
|
| 55 |
+
# from sklearn.metrics import (confusion_matrix, roc_auc_score, auc,
|
| 56 |
+
# precision_recall_fscore_support,
|
| 57 |
+
# precision_recall_curve, classification_report,
|
| 58 |
+
# roc_auc_score, average_precision_score,
|
| 59 |
+
# precision_score, recall_score, f1_score,
|
| 60 |
+
# accuracy_score)
|
| 61 |
+
# from sklearn.model_selection import StratifiedKFold
|
| 62 |
+
# from sklearn.utils import class_weight
|
| 63 |
+
# from scipy.stats import spearmanr, pearsonr
|
| 64 |
+
from torch import nn
|
| 65 |
+
from torch.nn import Linear
|
| 66 |
+
from torch.nn.utils.rnn import pad_sequence
|
| 67 |
+
from torch.utils.data import Dataset, DataLoader
|
| 68 |
+
from tqdm import tqdm, trange
|
| 69 |
+
|
| 70 |
+
# Set global variables
|
| 71 |
+
# matplotlib.rcParams.update({'font.size': 7})
|
| 72 |
+
seed = 19961231
|
| 73 |
+
random.seed(seed)
|
| 74 |
+
np.random.seed(seed)
|
| 75 |
+
torch.manual_seed(seed)
|
| 76 |
+
# torch.cuda.manual_seed(seed)
|
| 77 |
+
# torch.backends.cudnn.deterministic = True
|
| 78 |
+
# torch.backends.cudnn.benchmark = False
|
| 79 |
+
|
| 80 |
+
|
| 81 |
+
global idx_to_tok, prefix, epochs, layers, heads, fc_node, dropout_prob, embed_dim, batch_toks, device, repr_layers, evaluation, include, truncate, return_contacts, return_representation, mask_toks_id, finetune
|
| 82 |
+
|
| 83 |
+
epochs = 5
|
| 84 |
+
layers = 6
|
| 85 |
+
heads = 16
|
| 86 |
+
embed_dim = 128
|
| 87 |
+
batch_toks = 4096
|
| 88 |
+
fc_node = 64
|
| 89 |
+
dropout_prob = 0.5
|
| 90 |
+
folds = 10
|
| 91 |
+
repr_layers = [-1]
|
| 92 |
+
include = ["mean"]
|
| 93 |
+
truncate = True
|
| 94 |
+
finetune = False
|
| 95 |
+
return_contacts = False
|
| 96 |
+
return_representation = False
|
| 97 |
+
|
| 98 |
+
device = "cpu"
|
| 99 |
+
|
| 100 |
+
global tok_to_idx, idx_to_tok, mask_toks_id
|
| 101 |
+
alphabet = Alphabet(mask_prob = 0.15, standard_toks = 'AGCT')
|
| 102 |
+
assert alphabet.tok_to_idx == {'<pad>': 0, '<eos>': 1, '<unk>': 2, 'A': 3, 'G': 4, 'C': 5, 'T': 6, '<cls>': 7, '<mask>': 8, '<sep>': 9}
|
| 103 |
+
|
| 104 |
+
# tok_to_idx = {'<pad>': 0, '<eos>': 1, '<unk>': 2, 'A': 3, 'G': 4, 'C': 5, 'T': 6, '<cls>': 7, '<mask>': 8, '<sep>': 9}
|
| 105 |
+
tok_to_idx = {'-': 0, '&': 1, '?': 2, 'A': 3, 'G': 4, 'C': 5, 'T': 6, '!': 7, '*': 8, '|': 9}
|
| 106 |
+
idx_to_tok = {idx: tok for tok, idx in tok_to_idx.items()}
|
| 107 |
+
# st.write(tok_to_idx)
|
| 108 |
+
mask_toks_id = 8
|
| 109 |
+
|
| 110 |
+
global w1, w2, w3
|
| 111 |
+
w1, w2, w3 = 1, 1, 100
|
| 112 |
+
|
| 113 |
+
class CNN_linear(nn.Module):
|
| 114 |
+
def __init__(self):
|
| 115 |
+
super(CNN_linear, self).__init__()
|
| 116 |
+
|
| 117 |
+
self.esm2 = ESM2(num_layers = layers,
|
| 118 |
+
embed_dim = embed_dim,
|
| 119 |
+
attention_heads = heads,
|
| 120 |
+
alphabet = alphabet)
|
| 121 |
+
|
| 122 |
+
self.dropout = nn.Dropout(dropout_prob)
|
| 123 |
+
self.relu = nn.ReLU()
|
| 124 |
+
self.flatten = nn.Flatten()
|
| 125 |
+
self.fc = nn.Linear(in_features = embed_dim, out_features = fc_node)
|
| 126 |
+
self.output = nn.Linear(in_features = fc_node, out_features = 2)
|
| 127 |
+
|
| 128 |
+
def predict(self, tokens):
|
| 129 |
+
|
| 130 |
+
x = self.esm2(tokens, [layers], need_head_weights=False, return_contacts=False, return_representation = True)
|
| 131 |
+
x_cls = x["representations"][layers][:, 0]
|
| 132 |
+
|
| 133 |
+
o = self.fc(x_cls)
|
| 134 |
+
o = self.relu(o)
|
| 135 |
+
o = self.dropout(o)
|
| 136 |
+
o = self.output(o)
|
| 137 |
+
|
| 138 |
+
y_prob = torch.softmax(o, dim = 1)
|
| 139 |
+
y_pred = torch.argmax(y_prob, dim = 1)
|
| 140 |
+
|
| 141 |
+
if transform_type:
|
| 142 |
+
y_prob_transformed = prob_transform(y_prob[:,1])
|
| 143 |
+
return y_prob[:,1], y_pred, x['logits'], y_prob_transformed
|
| 144 |
+
else:
|
| 145 |
+
return y_prob[:,1], y_pred, x['logits'], o[:,1]
|
| 146 |
+
|
| 147 |
+
def forward(self, x1, x2):
|
| 148 |
+
logit_1, repr_1 = self.predict(x1)
|
| 149 |
+
logit_2, repr_2 = self.predict(x2)
|
| 150 |
+
return (logit_1, logit_2), (repr_1, repr_2)
|
| 151 |
+
|
| 152 |
+
def prob_transform(prob, **kwargs): # Logits
|
| 153 |
+
"""
|
| 154 |
+
Transforms probability values based on the specified method.
|
| 155 |
+
|
| 156 |
+
:param prob: torch.Tensor, the input probabilities to be transformed
|
| 157 |
+
:param transform_type: str, the type of transformation to be applied
|
| 158 |
+
:param kwargs: additional parameters for transformations
|
| 159 |
+
:return: torch.Tensor, transformed probabilities
|
| 160 |
+
"""
|
| 161 |
+
|
| 162 |
+
if transform_type == 'sigmoid':
|
| 163 |
+
x0 = kwget('x0', 0.5)
|
| 164 |
+
k = kwget('k', 10.0)
|
| 165 |
+
prob_transformed = 1 / (1 + torch.exp(-k * (prob - x0)))
|
| 166 |
+
|
| 167 |
+
elif transform_type == 'logit':
|
| 168 |
+
# Adding a small value to avoid log(0) and log(1)
|
| 169 |
+
prob_transformed = torch.log(prob + 1e-6) - torch.log(1 - prob + 1e-6)
|
| 170 |
+
|
| 171 |
+
elif transform_type == 'power_law':
|
| 172 |
+
gamma = kwget('gamma', 2.0)
|
| 173 |
+
prob_transformed = torch.pow(prob, gamma)
|
| 174 |
+
|
| 175 |
+
elif transform_type == 'tanh':
|
| 176 |
+
k = kwget('k', 2.0)
|
| 177 |
+
prob_transformed = torch.tanh(k * prob)
|
| 178 |
+
|
| 179 |
+
return prob_transformed
|
| 180 |
+
|
| 181 |
+
def random_replace(sequence, continuous_replace=False):
|
| 182 |
+
if end_nt_position == -1: end_nt_position = len(sequence)
|
| 183 |
+
if start_nt_position < 0 or end_nt_position > len(sequence) or start_nt_position > end_nt_position:
|
| 184 |
+
# raise ValueError("Invalid start/end positions")
|
| 185 |
+
st.write("Invalid start/end positions")
|
| 186 |
+
start_nt_position, end_nt_position = 0, -1
|
| 187 |
+
|
| 188 |
+
# 将序列切片成三部分:替换区域前、替换区域、替换区域后
|
| 189 |
+
pre_segment = sequence[:start_nt_position]
|
| 190 |
+
target_segment = list(sequence[start_nt_position:end_nt_position + 1]) # +1因为Python的切片是右开区间
|
| 191 |
+
post_segment = sequence[end_nt_position + 1:]
|
| 192 |
+
|
| 193 |
+
if not continuous_replace:
|
| 194 |
+
# 随机替换目标片段的mlm_tok_num个位置
|
| 195 |
+
indices = random.sample(range(len(target_segment)), mlm_tok_num)
|
| 196 |
+
for idx in indices:
|
| 197 |
+
target_segment[idx] = '*'
|
| 198 |
+
else:
|
| 199 |
+
# 在目标片段连续替换mlm_tok_num个位置
|
| 200 |
+
max_start_idx = len(target_segment) - mlm_tok_num # 确保从i开始的n_mut个元素不会超出目标片段的长度
|
| 201 |
+
if max_start_idx < 1: # 如果目标片段长度小于mlm_tok_num,返回原始序列
|
| 202 |
+
return target_segment
|
| 203 |
+
start_idx = random.randint(0, max_start_idx)
|
| 204 |
+
for idx in range(start_idx, start_idx + mlm_tok_num):
|
| 205 |
+
target_segment[idx] = '*'
|
| 206 |
+
|
| 207 |
+
# 合并并返回最终的序列
|
| 208 |
+
return ''.join([pre_segment] + target_segment + [post_segment])
|
| 209 |
+
|
| 210 |
+
|
| 211 |
+
def mlm_seq(seq):
|
| 212 |
+
seq_token, masked_sequence_token = [7],[7]
|
| 213 |
+
seq_token += [tok_to_idx[token] for token in seq]
|
| 214 |
+
|
| 215 |
+
masked_seq = random_replace(seq, n_mut) # 随机替换n_mut个元素为'*'
|
| 216 |
+
masked_seq_token += [tok_to_idx[token] for token in masked_seq]
|
| 217 |
+
|
| 218 |
+
return seq, masked_seq, torch.LongTensor(seq_token), torch.LongTensor(masked_seq_token)
|
| 219 |
+
|
| 220 |
+
def batch_mlm_seq(seq_list, continuous_replace = False):
|
| 221 |
+
batch_seq = []
|
| 222 |
+
batch_masked_seq = []
|
| 223 |
+
batch_seq_token_list = []
|
| 224 |
+
batch_masked_seq_token_list = []
|
| 225 |
+
|
| 226 |
+
for i, seq in enumerate(seq_list):
|
| 227 |
+
seq_token, masked_seq_token = [7], [7]
|
| 228 |
+
seq_token += [tok_to_idx[token] for token in seq]
|
| 229 |
+
|
| 230 |
+
masked_seq = random_replace(seq, continuous_replace) # 随机��换n_mut个元素为'*'
|
| 231 |
+
masked_seq_token += [tok_to_idx[token] for token in masked_seq]
|
| 232 |
+
|
| 233 |
+
batch_seq.append(seq)
|
| 234 |
+
batch_masked_seq.append(masked_seq)
|
| 235 |
+
|
| 236 |
+
batch_seq_token_list.append(seq_token)
|
| 237 |
+
batch_masked_seq_token_list.append(masked_seq_token)
|
| 238 |
+
|
| 239 |
+
return batch_seq, batch_masked_seq, torch.LongTensor(batch_seq_token_list), torch.LongTensor(batch_masked_seq_token_list)
|
| 240 |
+
|
| 241 |
+
def recovered_mlm_tokens(masked_seqs, masked_toks, esm_logits, exclude_low_prob = False):
|
| 242 |
+
# Only remain the AGCT logits
|
| 243 |
+
esm_logits = esm_logits[:,:,3:7]
|
| 244 |
+
# Get the predicted tokens using argmax
|
| 245 |
+
predicted_toks = (esm_logits.argmax(dim=-1)+3).tolist()
|
| 246 |
+
|
| 247 |
+
batch_size, seq_len, vocab_size = esm_logits.size()
|
| 248 |
+
if exclude_low_prob: min_prob = 1 / vocab_size
|
| 249 |
+
# Initialize an empty list to store the recovered sequences
|
| 250 |
+
recovered_sequences, recovered_toks = [], []
|
| 251 |
+
|
| 252 |
+
for i in range(batch_size):
|
| 253 |
+
recovered_sequence_i, recovered_tok_i = [], []
|
| 254 |
+
for j in range(seq_len):
|
| 255 |
+
if masked_toks[i][j] == 8:
|
| 256 |
+
st.write(i,j)
|
| 257 |
+
### Sample M recovery sequences using the logits
|
| 258 |
+
recovery_probs = torch.softmax(esm_logits[i, j], dim=-1)
|
| 259 |
+
recovery_probs[predicted_toks[i][j]-3] = 0 # Exclude the most probable token
|
| 260 |
+
if exclude_low_prob: recovery_probs[recovery_probs < min_prob] = 0 # Exclude tokens with low probs < min_prob
|
| 261 |
+
recovery_probs /= recovery_probs.sum() # Normalize the probabilities
|
| 262 |
+
|
| 263 |
+
### 有放回抽样
|
| 264 |
+
max_retries = 5
|
| 265 |
+
retries = 0
|
| 266 |
+
success = False
|
| 267 |
+
|
| 268 |
+
while retries < max_retries and not success:
|
| 269 |
+
try:
|
| 270 |
+
recovery_indices = list(np.random.choice(vocab_size, size=n_mlm_recovery_sampling, p=recovery_probs.cpu().detach().numpy(), replace=False))
|
| 271 |
+
success = True # 设置成功标志
|
| 272 |
+
except ValueError as e:
|
| 273 |
+
retries += 1
|
| 274 |
+
st.write(f"Attempt {retries} failed with error: {e}")
|
| 275 |
+
if retries >= max_retries:
|
| 276 |
+
st.write("Max retries reached. Skipping this iteration.")
|
| 277 |
+
|
| 278 |
+
### recovery to sequence
|
| 279 |
+
if retries < max_retries:
|
| 280 |
+
for idx in [predicted_toks[i][j]] + [3+i for i in recovery_indices]:
|
| 281 |
+
recovery_seq = deepcopy(list(masked_seqs[i]))
|
| 282 |
+
recovery_tok = deepcopy(masked_toks[i])
|
| 283 |
+
|
| 284 |
+
recovery_tok[j] = idx
|
| 285 |
+
recovery_seq[j-1] = idx_to_tok[idx]
|
| 286 |
+
|
| 287 |
+
recovered_tok_i.append(recovery_tok)
|
| 288 |
+
recovered_sequence_i.append(''.join(recovery_seq))
|
| 289 |
+
|
| 290 |
+
recovered_sequences.extend(recovered_sequence_i)
|
| 291 |
+
recovered_toks.extend(recovered_tok_i)
|
| 292 |
+
return recovered_sequences, torch.LongTensor(torch.stack(recovered_toks))
|
| 293 |
+
|
| 294 |
+
def recovered_mlm_multi_tokens(masked_seqs, masked_toks, esm_logits, exclude_low_prob = False):
|
| 295 |
+
# Only remain the AGCT logits
|
| 296 |
+
esm_logits = esm_logits[:,:,3:7]
|
| 297 |
+
# Get the predicted tokens using argmax
|
| 298 |
+
predicted_toks = (esm_logits.argmax(dim=-1)+3).tolist()
|
| 299 |
+
|
| 300 |
+
batch_size, seq_len, vocab_size = esm_logits.size()
|
| 301 |
+
if exclude_low_prob: min_prob = 1 / vocab_size
|
| 302 |
+
# Initialize an empty list to store the recovered sequences
|
| 303 |
+
recovered_sequences, recovered_toks = [], []
|
| 304 |
+
|
| 305 |
+
for i in range(batch_size):
|
| 306 |
+
recovered_sequence_i, recovered_tok_i = [], []
|
| 307 |
+
recovered_masked_num = 0
|
| 308 |
+
for j in range(seq_len):
|
| 309 |
+
if masked_toks[i][j] == 8:
|
| 310 |
+
### Sample M recovery sequences using the logits
|
| 311 |
+
recovery_probs = torch.softmax(esm_logits[i, j], dim=-1)
|
| 312 |
+
recovery_probs[predicted_toks[i][j]-3] = 0 # Exclude the most probable token
|
| 313 |
+
if exclude_low_prob: recovery_probs[recovery_probs < min_prob] = 0 # Exclude tokens with low probs < min_prob
|
| 314 |
+
recovery_probs /= recovery_probs.sum() # Normalize the probabilities
|
| 315 |
+
|
| 316 |
+
### 有放回抽样
|
| 317 |
+
max_retries = 5
|
| 318 |
+
retries = 0
|
| 319 |
+
success = False
|
| 320 |
+
|
| 321 |
+
while retries < max_retries and not success:
|
| 322 |
+
try:
|
| 323 |
+
recovery_indices = list(np.random.choice(vocab_size, size=n_mlm_recovery_sampling, p=recovery_probs.cpu().detach().numpy(), replace=False))
|
| 324 |
+
success = True # 设置成功标志
|
| 325 |
+
except ValueError as e:
|
| 326 |
+
retries += 1
|
| 327 |
+
st.write(f"Attempt {retries} failed with error: {e}")
|
| 328 |
+
if retries >= max_retries:
|
| 329 |
+
st.write("Max retries reached. Skipping this iteration.")
|
| 330 |
+
|
| 331 |
+
### recovery to sequence
|
| 332 |
+
|
| 333 |
+
if recovered_masked_num == 0:
|
| 334 |
+
if retries < max_retries:
|
| 335 |
+
for idx in [predicted_toks[i][j]] + [3+i for i in recovery_indices]:
|
| 336 |
+
recovery_seq = deepcopy(list(masked_seqs[i]))
|
| 337 |
+
recovery_tok = deepcopy(masked_toks[i])
|
| 338 |
+
|
| 339 |
+
recovery_tok[j] = idx
|
| 340 |
+
recovery_seq[j-1] = idx_to_tok[idx]
|
| 341 |
+
|
| 342 |
+
recovered_tok_i.append(recovery_tok)
|
| 343 |
+
recovered_sequence_i.append(''.join(recovery_seq))
|
| 344 |
+
|
| 345 |
+
elif recovered_masked_num > 0:
|
| 346 |
+
if retries < max_retries:
|
| 347 |
+
for idx in [predicted_toks[i][j]] + [3+i for i in recovery_indices]:
|
| 348 |
+
for recovery_seq, recovery_tok in zip(list(recovered_sequence_i), list(recovered_tok_i)): # 要在循环开始之前获取列表的副本来进行迭代。这样,在循环中即使我们修改了原始的列表,也不会影响迭代的行为。
|
| 349 |
+
|
| 350 |
+
recovery_seq_temp = list(recovery_seq)
|
| 351 |
+
recovery_tok[j] = idx
|
| 352 |
+
recovery_seq_temp[j-1] = idx_to_tok[idx]
|
| 353 |
+
|
| 354 |
+
recovered_tok_i.append(recovery_tok)
|
| 355 |
+
recovered_sequence_i.append(''.join(recovery_seq_temp))
|
| 356 |
+
|
| 357 |
+
recovered_masked_num += 1
|
| 358 |
+
recovered_indices = [i for i, s in enumerate(recovered_sequence_i) if '*' not in s]
|
| 359 |
+
recovered_tok_i = [recovered_tok_i[i] for i in recovered_indices]
|
| 360 |
+
recovered_sequence_i = [recovered_sequence_i[i] for i in recovered_indices]
|
| 361 |
+
|
| 362 |
+
recovered_sequences.extend(recovered_sequence_i)
|
| 363 |
+
recovered_toks.extend(recovered_tok_i)
|
| 364 |
+
|
| 365 |
+
recovered_sequences, recovered_toks = remove_duplicates_double(recovered_sequences, recovered_toks)
|
| 366 |
+
|
| 367 |
+
return recovered_sequences, torch.LongTensor(torch.stack(recovered_toks))
|
| 368 |
+
|
| 369 |
+
def mismatched_positions(s1, s2):
|
| 370 |
+
# 这个函数假定两个字符串的长度相同。
|
| 371 |
+
"""Return the number of positions where two strings differ."""
|
| 372 |
+
|
| 373 |
+
# The number of mismatches will be the sum of positions where characters are not the same
|
| 374 |
+
return sum(1 for c1, c2 in zip(s1, s2) if c1 != c2)
|
| 375 |
+
|
| 376 |
+
def remove_duplicates_triple(filtered_mut_seqs, filtered_mut_probs, filtered_mut_logits):
|
| 377 |
+
seen = {}
|
| 378 |
+
unique_seqs = []
|
| 379 |
+
unique_probs = []
|
| 380 |
+
unique_logits = []
|
| 381 |
+
|
| 382 |
+
for seq, prob, logit in zip(filtered_mut_seqs, filtered_mut_probs, filtered_mut_logits):
|
| 383 |
+
if seq not in seen:
|
| 384 |
+
unique_seqs.append(seq)
|
| 385 |
+
unique_probs.append(prob)
|
| 386 |
+
unique_logits.append(logit)
|
| 387 |
+
seen[seq] = True
|
| 388 |
+
|
| 389 |
+
return unique_seqs, unique_probs, unique_logits
|
| 390 |
+
|
| 391 |
+
def remove_duplicates_double(filtered_mut_seqs, filtered_mut_probs):
|
| 392 |
+
seen = {}
|
| 393 |
+
unique_seqs = []
|
| 394 |
+
unique_probs = []
|
| 395 |
+
|
| 396 |
+
for seq, prob in zip(filtered_mut_seqs, filtered_mut_probs):
|
| 397 |
+
if seq not in seen:
|
| 398 |
+
unique_seqs.append(seq)
|
| 399 |
+
unique_probs.append(prob)
|
| 400 |
+
seen[seq] = True
|
| 401 |
+
|
| 402 |
+
return unique_seqs, unique_probs
|
| 403 |
+
|
| 404 |
+
def mutated_seq(wt_seq, wt_label):
|
| 405 |
+
wt_seq = '!'+ wt_seq
|
| 406 |
+
wt_tok = torch.LongTensor([[tok_to_idx[token] for token in wt_seq]]).to(device)
|
| 407 |
+
wt_prob, wt_pred, _, wt_logit = model.predict(wt_tok)
|
| 408 |
+
|
| 409 |
+
st.write(f'Wild Type: Length = ', len(wt_seq), '\n', wt_seq)
|
| 410 |
+
st.write(f'Wild Type: Label = {wt_label}, Y_pred = {wt_pred.item()}, Y_prob = {wt_prob.item():.2%}')
|
| 411 |
+
|
| 412 |
+
# st.write(n_mut, mlm_tok_num, n_designs_ep, n_sampling_designs_ep, n_mlm_recovery_sampling, mutate2stronger)
|
| 413 |
+
# pbar = tqdm(total=n_mut)
|
| 414 |
+
mutated_seqs = []
|
| 415 |
+
i = 1
|
| 416 |
+
pbar = st.progress(i, text="mutated number of sequence")
|
| 417 |
+
while i <= n_mut:
|
| 418 |
+
if i == 1: seeds_ep = [wt_seq[1:]]
|
| 419 |
+
seeds_next_ep, seeds_probs_next_ep, seeds_logits_next_ep = [], [], []
|
| 420 |
+
for seed in seeds_ep:
|
| 421 |
+
seed_seq, masked_seed_seq, seed_seq_token, masked_seed_seq_token = batch_mlm_seq([seed] * n_designs_ep, continuous_replace = True) ### mask seed with 1 site to "*"
|
| 422 |
+
|
| 423 |
+
seed_prob, seed_pred, _, seed_logit = model.predict(seed_seq_token[0].unsqueeze_(0).to(device))
|
| 424 |
+
_, _, seed_esm_logit, _ = model.predict(masked_seed_seq_token.to(device))
|
| 425 |
+
mut_seqs, mut_toks = recovered_mlm_multi_tokens(masked_seed_seq, masked_seed_seq_token, seed_esm_logit)
|
| 426 |
+
mut_probs, mut_preds, mut_esm_logits, mut_logits = model.predict(mut_toks.to(device))
|
| 427 |
+
|
| 428 |
+
### Filter mut_seqs that mut_prob < seed_prob and mut_prob < wild_prob
|
| 429 |
+
filtered_mut_seqs = []
|
| 430 |
+
filtered_mut_probs = []
|
| 431 |
+
filtered_mut_logits = []
|
| 432 |
+
if mut_by_prob:
|
| 433 |
+
for z in range(len(mut_seqs)):
|
| 434 |
+
if mutate2stronger:
|
| 435 |
+
if mut_probs[z] >= seed_prob and mut_probs[z] >= wt_prob:
|
| 436 |
+
filtered_mut_seqs.append(mut_seqs[z])
|
| 437 |
+
filtered_mut_probs.append(mut_probs[z].cpu().detach().numpy())
|
| 438 |
+
filtered_mut_logits.append(mut_logits[z].cpu().detach().numpy())
|
| 439 |
+
else:
|
| 440 |
+
if mut_probs[z] < seed_prob and mut_probs[z] < wt_prob:
|
| 441 |
+
filtered_mut_seqs.append(mut_seqs[z])
|
| 442 |
+
filtered_mut_probs.append(mut_probs[z].cpu().detach().numpy())
|
| 443 |
+
filtered_mut_logits.append(mut_logits[z].cpu().detach().numpy())
|
| 444 |
+
else:
|
| 445 |
+
for z in range(len(mut_seqs)):
|
| 446 |
+
if mutate2stronger:
|
| 447 |
+
if mut_logits[z] >= seed_logit and mut_logits[z] >= wt_logit:
|
| 448 |
+
filtered_mut_seqs.append(mut_seqs[z])
|
| 449 |
+
filtered_mut_probs.append(mut_probs[z].cpu().detach().numpy())
|
| 450 |
+
filtered_mut_logits.append(mut_logits[z].cpu().detach().numpy())
|
| 451 |
+
else:
|
| 452 |
+
if mut_logits[z] < seed_logit and mut_logits[z] < wt_logit:
|
| 453 |
+
filtered_mut_seqs.append(mut_seqs[z])
|
| 454 |
+
filtered_mut_probs.append(mut_probs[z].cpu().detach().numpy())
|
| 455 |
+
filtered_mut_logits.append(mut_logits[z].cpu().detach().numpy())
|
| 456 |
+
|
| 457 |
+
|
| 458 |
+
|
| 459 |
+
### Save
|
| 460 |
+
seeds_next_ep.extend(filtered_mut_seqs)
|
| 461 |
+
seeds_probs_next_ep.extend(filtered_mut_probs)
|
| 462 |
+
seeds_logits_next_ep.extend(filtered_mut_logits)
|
| 463 |
+
seeds_next_ep, seeds_probs_next_ep, seeds_logits_next_ep = remove_duplicates_triple(seeds_next_ep, seeds_probs_next_ep, seeds_logits_next_ep)
|
| 464 |
+
|
| 465 |
+
### Sampling based on prob
|
| 466 |
+
if len(seeds_next_ep) > n_sampling_designs_ep:
|
| 467 |
+
seeds_probs_next_ep_norm = seeds_probs_next_ep / sum(seeds_probs_next_ep) # Normalize the probabilities
|
| 468 |
+
seeds_index_next_ep = np.random.choice(len(seeds_next_ep), n_sampling_designs_ep, p = seeds_probs_next_ep_norm, replace = False)
|
| 469 |
+
|
| 470 |
+
seeds_next_ep = np.array(seeds_next_ep)[seeds_index_next_ep]
|
| 471 |
+
seeds_probs_next_ep = np.array(seeds_probs_next_ep)[seeds_index_next_ep]
|
| 472 |
+
seeds_logits_next_ep = np.array(seeds_logits_next_ep)[seeds_index_next_ep]
|
| 473 |
+
seeds_mutated_num_next_ep = [mismatched_positions(wt_seq[1:], s) for s in seeds_next_ep]
|
| 474 |
+
|
| 475 |
+
mutated_seqs.extend(list(zip(seeds_next_ep, seeds_logits_next_ep, seeds_probs_next_ep, seeds_mutated_num_next_ep)))
|
| 476 |
+
|
| 477 |
+
seeds_ep = seeds_next_ep
|
| 478 |
+
i += 1
|
| 479 |
+
# pbar.update(1)
|
| 480 |
+
pbar.progress(i/n_mut, text="Mutating")
|
| 481 |
+
# pbar.close()
|
| 482 |
+
st.success('Done', icon="✅")
|
| 483 |
+
mutated_seqs.extend([(wt_seq[1:], wt_logit.item(), wt_prob.item(), 0)])
|
| 484 |
+
mutated_seqs = sorted(mutated_seqs, key=lambda x: x[2], reverse=True)
|
| 485 |
+
mutated_seqs = pd.DataFrame(mutated_seqs, columns = ['mutated_seq', 'predicted_logit', 'predicted_probability', 'mutated_num']).drop_duplicates('mutated_seq')
|
| 486 |
+
return mutated_seqs
|
| 487 |
+
|
| 488 |
+
def read_raw(raw_input):
|
| 489 |
+
ids = []
|
| 490 |
+
sequences = []
|
| 491 |
+
|
| 492 |
+
file = StringIO(raw_input)
|
| 493 |
+
for record in SeqIO.parse(file, "fasta"):
|
| 494 |
+
|
| 495 |
+
# 检查序列是否只包含A, G, C, T
|
| 496 |
+
sequence = str(record.seq.back_transcribe()).upper()
|
| 497 |
+
if not set(sequence).issubset(set("AGCT")):
|
| 498 |
+
st.write(f"Record '{record.description}' was skipped for containing invalid characters. Only A, G, C, T(U) are allowed.")
|
| 499 |
+
continue
|
| 500 |
+
|
| 501 |
+
# 将符合条件的序列添加到列表中
|
| 502 |
+
ids.append(record.id)
|
| 503 |
+
sequences.append(sequence)
|
| 504 |
+
|
| 505 |
+
return ids, sequences
|
| 506 |
+
|
| 507 |
+
def predict_raw(raw_input):
|
| 508 |
+
state_dict = torch.load('model.pt', map_location=torch.device(device))
|
| 509 |
+
new_state_dict = OrderedDict()
|
| 510 |
+
|
| 511 |
+
for k, v in state_dict.items():
|
| 512 |
+
name = k.replace('module.','')
|
| 513 |
+
new_state_dict[name] = v
|
| 514 |
+
|
| 515 |
+
model = CNN_linear().to(device)
|
| 516 |
+
model.load_state_dict(new_state_dict, strict = False)
|
| 517 |
+
model.eval()
|
| 518 |
+
st.write(model)
|
| 519 |
+
# st.write('====Parse Input====')
|
| 520 |
+
ids, seqs = read_raw(raw_input)
|
| 521 |
+
|
| 522 |
+
# st.write('====Predict====')
|
| 523 |
+
res_pd = pd.DataFrame()
|
| 524 |
+
for wt_seq, wt_id in zip(seqs, ids):
|
| 525 |
+
try:
|
| 526 |
+
st.write(wt_id, wt_seq)
|
| 527 |
+
res = mutated_seq(wt_seq, wt_id)
|
| 528 |
+
st.write(res)
|
| 529 |
+
res_pd.append(res)
|
| 530 |
+
except:
|
| 531 |
+
st.write('====Please Try Again this sequence: ', wt_id, wt_seq)
|
| 532 |
+
# st.write(pred)
|
| 533 |
+
return res_pd
|
| 534 |
+
|
| 535 |
+
# Run
|
| 536 |
+
if st.button("Predict and Mutate"):
|
| 537 |
+
if uploaded:
|
| 538 |
+
result = predict_raw(uploaded.getvalue().decode())
|
| 539 |
+
else:
|
| 540 |
+
result = predict_raw(seq)
|
| 541 |
+
|
| 542 |
+
result_file = result.to_csv(index=False)
|
| 543 |
+
st.download_button("Download", result_file, file_name=output_filename+".csv")
|
| 544 |
+
st.dataframe(result)
|
| 545 |
+
|
| 546 |
+
|
| 547 |
+
|
| 548 |
+
|
| 549 |
+
|
| 550 |
+
|
git.sh
ADDED
|
@@ -0,0 +1,3 @@
|
|
|
|
|
|
|
|
|
|
|
|
|
| 1 |
+
git add .
|
| 2 |
+
git commit -m "app"
|
| 3 |
+
git push
|