使用逻辑回归碱基特异性贡献

在基因编辑过程中，有很多统计学方法用来探索 gRNA 碱基特异性贡献。主要目的是看到底是哪些位置的哪些碱基影响了基因编辑的效率。今天主要介绍的方法是使用逻辑回归来计算碱基特异性贡献。

逻辑回归主要原理

数据准备

数据准备步骤很简单，只需要准备两列的表格，第一列是 gRNA 序列，第二列是编辑效率(最好是 Count,我们在 Python 中计算效率，这样有利于浮点数的准确性)。

依赖的第三方包

pandas
numpy
sklearn
scipy
logomaker
matplotlib

分析流程详述

数据预处理

首先，导入必须的包

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import Ridge
from sklearn.preprocessing import OneHotEncoder
from scipy.stats import pearsonr,spearmanr
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Ridge
from sklearn.linear_model import Lasso
from sklearn.model_selection import cross_val_score
from sklearn import metrics
import pickle
import logomaker
import matplotlib.pyplot as plt

然后，读取数据，并预处理数据。需要注意的是，gRNA 的编辑效率可能为 0，当我们做logit 转换的时候，会出现 Inf 的情况, 我们这里也进行一定的转化。

## 读取数据
df=pd.read_csv("train_data.tsv",sep="\t")

### this is for indel
ratio = df["indel"] / df["total"]
# 解决 logit 的 0/1 问题
epsilon = 1e-6
ratio = ratio.clip(epsilon, 1 - epsilon)

# logit transform
df["logit_ratio"] = np.log(ratio / (1 - ratio))

特征提取

我们使用 onehot 编码来提取特征,

### one-hot
spacers = df["spacer"].apply(list).tolist()
L = len(spacers[0])  # spacer 长度

# 为每个位置设置 A C G T
encoder = OneHotEncoder(categories=[['A','C','G','T']] * L,
                        sparse_output=False)

X = encoder.fit_transform(spacers)
y = df["logit_ratio"].values

划分训练集和测试集

#######################
# 2. 划分训练 & 测试集
#######################
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)

模型训练

我们这里使用 5 折交叉验证来评估模型的性能

##########################
# 线性回归模型
##########################
lr = LinearRegression()

scores_lr = cross_val_score(
    lr,
    X_train,
    y_train,
    scoring="neg_mean_squared_error",
    cv=5
)

print("LinearRegression 5-fold MSE (each fold):", -scores_lr)
print("Mean MSE:", -scores_lr.mean())

模型评估

一旦交叉验证的结果比较稳定，我们就可以训练模型，并在测试集上进行评估

#######################
# 4. 测试集预测相关性
#######################

# 用最佳模型（lasso 更稳定）
lr.fit(X_train, y_train)
test_pred = lr.predict(X_test)

corr, pval = pearsonr(y_test, test_pred)
print("Test Pearson correlation:", corr)
print("P-value:", pval)

模型保存

#######################
# 5. 模型保存
#######################
pickle.dump(lr, open("lasso_model.pkl", "wb"))

参数的可视化

我们可以使用 Sequence Logo 来可视化碱基特异性贡献,首先提取权重.

### 提取权重

L = 20
coef = lr.coef_                 # shape = (L*4,)
W = coef.reshape(L, 4)             # shape = (L,4)

# 导出权重 DataFrame
weight_df = pd.DataFrame(
    W,
    columns=["A", "C", "G", "T"]
)

weight_df["position"] = np.arange(1, L+1)

weight_df.head()

然后使用 logomaker 来绘制 logo

#######################
# 6. 绘制 sequence logo
#######################

color_scheme = {
    'A': "#CA5E4E",
    'C': "#E0BC5E",
    'G': "#99B85E",
    'T': "#5F9BD1",
}

plt.figure(figsize=(12, 3))

logo = logomaker.Logo(
    weight_df[["A","C","G","T"]],
    color_scheme=color_scheme,
    flip_below=False
)

# 只保留 X 轴; 去除 Y 轴与外框
logo.ax.spines["top"].set_visible(False)
logo.ax.spines["right"].set_visible(False)
logo.ax.spines["left"].set_visible(False)
logo.ax.spines["bottom"].set_visible(True)  # 只保留 bottom

logo.ax.set_xticks(range(0, L))
logo.ax.set_xticklabels(np.arange(1, L+1))

logo.ax.set_ylabel("Weight (logit-space)")
plt.title("Sequence Logo of Model Weights")

# Get the current figure to add a subtitle
fig = plt.gcf()
#fig.text(0.5, 0.85, "Logistic Regression (r = 0.58)", ha='center', va='top', fontsize=10)

plt.savefig("Indel_motif.pdf")
plt.show()

Indel_motif