本文尝试利用机器学习模型中对特征值的筛选方法，来筛选单细胞中不同亚型中关键的gene。
1. 数据获取：
为方便我们从单细胞分析seurat对象过渡到python分析环境，以下提供seurat 对象转换成H5AD的代码，请参考使用。
# 数据文件格式转换
library(SeuratDisk)
library(Seurat)
SaveH5Seurat(seurat.obj,filename="test.h5seurat", overwrite = TRUE)
Convert("test.h5seurat", dest = "h5ad", overwrite = TRUE)
## 数据读取和切分
df = anndata.read("test.h5ad")
X_train, X_test, y_train, y_test = train_test_split(df.X, df.obs.celltype4, random_state=1)
2. 利用随机森林模型选择基因的排序重要性
2.1. 理论背景：随机排列每个特征的值，然后监控模型性能下降的程度。如果获得了更大的下降意味着特征更重要
from sklearn.ensemble import RandomForestClassifier
from sklearn.inspection import permutation_importance
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import anndata
rf = RandomForestClassifier(n_estimators=100, random_state=1)
rf.fit(X_train, y_train)
baseline = rf.score(X_test, y_test)
result = permutation_importance(rf, X_test, y_test, n_repeats=10, random_state=1, scoring='accuracy')
importances = result.importances_mean
### 同时，rf模型可以单独输出重要性结果，更加方便
rff = RandomForestClassifier(n_estimators=100, random_state=1)
rff.fit(df.X, df.obs.celltype4)
importances2 = rff.feature_importances_
### 数据保存
import pandas as pd
dff = pd.DataFrame({"a":list(df.var.index),"b" : list(importances)})
dff = dff.sort_values(by="b",ascending=False)
dff.to_csv("frist_importance.csv")
### 可视化
plt.bar(range(len(importances)), importances)
plt.xlabel('Feature Index')
plt.ylabel('Permutation Importance')
plt.show()
rf1.png
2.2. 理论基础：迭代的删除每一个特征，查看模型的准确性下降情况，以推断其特征解释性。
from sklearn.metrics import accuracy_score
import numpy as np
rf3 = RandomForestClassifier(n_estimators=100, random_state=1)
rf3.fit(X_train, y_train)
base_acc = accuracy_score(y_test, rf3.predict(X_test))
importances3 = []
for i in range(X_train.shape[1]):
X_temp = np.delete(X_train, i, axis=1)
rf3.fit(X_temp, y_train)
acc = accuracy_score(y_test, rf3.predict(np.delete(X_test, i, axis=1)))
importances3.append(base_acc - acc)
plt.bar(range(len(importances3)), importances3)
plt.show()
3. 递归去除不重要的特征，以获取到需要的特征数量。
from sklearn.feature_selection import RFE
rf = RandomForestClassifier()
rfe = RFE(rf, n_features_to_select=10)
rfe.fit(df.X, df.obs.celltype4)
print(rfe.ranking_)
4. XGBoost特性重要性
理论基础：计算一个特性用于跨所有树拆分数据的次数。更多的分裂意味着更重要。
import xgboost as xgb
model = xgb.XGBClassifier()
model.fit(df.X, df.obs.celltype4)
importances = model.feature_importances_
importances = pd.Series(importances, index=range(df.X.shape[1]))
importances.plot.bar()
5. 主成分分析 PCA
理论基础：对特征进行主成分分析，并查看每个主成分的解释方差比。在前几个组件上具有较高负载的特性更为重要。注意，该方法不能获取到对应的gene name信息。
from sklearn.decomposition import PCA
pca = PCA()
pca.fit(df.X)
plt.bar(range(pca.n_components_), pca.explained_variance_ratio_)
plt.xlabel('PCA components')
plt.ylabel('Explained Variance')
6. 方差分析 ANOVA
使用f_classif()获得每个特征的方差分析f值。f值越高，表明特征与目标的相关性越强。
from sklearn.feature_selection import f_classif
fval = f_classif(df.X, df.obs.celltype4)
fval = pd.Series(fval[0], index=range(df.X.shape[1]))
fval.plot.bar()
7. 卡方检验
使用chi2()获得每个特征的卡方统计信息。得分越高的特征越有可能独立于目标。
from sklearn.feature_selection import chi2
chi_scores = chi2(df.X, df.obs.celltype4)
chi_scores = pd.Series(chi_scores[0], index=range(df.X.shape[1]))
chi_scores.plot.bar()
