TMB背景知识介绍
免疫疗法是近年来肿瘤治疗的热点，肿瘤突变负荷(Tumor Mutational Burden，TMB)可以作为免疫疗法的biomarker，TMB高的患者更有可能在免疫治疗中获益。
全外显子测序(WES)是评价TMB的金标准，但是价格还较高，限制了其在临床上的应用。目前已有研究和商品基于panel内基因，发现其和WES水平的TMB具有高度的一致性，可以用来评价临床的TMB，但是这些报道均没有提及如何从人类近30000个基因中筛选出几百个基因组成panel。我们的目的就在于建立一种有效的基因筛选方法，使得到的panel能够很好的评价TMB。
TMB计算
从TCGA数据库下载了555例肺癌患者的WES突变数据，整理后数据结构如下：
| Sample | Gene1 | Gene2 | ... | Genep | total_gene_num |
|---|---|---|---|---|---|
| sample1 | 3 | 2 | ... | 0 | 113 |
| sample2 | 2 | 0 | ... | 4 | 96 |
| ... | ... | ... | ... | ... | ... |
| samplen | 1 | 5 | ... | 9 | 1412 |
统计每一个样本中每一个基因的突变数目，它在WES水平上的突变数目就是所有基因突变数目求和，即total_gene_num。
因此我们的任务就是：从p个基因中挑选出一定数目的基因，通过这些基因能够很好地预测total_gene_num，现有的研究一般是选择两者的相关系数作为衡量的指标。比如：选择了Gene1、Gene2和Genen，对这3个基因的突变数目求和，然后计算r2_score([5, 6, 15], [113, 96, 1412])，如果值高，这说明选择的基因好。
机器学习
对于一组data，可以看作是以下形式：
其中是响应变量，
是自变量，
是随机误差，
是真实的关系
而机器学习能做的是：
其中是
的预测，
是假设的模型
我们要做的就是找到一个合适的，来进行学习。
Lasso
一般来说，是复杂的，未知的，不管是为了求解方便，还是我们的能力有限，
都和真实的
有差距。“所有的模型都是错的，但有的模型是有用的”，我们不一定能真正地找到一个正确的
，但是我们能够做出较好地预测，这就够了。
幸运的是，在我们要研究的TMB panel基因筛选的问题中，真实的是已知的！！！
在上面对TMB计算的描述中：
真实的模型就是一个截距项(intercept)为0，所有系数(coefficient)为1的线性模型，由于特征数p >> 样本数n，所以使用最小二乘法是无法对方程组求解的。
我们的目的是对这个模型做特征选择，尽可能多地去掉特征，仅保留几百个基因。由此自然就想到了正则化的
模型。通过调节
，它可以尽可能地将更多的coef压缩为0，达到特征筛选的目的。
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn import linear_model
from sklearn.metrics import r2_score, mean_squared_error
from sklearn.model_selection import train_test_split, StratifiedShuffleSplit
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import GridSearchCV
读入数据
file = "E:/patent/mutation.num.stat.xls"
df = pd.read_csv(file, sep="\t")
分层划分数据集
new_df = df.copy()
new_df["cat"] = np.ceil(new_df["MutNum_exome"] / 1.5)
new_df["cat"].where(new_df["cat"] < 5, 5.0, inplace=True)
split = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=1412)
for train_index, test_index in split.split(new_df, new_df["cat"]):
strat_train_set = new_df.iloc[train_index]
strat_test_set = new_df.iloc[test_index]
将数据集转换为ndarray
X_train = strat_train_set[df.columns[3:]]
y_train = strat_train_set["MutNum_exome"]
X_test = strat_test_set[df.columns[3:]]
y_test = strat_test_set["MutNum_exome"]
网格搜索确定参数
param_grid = [
{
"fit_intercept": [False, True],
"positive": [True, False],
"selection": ["random"],
"alpha": [0.1, 0.2, 0.3, 0.4, 0.5],
"max_iter": [3000],
}
]
model = linear_model.Lasso()
grid_search = GridSearchCV(
model,
param_grid,
cv=10,
scoring="neg_mean_squared_error"
)
grid_search.fit(X_train, y_train)
F:\Program_Files\Anaconda\lib\site-packages\sklearn\model_selection\_search.py:841: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
DeprecationWarning)
GridSearchCV(cv=10, error_score='raise-deprecating',
estimator=Lasso(alpha=1.0, copy_X=True, fit_intercept=True, max_iter=1000,
normalize=False, positive=False, precompute=False, random_state=None,
selection='cyclic', tol=0.0001, warm_start=False),
fit_params=None, iid='warn', n_jobs=None,
param_grid=[{'fit_intercept': [False, True], 'positive': [True, False], 'selection': ['random'], 'alpha': [0.1, 0.2, 0.3, 0.4, 0.5], 'max_iter': [3000]}],
pre_dispatch='2*n_jobs', refit=True, return_train_score='warn',
scoring='neg_mean_squared_error', verbose=0)
最优参数
grid_search.best_params_
{'alpha': 0.1,
'fit_intercept': False,
'max_iter': 3000,
'positive': True,
'selection': 'random'}
从最优参数可以看出，模型是符合我们预期的：不需要拟合截距，所有的系数都为正值。
拟合结果
y_test_pred = grid_search.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
best_alpha = grid_search.best_estimator_.alpha
nonzero_coef = sum(grid_search.best_estimator_.coef_ != 0)
print(
f"最优alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"非零特征数: {nonzero_coef}"
)
最优alpha: 0.1
RMSE: 60.47737241913124
非零特征数: 273
new_df["MutNum_exome"].mean()
121.3963963963964
可以看到，整个数据集中total_gene_num的均值为121.39，而模型的RMSE达到60.47，几乎是均值的一半了，这是不太能接受的结果。
尝试调整alpha获得更优的结果，我们注意到目前最优alpha为0.1，是我们候选的alpha中的端点值，所以继续向下筛选。
LassoCV确定最优alpha
cadidate_alphas = np.linspace(0.01, 0.1, 1000)
model = linear_model.LassoCV(
alphas=cadidate_alphas,
fit_intercept=False,
positive=True,
selection="random",
max_iter=3000,
cv=10
)
model.fit(X_train, y_train)
LassoCV(alphas=array([0.01 , 0.01009, ..., 0.09991, 0.1 ]), copy_X=True,
cv=10, eps=0.001, fit_intercept=False, max_iter=3000, n_alphas=100,
n_jobs=None, normalize=False, positive=True, precompute='auto',
random_state=None, selection='random', tol=0.0001, verbose=False)
y_test_pred = model.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
best_alpha = model.alpha_
r2 = r2_score(y_test, y_test_pred)
nonzero_coef = sum(model.coef_ != 0)
print(
f"最优alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"R2: {r2}\n"
f"非零特征数: {nonzero_coef}"
)
最优alpha: 0.01
RMSE: 57.024710487728775
R2: 0.9119547827475557
非零特征数: 405
RMSE有所改善，但是选出的alpha依然是端点值，尝试继续缩小alpha
继续缩小alpha，减对模型的限制
cadidate_alphas = np.linspace(0.001, 0.01, 1000)
# max_iter should increase
model2 = linear_model.LassoCV(
alphas=cadidate_alphas,
fit_intercept=False,
positive=True,
selection="random",
max_iter=5000,
cv=10
)
model2.fit(X_train, y_train)
LassoCV(alphas=array([0.001 , 0.00101, ..., 0.00999, 0.01 ]), copy_X=True,
cv=10, eps=0.001, fit_intercept=False, max_iter=5000, n_alphas=100,
n_jobs=None, normalize=False, positive=True, precompute='auto',
random_state=None, selection='random', tol=0.0001, verbose=False)
y_test_pred = model2.predict(X_test)
rmse = np.sqrt(mean_squared_error(y_test, y_test_pred))
best_alpha = model2.alpha_
r2 = r2_score(y_test, y_test_pred)
nonzero_coef = sum(model2.coef_ != 0)
print(
f"最优alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"R2: {r2}\n"
f"非零特征数: {nonzero_coef}"
)
最优alpha: 0.001
RMSE: 56.6309058203357
R2: 0.913166639680416
非零特征数: 523
越小的alpha越接近普通的线性模型，但是我们注意到在alpha=0.01时，非零的特征数为405，我们的训练集中样本数为444，已经非常接近了，如果再减小alpha的话会出现n < p的情况，普通的线性模型已经不适合了。
所以在当前样本量的情况下，alpha=0.01已经基本是极限情况了，57.02的RMSE也基本是最优的情况，想要改善的话需要增加样本量。
如果我们持续放松对模型的限制，也就失去了特征筛选的意义。
基因数量与TMB的关系
从直观的角度就能理解，基因数目越大，panel和WES的就越接近1。
下图展示了随机选取一定数目基因，二者相关性的情况：
可以看出随着基因数目增加，相关性也稳定增大，我们似乎不需要做特征筛选，只需要保证一定数目的基因就好了。
这篇文章讨论了基因数目对的影响，微信公众号_从基因Panel测序数据评估15种肿瘤的TMB分布及与基因间的关系也做了描述了该种现象：做TMB评价要选择大panel，但他们没有阐述背后的原因。
这种情况可以从“总体均值和样本均值”的角度来解释：
TMB的计算方法是“体细胞突变数目/对应区域的长度”，即每Mb碱基的体细胞突变数，是一个均值。对于WES水平的TMB就是总体均值，而panel水平的TMB就是样本均值。我们知道，样本均值是总体均值的无偏估计，随着样本数量增加，样本均值的也越接近总体均值(标准误差减小)。
结论
如果以机器学习的观点来看TMB panel基因筛选，它就是一个真实模型为线性模型的学习问题，且真实模型没有intercept，coef均为正值。
由于p >> n，所以如果要做特征选择，可以选用Lasso模型。
受限于样本量，alpha只能松弛到一定程度：使得非零系数样本数，得到RMSE是有一个最小值的
因此，基因筛选主要取决于两点：
- 样本量
样本量越大，可以筛选出的基因数越多 - 基因数目
基因数肯定是越大越准确，但也就失去了筛选的意义，所以可以在RMSE无明显增大的情况下，增大alpha，获得更小的基因数。
基因数目的影响远大于具体选哪些基因的影响
在该问题的研究中，应该转换思路：不是从所有基因中选出一定数量的基因组成panel，这个panel可以评价TMB；而是，我已经有一个panel，再来评价该panel计算得到TMB是否和WES的一致
