TMB背景知識(shí)介紹
免疫療法是近年來(lái)腫瘤治療的熱點(diǎn)百框，腫瘤突變負(fù)荷(Tumor Mutational Burden待讳，TMB)可以作為免疫療法的biomarker魁莉，TMB高的患者更有可能在免疫治療中獲益。
全外顯子測(cè)序(WES)是評(píng)價(jià)TMB的金標(biāo)準(zhǔn)馅而，但是價(jià)格還較高逾冬，限制了其在臨床上的應(yīng)用黍聂。目前已有研究和商品基于panel內(nèi)基因，發(fā)現(xiàn)其和WES水平的TMB具有高度的一致性身腻，可以用來(lái)評(píng)價(jià)臨床的TMB产还，但是這些報(bào)道均沒(méi)有提及如何從人類近30000個(gè)基因中篩選出幾百個(gè)基因組成panel。我們的目的就在于建立一種有效的基因篩選方法嘀趟，使得到的panel能夠很好的評(píng)價(jià)TMB雕沉。
TMB計(jì)算
從TCGA數(shù)據(jù)庫(kù)下載了555例肺癌患者的WES突變數(shù)據(jù)，整理后數(shù)據(jù)結(jié)構(gòu)如下：
| Sample | Gene1 | Gene2 | ... | Genep | total_gene_num |
|---|---|---|---|---|---|
| sample1 | 3 | 2 | ... | 0 | 113 |
| sample2 | 2 | 0 | ... | 4 | 96 |
| ... | ... | ... | ... | ... | ... |
| samplen | 1 | 5 | ... | 9 | 1412 |
統(tǒng)計(jì)每一個(gè)樣本中每一個(gè)基因的突變數(shù)目去件，它在WES水平上的突變數(shù)目就是所有基因突變數(shù)目求和坡椒，即total_gene_num扰路。
因此我們的任務(wù)就是：從p個(gè)基因中挑選出一定數(shù)目的基因，通過(guò)這些基因能夠很好地預(yù)測(cè)total_gene_num倔叼，現(xiàn)有的研究一般是選擇兩者的相關(guān)系數(shù)作為衡量的指標(biāo)汗唱。比如：選擇了Gene1、Gene2和Genen丈攒，對(duì)這3個(gè)基因的突變數(shù)目求和哩罪，然后計(jì)算r2_score([5, 6, 15], [113, 96, 1412])，如果值高巡验，這說(shuō)明選擇的基因好际插。
機(jī)器學(xué)習(xí)
對(duì)于一組data，可以看作是以下形式：
其中是響應(yīng)變量显设，
是自變量框弛，
是隨機(jī)誤差，
是真實(shí)的關(guān)系
而機(jī)器學(xué)習(xí)能做的是：
其中是
的預(yù)測(cè)捕捂，
是假設(shè)的模型
我們要做的就是找到一個(gè)合適的瑟枫，來(lái)進(jìn)行學(xué)習(xí)。
Lasso
一般來(lái)說(shuō)指攒，是復(fù)雜的慷妙，未知的，不管是為了求解方便允悦，還是我們的能力有限膝擂，
都和真實(shí)的
有差距∠冻冢“所有的模型都是錯(cuò)的架馋，但有的模型是有用的”，我們不一定能真正地找到一個(gè)正確的
驶鹉，但是我們能夠做出較好地預(yù)測(cè)绩蜻，這就夠了铣墨。
幸運(yùn)的是室埋，在我們要研究的TMB panel基因篩選的問(wèn)題中，真實(shí)的是已知的Ｒ猎肌Ｒο！
在上面對(duì)TMB計(jì)算的描述中：
真實(shí)的模型就是一個(gè)截距項(xiàng)(intercept)為0屡律，所有系數(shù)(coefficient)為1的線性模型腌逢，由于特征數(shù)p >> 樣本數(shù)n，所以使用最小二乘法是無(wú)法對(duì)方程組求解的超埋。
我們的目的是對(duì)這個(gè)模型做特征選擇搏讶，盡可能多地去掉特征佳鳖，僅保留幾百個(gè)基因。由此自然就想到了正則化的
模型媒惕。通過(guò)調(diào)節(jié)
系吩，它可以盡可能地將更多的coef壓縮為0，達(dá)到特征篩選的目的妒蔚。
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
讀入數(shù)據(jù)
file = "E:/patent/mutation.num.stat.xls"
df = pd.read_csv(file, sep="\t")
分層劃分?jǐn)?shù)據(jù)集
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]
將數(shù)據(jù)集轉(zhuǎn)換為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"]
網(wǎng)格搜索確定參數(shù)
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)
最優(yōu)參數(shù)
grid_search.best_params_
{'alpha': 0.1,
'fit_intercept': False,
'max_iter': 3000,
'positive': True,
'selection': 'random'}
從最優(yōu)參數(shù)可以看出穿挨，模型是符合我們預(yù)期的：不需要擬合截距，所有的系數(shù)都為正值肴盏。
擬合結(jié)果
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"最優(yōu)alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"非零特征數(shù): {nonzero_coef}"
)
最優(yōu)alpha: 0.1
RMSE: 60.47737241913124
非零特征數(shù): 273
new_df["MutNum_exome"].mean()
121.3963963963964
可以看到科盛，整個(gè)數(shù)據(jù)集中total_gene_num的均值為121.39，而模型的RMSE達(dá)到60.47菜皂，幾乎是均值的一半了贞绵，這是不太能接受的結(jié)果。
嘗試調(diào)整alpha獲得更優(yōu)的結(jié)果幌墓，我們注意到目前最優(yōu)alpha為0.1但壮，是我們候選的alpha中的端點(diǎn)值，所以繼續(xù)向下篩選常侣。
LassoCV確定最優(yōu)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"最優(yōu)alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"R2: {r2}\n"
f"非零特征數(shù): {nonzero_coef}"
)
最優(yōu)alpha: 0.01
RMSE: 57.024710487728775
R2: 0.9119547827475557
非零特征數(shù): 405
RMSE有所改善蜡饵，但是選出的alpha依然是端點(diǎn)值，嘗試?yán)^續(xù)縮小alpha
繼續(xù)縮小alpha胳施，減對(duì)模型的限制
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"最優(yōu)alpha: {best_alpha}\n"
f"RMSE: {rmse}\n"
f"R2: {r2}\n"
f"非零特征數(shù): {nonzero_coef}"
)
最優(yōu)alpha: 0.001
RMSE: 56.6309058203357
R2: 0.913166639680416
非零特征數(shù): 523
越小的alpha越接近普通的線性模型溯祸，但是我們注意到在alpha=0.01時(shí)，非零的特征數(shù)為405舞肆，我們的訓(xùn)練集中樣本數(shù)為444焦辅，已經(jīng)非常接近了，如果再減小alpha的話會(huì)出現(xiàn)n < p的情況椿胯，普通的線性模型已經(jīng)不適合了筷登。
所以在當(dāng)前樣本量的情況下，alpha=0.01已經(jīng)基本是極限情況了哩盲，57.02的RMSE也基本是最優(yōu)的情況前方，想要改善的話需要增加樣本量。
如果我們持續(xù)放松對(duì)模型的限制廉油，也就失去了特征篩選的意義惠险。
基因數(shù)量與TMB的關(guān)系
從直觀的角度就能理解，基因數(shù)目越大抒线，panel和WES的就越接近1班巩。
下圖展示了隨機(jī)選取一定數(shù)目基因，二者相關(guān)性的情況：
可以看出隨著基因數(shù)目增加嘶炭，相關(guān)性也穩(wěn)定增大抱慌，我們似乎不需要做特征篩選逊桦，只需要保證一定數(shù)目的基因就好了。
這篇文章討論了基因數(shù)目對(duì)的影響抑进，微信公眾號(hào)_從基因Panel測(cè)序數(shù)據(jù)評(píng)估15種腫瘤的TMB分布及與基因間的關(guān)系也做了描述了該種現(xiàn)象：做TMB評(píng)價(jià)要選擇大panel卫袒，但他們沒(méi)有闡述背后的原因。
這種情況可以從“總體均值和樣本均值”的角度來(lái)解釋：
TMB的計(jì)算方法是“體細(xì)胞突變數(shù)目/對(duì)應(yīng)區(qū)域的長(zhǎng)度”单匣，即每Mb堿基的體細(xì)胞突變數(shù)夕凝，是一個(gè)均值。對(duì)于WES水平的TMB就是總體均值户秤，而panel水平的TMB就是樣本均值码秉。我們知道，樣本均值是總體均值的無(wú)偏估計(jì)鸡号，隨著樣本數(shù)量增加转砖，樣本均值的也越接近總體均值(標(biāo)準(zhǔn)誤差減小)。
結(jié)論
如果以機(jī)器學(xué)習(xí)的觀點(diǎn)來(lái)看TMB panel基因篩選鲸伴，它就是一個(gè)真實(shí)模型為線性模型的學(xué)習(xí)問(wèn)題府蔗，且真實(shí)模型沒(méi)有intercept，coef均為正值汞窗。
由于p >> n姓赤，所以如果要做特征選擇，可以選用Lasso模型仲吏。
受限于樣本量不铆，alpha只能松弛到一定程度：使得非零系數(shù)樣本數(shù)，得到RMSE是有一個(gè)最小值的
因此裹唆，基因篩選主要取決于兩點(diǎn)：
- 樣本量
樣本量越大誓斥，可以篩選出的基因數(shù)越多 - 基因數(shù)目
基因數(shù)肯定是越大越準(zhǔn)確，但也就失去了篩選的意義许帐，所以可以在RMSE無(wú)明顯增大的情況下劳坑，增大alpha，獲得更小的基因數(shù)成畦。
基因數(shù)目的影響遠(yuǎn)大于具體選哪些基因的影響
在該問(wèn)題的研究中距芬，應(yīng)該轉(zhuǎn)換思路：不是從所有基因中選出一定數(shù)量的基因組成panel，這個(gè)panel可以評(píng)價(jià)TMB羡鸥；而是蔑穴，我已經(jīng)有一個(gè)panel忠寻，再來(lái)評(píng)價(jià)該panel計(jì)算得到TMB是否和WES的一致
