<jupyter_start><jupyter_text>First pass with no filtering<jupyter_code>#Cys Shuffle fitting with filtering
INFILE = "./UstD_CysShuffDataset.csv" #Which dataset file to load
FILTER = 0.0 #Disregard measurements below this filter value

data = [[],[]]
numlost = 0
#Load data and generate RoF curve
knownacts = {}
with open(INFILE,'r') as f:
    line = f.readline()
    while len(line) > 0:
        temp = line.split(',')
        if len(temp[0]) > 0:
            if (FILTER < float(temp[1])):
                data[0].append(temp[0])
                data[1].append(float(temp[1]))
                knownacts[data[0][-1]] = data[1][-1]
            else:
                numlost += 1
        line = f.readline()

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
ax.scatter(list(range(len(data[0]))), sorted(data[1],reverse=True))
ax.set_title("Retention of Function Curve (%s removed)"%numlost)
ax.set_xlabel("Rank")
ax.set_ylabel("Fold Activity")
fig.show()

#Get number of unique sequences
print("Total # measurements: %s"%len(data[0]))
print("Unique sequences: %s"%len(set(data[0])))
#Encode sequences
data[1] = np.array(data[1])
log_measurements = np.log(data[1])
encoder = SequenceEncoder(data[0])
encseqs, mapping = encoder.OneHotOrthoEncoder(givemap=True)
if len(encseqs) > len(log_measurements):
    print("WARN: Less measurements than vector size!  Watch for overfitting")
#Perform LOOCV to get statistics
predictions = []
for i in range(len(encseqs)):
    reg = sklm.LinearRegression(fit_intercept=False).fit(encseqs[:i] + encseqs[i+1:], np.delete(log_measurements, i))
    predictions.append(reg.predict([encseqs[i]]))
predictions = np.exp(predictions).flatten()

mae = np.average(np.abs(data[1] - predictions))#mean absolute error
rsq = r2_score(data[1], predictions)

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
x = data[1]
y = predictions
min_p = np.min([np.min(x), np.min(y)])
max_p = np.max([np.max(x), np.max(y)])

ax.plot([min_p, max_p],[min_p, max_p],color='orange')
ax.scatter(x, y)
ax.set_title("Leave-One-Out Cross-Validation (MAE: %.2f, R^2: %.2f)"%(mae, rsq))
ax.set_xlabel("Measured Fold-Activity")
ax.set_ylabel("Predicted Fold-Activity")
fig.show()
#Retrain on full dataset
reg = sklm.LinearRegression(fit_intercept=False).fit(encseqs, log_measurements)
predictions = reg.predict(encseqs)
predictions = np.exp(predictions)

mae = np.average(np.abs(data[1] - predictions))#mean absolute error
rsq = r2_score(data[1], predictions)

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
x = data[1]
y = predictions
min_p = np.min([np.min(x), np.min(y)])
max_p = np.max([np.max(x), np.max(y)])

ax.plot([min_p, max_p],[min_p, max_p],color='orange')
ax.scatter(x, y)
ax.set_title("Full Dataset Model (MAE: %.2f, R^2: %.2f)"%(mae, rsq))
ax.set_xlabel("Measured Fold-Activity")
ax.set_ylabel("Predicted Fold-Activity")
fig.show()<jupyter_output>/home/jon/anaconda3/envs/pytorch/lib/python3.7/site-packages/ipykernel_launcher.py:21: UserWarning: Matplotlib is currently using module://ipykernel.pylab.backend_inline, which is a non-GUI backend, so cannot show the figure.<jupyter_text>Second pass with reduced dataset<jupyter_code>#Cys Shuffle fitting with filtering
INFILE = "./UstD_CysShuffDataset.csv" #Which dataset file to load
FILTER = 0.1 #Disregard measurements below this filter value

data = [[],[]]
numlost = 0
#Load data and generate RoF curve
knownacts = {}
with open(INFILE,'r') as f:
    line = f.readline()
    while len(line) > 0:
        temp = line.split(',')
        if len(temp[0]) > 0:
            if (FILTER < float(temp[1])):
                data[0].append(temp[0])
                data[1].append(float(temp[1]))
                knownacts[data[0][-1]] = data[1][-1]
            else:
                numlost += 1
        line = f.readline()

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
ax.scatter(list(range(len(data[0]))), sorted(data[1],reverse=True))
ax.set_title("Retention of Function Curve (%s removed)"%numlost)
ax.set_xlabel("Rank")
ax.set_ylabel("Fold Activity")
fig.show()

#Get number of unique sequences
print("Total # measurements: %s"%len(data[0]))
print("Unique sequences: %s"%len(set(data[0])))
#Encode sequences
data[1] = np.array(data[1])
log_measurements = np.log(data[1])
encoder = SequenceEncoder(data[0])
encseqs, mapping = encoder.OneHotOrthoEncoder(givemap=True)
if len(encseqs) > len(log_measurements):
    print("WARN: Less measurements than vector size!  Watch for overfitting")
#Perform LOOCV to get statistics
predictions = []
for i in range(len(encseqs)):
    reg = sklm.LinearRegression(fit_intercept=False).fit(encseqs[:i] + encseqs[i+1:], np.delete(log_measurements, i))
    predictions.append(reg.predict([encseqs[i]]))
predictions = np.exp(predictions).flatten()

mae = np.average(np.abs(data[1] - predictions))#mean absolute error
rsq = r2_score(data[1], predictions)

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
x = data[1]
y = predictions
min_p = np.min([np.min(x), np.min(y)])
max_p = np.max([np.max(x), np.max(y)])

ax.plot([min_p, max_p],[min_p, max_p],color='orange')
ax.scatter(x, y)
ax.set_title("Leave-One-Out Cross-Validation (MAE: %.2f, R^2: %.2f)"%(mae, rsq))
ax.set_xlabel("Experimental Fold-Activity")
ax.set_ylabel("Predicted Fold-Activity")
fig.show()
#Retrain on full dataset
reg = sklm.LinearRegression(fit_intercept=False).fit(encseqs, log_measurements)
predictions = reg.predict(encseqs)
predictions = np.exp(predictions)

mae = np.average(np.abs(data[1] - predictions))#mean absolute error
rsq = r2_score(data[1], predictions)

fig, ax = plt.subplots()
fig.patch.set_facecolor('white')
x = data[1]
y = predictions
min_p = np.min([np.min(x), np.min(y)])
max_p = np.max([np.max(x), np.max(y)])

ax.plot([min_p, max_p],[min_p, max_p],color='orange')
ax.scatter(x, y)
ax.set_title("Full Dataset Model (MAE: %.2f, R^2: %.2f)"%(mae, rsq))
ax.set_xlabel("Measured Fold-Activity")
ax.set_ylabel("Predicted Fold-Activity")
fig.show()
#Reconstruct weight matrix in a readable format
num_pos = len(mapping)
possible_aas = []
pos_list = []
counter = 1
for pos in mapping:
    temp = list(pos)[1:]
    possible_aas += temp
    pos_list += ["Pos%s"%(counter)] * len(temp)
    counter += 1
unique_aas = list(set(possible_aas))
unique_aas.sort()
unique_pos = list(set(pos_list))
unique_pos.sort()
weights = np.full((len(unique_aas), num_pos), '-', dtype=object)
weights = pd.DataFrame(weights, index=unique_aas, columns=unique_pos)
model_weights = np.exp(reg.coef_)
for i in range(len(possible_aas)):
    weights.at[possible_aas[i], pos_list[i]] = "%.1f"%model_weights[i+1]
weights.columns = ['C122X', 'C139X', 'C227X', 'C236X', 'C428X']
print("Parent Bias: %.1f"%(model_weights[0]))
display(weights)
#Generate all permutations of library
permutationseqs = [[x] for x in list(mapping[0].keys())]
for i in range(1, len(mapping)):
    temp = []
    for j in permutationseqs:
        for k in mapping[i]:
            temp.append(j + [k])
    permutationseqs = temp[:]
permutationencs = {}
for seq in permutationseqs:
    enc = [1]
    for pos in range(len(seq)):
        enc += mapping[pos][seq[pos]]
    permutationencs["".join(seq)] = enc

#Create library prediction and report sorted activities
permutationpreds = reg.predict([permutationencs[x] for x in permutationencs])
permutationpreds = np.exp(permutationpreds)
permutationpreds, permutationseqs = zip(*sorted(zip(permutationpreds, permutationseqs), reverse=True))
print("TOP PREDICTED ACTIVITIES")
for i in range(10):
    key = "".join(permutationseqs[i])
    if key in knownacts:
        print("%s: %.2f(%.2f)"%(key, permutationpreds[i],knownacts[key]))
    else:
        print("%s: %.2f"%(key, permutationpreds[i]))<jupyter_output>TOP PREDICTED ACTIVITIES
ACCSC: 3.11
ASCSC: 3.07
ACASC: 2.90
ASASC: 2.87
SCCSC: 2.86
SSCSC: 2.82
SCASC: 2.67(2.71)
SSASC: 2.63
GCCSC: 1.92(2.60)
GSCSC: 1.89
{'fraction_non_alphanumeric': 0.16935381490415133, 'fraction_numerical': 0.018915830900088866, 'mean_word_length': 4.8967065868263475, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '331321', 'license': 'MIT', 'n_tokens_mistral': 3070, 'n_tokens_neox': 2957, 'n_words': 617, 'path': 'UstDCysShuffLRM.ipynb', 'repo_name': 'bullerlaboratory/NatCatalUstDPaper', 'size': '7877'}