[PYTHON] Apprentissage automatique à l'aide de données d'expression génique
Apprentissage automatique à l'aide de données d'expression génique
(2017.02.02 OSX El Capitan)
À partir de profils d'expression génique, en utilisant des données d'expression génique provenant de divers tissus humains, Apprenons et déduisons de quel tissu humain provient chaque échantillon. J'ai essayé de l'exécuter en me référant à Python Machine Learning Programming.
Préparation des données d'expression génique pour l'apprentissage automatique
Obtenez des données du portail GTEx
GTEx Portal publie des données d'expression génique provenant de divers tissus humains. Dans la dernière version de l'analyse GTEx V6p, de 53 tissus humains et 544 donneurs Des données d'expression génique pour un total de 8555 échantillons peuvent être obtenues.
En enregistrant un compte sur le site, vous pourrez télécharger des données d'expression génique.
Il y a un lien de téléchargement sur la page qui peut être atteint depuis Datasets-> Download sur la barre de navigation du site.
RPKM (lectures par kilobase de transcription par million de lectures mappées) Corrigé Le fichier de données d'expression génique porte les noms de fichier suivants.
- GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct.gz
Il est fourni au format GCT (Gene Cluster Text file format). Les colonnes sont des échantillons (les noms de colonnes sont des noms d'échantillons) et les lignes sont des quantités de caractéristiques (les noms de lignes sont des noms de gènes).
L'annotation (informations d'étiquette) de chaque échantillon est Il est fourni dans la version GTEx Analysis V6 et porte les noms de fichier suivants.
- GTEx_Data_V6_Annotations_SampleAttributesDS.txt
L'ID d'échantillon se trouve dans la première colonne et le tissu à partir duquel est dérivé dans les sixième et septième colonnes.
GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct.gz et GTEx_Data_V6_Annotations_SampleAttributesDS.txt Cliquez sur le lien de téléchargement de chaque fichier à télécharger. Ensuite, double-cliquez sur le fichier .gz pour le décompresser.
Vérifier avec la commande Unix / Linux
Utilisez les commandes Linux pour examiner les 10 premières lignes et 3 colonnes du fichier de données d'expression génique. À partir des deux premières lignes, vous pouvez voir que le niveau d'expression (niveau de fonctionnalité) du gène 56238 est de 8555 échantillons. De plus, à partir du nom du gène dans la première colonne, vous pouvez voir que Ensembl Gene ID est utilisé. Les données de niveau d'expression elles-mêmes commencent à partir de la 3ème colonne de la 4ème ligne.
cut -f1,2,3 GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct | head
#1.2
56238 8555
Name Description GTEX-111CU-1826-SM-5GZYN
ENSG00000223972.4 DDX11L1 0
ENSG00000227232.4 WASH7P 6.50895977020264
ENSG00000243485.2 MIR1302-11 0
ENSG00000237613.2 FAM138A 0
ENSG00000268020.2 OR4G4P 0
ENSG00000240361.1 OR4G11P 0
ENSG00000186092.4 OR4F5 0
Lecture des données d'expression génique
Utilisez le [package pandas read_csv] de python (http://pandas.pydata.org/pandas-docs/stable/generated/pandas.read_csv.html) pour lire les données.
Puisqu'il est délimité par des tabulations, c'est sep = '\ t', et je veux sauter les deux premières lignes, donc je le règle sur skiprows = 2.
La première ligne lue en premier est l'en-tête (nom de la colonne).
Parce que je veux lire 8555 exemples de données (jusqu'à 8557 colonnes) à l'exception de la description dans la deuxième colonne.
ʻUsecols = [0] + list (range (2,8557)) . Ensuite, définissez ʻindex_clo = 0 et spécifiez la première colonne comme nom de ligne.
import pandas as pd
usecols = [0] + list(range(2,8557))
df1 = pd.read_csv('GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct', sep='\t', skiprows=2, usecols=usecols, index_col=0)
>>> df1.index
Index(['ENSG00000223972.4', 'ENSG00000227232.4', 'ENSG00000243485.2',
'ENSG00000237613.2', 'ENSG00000268020.2', 'ENSG00000240361.1',
'ENSG00000186092.4', 'ENSG00000238009.2', 'ENSG00000233750.3',
'ENSG00000237683.5',
...
'ENSG00000198886.2', 'ENSG00000210176.1', 'ENSG00000210184.1',
'ENSG00000210191.1', 'ENSG00000198786.2', 'ENSG00000198695.2',
'ENSG00000210194.1', 'ENSG00000198727.2', 'ENSG00000210195.2',
'ENSG00000210196.2'],
dtype='object', name='Name', length=56238)
>>> df1.columns
Index(['GTEX-111CU-1826-SM-5GZYN', 'GTEX-111FC-0226-SM-5N9B8',
'GTEX-111VG-2326-SM-5N9BK', 'GTEX-111YS-2426-SM-5GZZQ',
'GTEX-1122O-2026-SM-5NQ91', 'GTEX-1128S-2126-SM-5H12U',
'GTEX-113IC-0226-SM-5HL5C', 'GTEX-117YX-2226-SM-5EGJJ',
'GTEX-11DXW-0326-SM-5H11W', 'GTEX-11DXX-2326-SM-5Q5A2',
...
'GTEX-ZVE2-0006-SM-51MRW', 'GTEX-ZVP2-0005-SM-51MRK',
'GTEX-ZVT2-0005-SM-57WBW', 'GTEX-ZVT3-0006-SM-51MT9',
'GTEX-ZVT4-0006-SM-57WB8', 'GTEX-ZVTK-0006-SM-57WBK',
'GTEX-ZVZP-0006-SM-51MSW', 'GTEX-ZVZQ-0006-SM-51MR8',
'GTEX-ZXES-0005-SM-57WCB', 'GTEX-ZXG5-0005-SM-57WCN'],
dtype='object', length=8555)
transformation log, tracé de la distribution de densité
Ajoutez le pseudo-comptage 1 au niveau d'expression génique lu (RPKM) pour la conversion log2.
import numpy as np
df1pc = df1 + 1
df1log = np.log2(df1pc)
Sélectionnez au hasard 10 échantillons et tracez la distribution de densité des niveaux d'expression génique. Utilisez pandas DataFrame.plot.density pour le traçage.
import matplotlib.pyplot as plt
random_cols = np.random.choice(8555, 10)
df1log.iloc[:, random_cols].plot.density(fontsize=10)
plt.savefig('gtex_log_rpkm_density_random10sample.png', dpi=150)
plt.show()
Étant donné que la distribution varie considérablement d'un échantillon à l'autre, recherchez une distribution de référence à partir de l'ensemble de données d'apprentissage. Il peut être préférable d'effectuer une correction quantile, etc.
Chargement des informations sur l'organisation de l'échantillon (étiquette de classe)
À partir de GTEx_Data_V6_Annotations_SampleAttributesDS.txt, qui contient des exemples d'informations sur l'organisation, Lisez l'ID de l'échantillon dans la 1ère colonne et les informations sur l'organisation dans la 7ème colonne.
usecols = [0, 6]
df2 = pd.read_csv('GTEx_Data_V6_Annotations_SampleAttributesDS.txt', sep='\t', usecols=usecols, index_col=0)
En transposant les lignes et colonnes des données d'expression génique et en les combinant avec la base de données des informations sur les tissus de l'échantillon Donne des informations sur l'organisation. Utilisez pandas concat.
df = pd.concat([df1log.T, df2], axis=1, join_axes=[df1log.T.index])
Depuis Python Machine Learning Programming 4.2.2 (P99) Comme il est difficile de gérer les informations d'organisation telles quelles, convertissez la chaîne de caractères en entier et remplacez-la. Puisque les données d'expression génique ont été transposées et combinées, les lignes sont les échantillons et les colonnes les quantités de caractéristiques et les informations sur les tissus.
from sklearn.preprocessing import LabelEncoder
class_le = LabelEncoder()
label = class_le.fit_transform(df['SMTSD'].values)
label_r = class_le.inverse_transform(label)
df['SMTSD'] = label
Divisé en données d'entraînement et données de test
Depuis Python Machine Learning Programming 4.3 (P102) Divisez les données d'expression génique marquées en quantité de caractéristiques (X) et étiquette de classe (y).
smtsd_index = len(df.columns) - 1
X, y = df.iloc[:, :smtsd_index].values, df.iloc[:, smtsd_index].values
Ensuite, il est divisé en données d'entraînement et données de test.
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = \
train_test_split(X, y, test_size=0.3, random_state=0)
Apprentissage et raisonnement
Essayons d'apprendre et de déduire par la méthode de voisinage k. Utilisez sklearn.neighbors.KNeighborsClassifier.
Il est temps de supprimer les variables inutiles pour faire de la mémoire.
del df
del df1
del df1pc
del df1log
del df2
k = 3, 4 Exécuter en parallèle. Il a fallu environ 10 minutes pour calculer sur le Intel Core i7 4GHz.
from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=3, n_jobs=4)
knn.fit(X_train, y_train)
print('Training accuracy: ', knn.score(X_train, y_train))
print('Test accuracy: ', knn.score(X_test, y_test))
Comme indiqué ci-dessous, le taux de réponse correcte dans les données de test était d'environ 91%.
>>> print('Training accuracy: ', knn.score(X_train, y_train))
Training accuracy: 0.95641282565130259
>>> print('Test accuracy: ', knn.score(X_test, y_test))
Test accuracy: 0.90962212699649392
Sélection de la quantité de caractéristiques en fonction du niveau d'expression génique
Sélectionnez le montant de la fonction. Les gènes dont l'expression génique totale ne dépasse pas 500 seront supprimés. Programmation Python Machine Learning 4.5.2 (P113) Sélection de fonctionnalités séquentielles et Créez une classe pour sélectionner des fonctionnalités en vous référant au Code source de scikit-learn.
from sklearn.base import TransformerMixin
import numpy as np
class GexSumFS(TransformerMixin):
"""
Feature selection based on the sum of gene expression values
"""
def __init__(self, sum):
self.sum = sum
def fit(self, X, y=None):
self.indices_ = np.where(np.sum(X, axis=0) > self.sum)[0]
return self
def transform(self, X, y=None):
return X[:, self.indices_]
Utilisez la classe ci-dessus pour filtrer par niveau d'expression total.
sumfs = GexSumFS(sum = 500)
X_train_sumfs = sumfs.fit_transform(X_train)
X_test_sumfs = sumfs.transform(X_test)
Le nombre de gènes (nombre de caractéristiques) a diminué à 27754.
>>> X_train_sumfs.shape[1]
27754
Fonctionnalités Apprenons et déduisons avec les données sélectionnées.
knn.fit(X_train_sumfs, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs, y_test))
Le taux de réponse correcte était d'environ 91%, ce qui était presque le même qu'avant la sélection des fonctionnalités. Le taux de précision a légèrement diminué, La sélection des fonctionnalités avait l'avantage de réduire un peu le temps de calcul.
>>> print('Training accuracy: ', knn.score(X_train_sumfs, y_train))
Training accuracy: 0.955577822311
>>> print('Test accuracy: ', knn.score(X_test_sumfs, y_test))
Test accuracy: 0.908843007402
Lorsque le niveau d'expression génique est standardisé (moyenne 0, dispersion 1), les gènes à faible niveau d'expression ont tendance à être bruyants. La sélection des fonctionnalités en fonction du niveau d'expression devrait être plus efficace.
Quantile normalization
À partir du graphique de distribution de densité, la distribution du niveau d'expression était assez différente d'un échantillon à l'autre. Essayez la correction quantile, puis apprenez et inférez. Créez une classe qui effectue une correction quantile comme indiqué ci-dessous.
class QuantileNormalization(TransformerMixin):
"""
Quantile normalization
"""
def __init__(self):
pass
def fit(self, X, y=None):
self.ref_dist_ = np.zeros(X.shape[1])
X_rank = pd.DataFrame(X).rank(axis=1, method='min').astype(int)
for i in range(X.shape[1]):
indices = np.where(X_rank == i+1)
if len(indices[0]) > 0:
self.ref_dist_[i] = np.mean(X[indices])
else:
self.ref_dist_[i] = self.ref_dist_[i-1]
return self
def transform(self, X, y=None):
X_norm = pd.DataFrame(X).rank(axis=1, method='min').values
for i in range(X.shape[1]):
X_norm[X_norm == i+1] = self.ref_dist_[i]
return X_norm
Utilisez la classe ci-dessus pour effectuer une correction quantile. Il a fallu une journée entière pour calculer, probablement parce que la méthode de mise en œuvre était mauvaise.
q_norm = QuantileNormalization()
X_train_sumfs_qnorm = q_norm.fit_transform(X_train_sumfs)
X_test_sumfs_qnorm = q_norm.transform(X_test_sumfs)
Regardez la distribution de densité du profil d'expression génique corrigé et confirmez que la distribution de densité est uniforme.
random_rows = np.random.choice(range(X_train_sumfs_qnorm.shape[0]), 10)
pd.DataFrame(X_train_sumfs_qnorm[random_rows, :]).T.plot.density(fontsize=10)
plt.savefig('gtex_qnorm_density_random10sample.png', dpi=150)
plt.show()
Apprenons et déduisons avec des données corrigées quantiles
knn.fit(X_train_sumfs_qnorm, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm, y_test))
Le taux de réponse correcte de l'inférence a augmenté d'environ 0,6%.
>>> print('Training accuracy: ', knn.score(X_train_sumfs_qnorm, y_train))
Training accuracy: 0.956078824315
>>> print('Test accuracy: ', knn.score(X_test_sumfs_qnorm, y_test))
Test accuracy: 0.915855083755
Standardisation
Essayons d'apprendre et d'inférer après avoir standardisé les données avec une distribution uniforme par correction quantile. Programmation Python Machine Learning 4.4 Reportez-vous à Alignement des échelles de fonctionnalités (P106).
from sklearn.preprocessing import StandardScaler
stdsc = StandardScaler()
X_train_sumfs_qnorm_std = stdsc.fit_transform(X_train_sumfs_qnorm)
X_test_sumfs_qnorm_std = stdsc.transform(X_test_sumfs_qnorm)
knn.fit(X_train_sumfs_qnorm_std, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_std, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_std, y_test))
>>> print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_std, y_train))
Training accuracy: 0.951736806947
>>> print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_std, y_test))
Test accuracy: 0.906505648617
Dans ce cas, la normalisation semble réduire le taux de précision. La standardisation de gènes avec de très petites variances peut entraîner de grandes valeurs. Cela peut être la raison de la diminution du taux de précision.
Normalisation
Au lieu de normaliser, essayez de normaliser (mise à l'échelle de 0 à 1) avant d'apprendre et de déduire. Programmation Python Machine Learning 4.4 Reportez-vous à Alignement des échelles de fonctionnalités (P105).
from sklearn.preprocessing import MinMaxScaler
mms = MinMaxScaler()
X_train_sumfs_qnorm_norm = mms.fit_transform(X_train_sumfs_qnorm)
X_test_sumfs_qnorm_norm = mms.transform(X_test_sumfs_qnorm)
knn.fit(X_train_sumfs_qnorm_norm, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_norm, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_norm, y_test))
>>> print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_norm, y_train))
Training accuracy: 0.952571810287
>>> print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_norm, y_test))
Test accuracy: 0.910011686794
Dans ce cas, la normalisation semble également réduire le taux de précision. Cela peut être dû au fait que la normalisation augmente l'influence du gène avec moins de variation.
Sélection des caractéristiques en fonction de la taille de la dispersion
L'élimination préliminaire des gènes moins dispersés peut aider à la standardisation. Créez une classe qui élimine les entités avec de petites variations, comme indiqué ci-dessous.
from sklearn.base import TransformerMixin
import numpy as np
class GexVarFS(TransformerMixin):
"""
Feature selection based on the variance of gene expression values
"""
def __init__(self, var):
self.var = var
def fit(self, X, y=None):
self.indices_ = np.where(np.var(X, axis=0) > self.var)[0]
return self
def transform(self, X, y=None):
return X[:, self.indices_]
En utilisant la classe ci-dessus, essayez d'apprendre et de déduire après la sélection de la fonctionnalité (distribuer 1 ou moins est supprimé).
varfs = GexVarFS(var = 1)
X_train_sumfs_qnorm_varfs = varfs.fit_transform(X_train_sumfs_qnorm)
X_test_sumfs_qnorm_varfs = varfs.transform(X_test_sumfs_qnorm)
knn.fit(X_train_sumfs_qnorm_varfs, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_varfs, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_varfs, y_test))
Bien que le nombre de gènes ait diminué et que le temps de calcul ait considérablement diminué, Le taux de précision a légèrement diminué en raison de la sélection des fonctionnalités comme indiqué ci-dessous.
>>> print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_varfs, y_train))
Training accuracy: 0.957915831663
>>> print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_varfs, y_test))
Test accuracy: 0.913128165173
Ensuite, standardisez les données sélectionnées par l'entité en fonction de la taille de la distribution, puis apprenez et inférez.
X_train_sumfs_qnorm_varfs_std = stdsc.fit_transform(X_train_sumfs_qnorm_varfs)
X_test_sumfs_qnorm_varfs_std = stdsc.transform(X_test_sumfs_qnorm_varfs)
knn.fit(X_train_sumfs_qnorm_varfs_std, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_varfs_std, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_varfs_std, y_test))
Dans ce cas, la normalisation a augmenté le taux de précision. De plus, le taux de précision a légèrement augmenté en raison de la sélection et de la normalisation des fonctionnalités par rapport à la sélection des fonctionnalités précédentes.
>>> print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_varfs_std, y_train))
Training accuracy: 0.958750835003
>>> print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_varfs_std, y_test))
Test accuracy: 0.918971562135
Résumé
En fin de compte, le taux de réponse correcte de l'inférence est d'environ 92%, par rapport au cas où la valeur RPKM log-convertie est utilisée telle quelle. Il s'est amélioré d'environ 1%. Réduction de dimension par analyse en composantes principales, De meilleurs résultats peuvent être obtenus en changeant la méthode d'apprentissage automatique en quelque chose comme SVM.
Enfin, vérifions le taux de réponse correct pour chaque organisation.
y_pred = knn.predict(X_test_sumfs_qnorm_varfs_std)
inv_dict = dict(zip(y, class_le.inverse_transform(y)))
for i in np.unique(label):
accuracy_score = np.sum((y_pred[y_test == i] == i).astype(int)) / len(y_pred[y_test == i])
print(inv_dict[i] + ': ' + str(accuracy_score))
Adipose - Subcutaneous: 1.0
Adipose - Visceral (Omentum): 0.984126984127
Adrenal Gland: 1.0
Artery - Aorta: 0.985294117647
Artery - Coronary: 0.911764705882
Artery - Tibial: 0.988636363636
Bladder: 1.0
Brain - Amygdala: 0.95
Brain - Anterior cingulate cortex (BA24): 0.76
Brain - Caudate (basal ganglia): 0.709677419355
Brain - Cerebellar Hemisphere: 0.838709677419
Brain - Cerebellum: 0.883720930233
Brain - Cortex: 0.864864864865
Brain - Frontal Cortex (BA9): 0.53125
Brain - Hippocampus: 0.8
Brain - Hypothalamus: 0.964285714286
Brain - Nucleus accumbens (basal ganglia): 0.705882352941
Brain - Putamen (basal ganglia): 0.344827586207
Brain - Spinal cord (cervical c-1): 0.928571428571
Brain - Substantia nigra: 0.666666666667
Breast - Mammary Tissue: 0.823529411765
Cells - EBV-transformed lymphocytes: 1.0
Cells - Transformed fibroblasts: 1.0
Cervix - Ectocervix: 1.0
Cervix - Endocervix: 0.0
Colon - Sigmoid: 0.725
Colon - Transverse: 0.796875
Esophagus - Gastroesophageal Junction: 0.509803921569
Esophagus - Mucosa: 1.0
Esophagus - Muscularis: 0.876543209877
Fallopian Tube: 0.5
Heart - Atrial Appendage: 0.98275862069
Heart - Left Ventricle: 1.0
Kidney - Cortex: 1.0
Liver: 1.0
Lung: 1.0
Minor Salivary Gland: 1.0
Muscle - Skeletal: 1.0
Nerve - Tibial: 1.0
Ovary: 1.0
Pancreas: 1.0
Pituitary: 1.0
Prostate: 1.0
Skin - Not Sun Exposed (Suprapubic): 0.807692307692
Skin - Sun Exposed (Lower leg): 0.925619834711
Small Intestine - Terminal Ileum: 0.96
Spleen: 1.0
Stomach: 0.903225806452
Testis: 1.0
Thyroid: 1.0
Uterus: 0.95
Vagina: 0.909090909091
Whole Blood: 1.0
Résumé du code
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import LabelEncoder
from sklearn.cross_validation import train_test_split
from sklearn.base import TransformerMixin
from sklearn.preprocessing import StandardScaler
from sklearn.neighbors import KNeighborsClassifier
### read gene expression matrix
usecols = [0] + list(range(2,8557))
df1 = pd.read_csv('GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct', sep='\t', skiprows=2, usecols=usecols, index_col=0)
df1pc = df1 + 1
df1log = np.log2(df1pc)
### plot gene expression profiles
random_cols = np.random.choice(8555, 10)
df1log.iloc[:, random_cols].plot.density(fontsize=10)
plt.savefig('gtex_log_rpkm_density_random10sample.png', dpi=150)
#plt.show()
### read class labels
usecols = [0, 6]
df2 = pd.read_csv('GTEx_Data_V6_Annotations_SampleAttributesDS.txt', sep='\t', usecols=usecols, index_col=0)
### attach class labels
df = pd.concat([df1log.T, df2], axis=1, join_axes=[df1log.T.index])
### convert string labels
class_le = LabelEncoder()
label = class_le.fit_transform(df['SMTSD'].values)
label_r = class_le.inverse_transform(label)
df['SMTSD'] = label
### split into training and test data
smtsd_index = len(df.columns) - 1
X, y = df.iloc[:, :smtsd_index].values, df.iloc[:, smtsd_index].values
X_train, X_test, y_train, y_test = \
train_test_split(X, y, test_size=0.3, random_state=0)
### save memory usage
del df
del df1
del df1pc
del df1log
del df2
### feature selection based on gene expression values
class GexSumFS(TransformerMixin):
"""
Feature selection based on the sum of gene expression values
"""
def __init__(self, sum):
self.sum = sum
def fit(self, X, y=None):
self.indices_ = np.where(np.sum(X, axis=0) > self.sum)[0]
return self
def transform(self, X, y=None):
return X[:, self.indices_]
sumfs = GexSumFS(sum = 500)
X_train_sumfs = sumfs.fit_transform(X_train)
X_test_sumfs = sumfs.transform(X_test)
### quantile normalization
class QuantileNormalization(TransformerMixin):
"""
Quantile normalization
"""
def __init__(self):
pass
def fit(self, X, y=None):
self.ref_dist_ = np.zeros(X.shape[1])
X_rank = pd.DataFrame(X).rank(axis=1, method='min').astype(int)
for i in range(X.shape[1]):
indices = np.where(X_rank == i+1)
if len(indices[0]) > 0:
self.ref_dist_[i] = np.mean(X[indices])
else:
self.ref_dist_[i] = self.ref_dist_[i-1]
return self
def transform(self, X, y=None):
X_norm = pd.DataFrame(X).rank(axis=1, method='min').values
for i in range(X.shape[1]):
X_norm[X_norm == i+1] = self.ref_dist_[i]
return X_norm
q_norm = QuantileNormalization()
X_train_sumfs_qnorm = q_norm.fit_transform(X_train_sumfs)
X_test_sumfs_qnorm = q_norm.transform(X_test_sumfs)
random_rows = np.random.choice(range(X_train_sumfs_qnorm.shape[0]), 10)
pd.DataFrame(X_train_sumfs_qnorm[random_rows, :]).T.plot.density(fontsize=10)
plt.savefig('gtex_qnorm_density_random10sample.png', dpi=150)
#plt.show()
### feature selection based on gene expression variance
from sklearn.base import TransformerMixin
import numpy as np
class GexVarFS(TransformerMixin):
"""
Feature selection based on the variance of gene expression values
"""
def __init__(self, var):
self.var = var
def fit(self, X, y=None):
self.indices_ = np.where(np.var(X, axis=0) > self.var)[0]
return self
def transform(self, X, y=None):
return X[:, self.indices_]
varfs = GexVarFS(var = 1)
X_train_sumfs_qnorm_varfs = varfs.fit_transform(X_train_sumfs_qnorm)
X_test_sumfs_qnorm_varfs = varfs.transform(X_test_sumfs_qnorm)
### standardization
stdsc = StandardScaler()
X_train_sumfs_qnorm_varfs_std = stdsc.fit_transform(X_train_sumfs_qnorm_varfs)
X_test_sumfs_qnorm_varfs_std = stdsc.transform(X_test_sumfs_qnorm_varfs)
### learning and prediction
knn = KNeighborsClassifier(n_neighbors=3, n_jobs=4)
knn.fit(X_train_sumfs_qnorm_varfs_std, y_train)
print('Training accuracy: ', knn.score(X_train_sumfs_qnorm_varfs_std, y_train))
print('Test accuracy: ', knn.score(X_test_sumfs_qnorm_varfs_std, y_test))
### accuracy score for each class
y_pred = knn.predict(X_test_sumfs_qnorm_varfs_std)
inv_dict = dict(zip(y, class_le.inverse_transform(y)))
for i in np.unique(label):
accuracy_score = np.sum((y_pred[y_test == i] == i).astype(int)) / len(y_pred[y_test == i])
print(inv_dict[i] + ': ' + str(accuracy_score))
Recommended Posts