J'ai essayé l'analyse de données IRMf avec python (Introduction au décodage des informations cérébrales)
Qiita n'a pas d'exemple d'analyse des données IRMf ... !!
L'autre jour, j'ai pensé: "Il existe divers articles sur qiita ... Je pense qu'il n'y a pas beaucoup d'exemples d'analyse de ** données IRMf ** en utilisant python ou l'apprentissage automatique. Par conséquent, j'ai cherché dans qiita un exemple d'analyse des données IRMf de l'IRMf.
cette? Pas vraiment ...? !! Il existe des articles de commentaires dans les domaines des outils et des neurosciences, mais je n'ai trouvé aucune analyse utilisant des données d'IRMf réelles.
Donc ** cette fois, nous analyserons l'apprentissage automatique à l'aide de données IRMf . Même s'il s'agit d'une image du cerveau, c'est une image en trois dimensions après tout. Une vue en coupe d'une image bidimensionnelle obtenue en découpant le cerveau en deux est formée par chevauchement, et chaque pixel a des informations dans une unité appelée « boxel **». Je suis juste là. En ce qui concerne les images cérébrales 4D, c'est juste que chaque information en forme de boîte en trois dimensions contient des informations de série chronologique (lors de la numérisation avec l'IRMf).
Les images du cerveau sont souvent publiées dans un format spécial appelé format NIfTY, qui doit être lu. Cette fois, au lieu de Matlab, qui serait souvent utilisé par les étudiants et les chercheurs, nous utiliserons la bibliothèque ** Nilearn ** de python pour lire et analyser les images du cerveau (gardez à l'esprit que tout le monde peut les analyser). Par conséquent, j'ai arrêté de payer Matlab). L'analyse d'apprentissage automatique elle-même est apprise et prédite par le SVM ** scikit-learn ** que toute personne ayant effectué un apprentissage automatique l'aura utilisée ou en aura entendu parler.
De plus, bien qu'il s'agisse d'une approche utilisant l'apprentissage automatique, cette fois nous utiliserons des données d'IRMf qui masquent une région cérébrale spécifique pour effectuer une analyse qui prédit "ce qui était vu par l'IRMf à ce moment-là à partir de l'activité cérébrale" * *. Dans l'industrie de l'analyse cérébrale, cela semble s'appeler ** décodage des informations cérébrales **, mais il ne s'agit que d'apprentissage automatique (désolé pour les gens de l'industrie). Dans le passé, nous avons utilisé une approche telle que le codage, dans lequel le comportement effectué dans une tâche expérimentale était utilisé comme variable explicative, et la région du cerveau qui répondait de manière statistiquement significative lorsqu'elle était effectuée par IRMf était détectée comme une variable indépendante. Je pense que cela venait du fait que cela s'appelait décodage parce qu'il adoptait l'approche opposée (c'est-à-dire au lieu de regarder la région du cerveau qui était statistiquement significativement excitée par la tâche, mais l'activité d'une région du cerveau. Prédire la tâche à partir de).
À propos des données réellement utilisées
Cette fois, je l'ai analysé en me référant au tutoriel d'analyse de décodage de Nilearn: Machine learning for Neuro-Imaging in Python, qui est utilisé pour la lecture de données IRMf, etc. ** L'expérience Haxby 2001 ** utilisée ici Utilisez l'image du cerveau publiée dans.
A introduction tutorial to fMRI decoding
Si vous suivez le didacticiel ci-dessus, il indique que vous pouvez facilement collecter un ensemble de données de sujets à partir de l'expérience Haxby 2001 avec le code suivant.
from nilearn import datasets
#Par défaut, les données du deuxième sujet sont récupérées
haxby_dataset = datasets.fetch_haxby()
Cependant, il semble qu'il y ait un problème avec la base de données NIRTC qui fournit les données, mais je n'ai pas pu la télécharger, j'ai donc pris le moyen de la télécharger directement et de la sauvegarder localement. Vous pouvez télécharger subj2-2010.01.14.tar.gz contenant les données du deuxième sujet à partir de cette URL, donc si vous souhaitez analyser les données, veuillez le télécharger à partir d'ici (même les données d'autres numéros de sujet) n'a pas d'importance).
http://data.pymvpa.org/datasets/haxby2001/
Dans cette expérience, les sujets sont présentés avec diverses catégories de stimuli visuels (tels que des ciseaux, des visages humains et des chats) dans un scanner IRMf. Dans ce tutoriel, nous allons prédire quelle catégorie un sujet regarde avec l'IRMf à partir de l'activité cérébrale IRMf enregistrée dans la zone du tractus visuel cortical ventral (qui est impliqué dans la représentation de la couleur et de la forme).
Procédure d'analyse
- Dans le code source, le nom du dossier téléchargé a été changé en subj2.
Charger et visualiser des images
Tout d'abord, lisons les données IRMf 4D. Comme mentionné ci-dessus, il s'agit d'un format qui contient des informations de série chronologique lors de l'exécution d'une tâche expérimentale dans une image cérébrale qui est des données boxel (données tridimensionnelles).
Cependant, puisqu'il s'agit de données IRMf 4D, elles ne peuvent pas être visualisées par un traitement correspondant aux données IRMf 3D (sans informations de série temporelle) telles que nilearn.plotting.plot_epi. Par conséquent, les données d'IRMf 4D visualisées sont moyennées et converties en une seule donnée d'IRMf 3D. Ceci peut être réalisé en important depuis nilearn.image import mean_img.
fmri_filename = "~/Desktop/nilearn/subj2/bold.nii.gz"
from nilearn import plotting
from nilearn.image import mean_img
plotting.view_img(mean_img(fmri_filename), threshold=None)
résultat
Tout d'abord, j'ai pu visualiser l'image du cerveau. C'est formidable de pouvoir visualiser en quelques lignes seulement.
Extraire les fonctionnalités des données IRMf et les convertir en données au format numpy
Peut-être que certaines personnes étaient soulagées de voir le mot engourdi. Ici, nous utilisons une fonction appelée nilearn.input_data.NiftiMasker pour masquer les données IRMf du tractus visuel cortical ventral et les convertir au format numpy.
Tout d'abord, visualisons les données IRMf de l'appareil visuel cortical ventral uniquement.
mask_filename = "~/Desktop/nilearn/subj2/mask4_vt.nii.gz"
anat_filename = "~/Desktop/nilearn/subj2/anat.nii.gz"
#Essayez d'afficher la plage du masque sur l'arrière-plan de l'image anatomique du sujet
plotting.plot_roi(mask_filename, bg_img=anat_filename,
cmap='Paired')
résultat
Ensuite, créons un masque. Je souhaite normaliser les données afin d'améliorer la précision de l'apprentissage automatique. Le masque sera extrait sous la forme d'un tableau 2D lorsqu'il sera prêt pour l'apprentissage automatique avec nilearn.
from nilearn.input_data import NiftiMasker
masker = NiftiMasker(mask_img=mask_filename, standardize=True) #Standardiser le masque
fmri_masked = masker.fit_transform(fmri_filename) #Appliquer le signal BOLD au masque normalisé
En passant, en exécutant masker.generate_report (), vous pouvez afficher l'image NIfTY du masque et vérifier divers paramètres et ainsi de suite.
Vous avez maintenant extrait la matrice de données masquée sous la forme fmri_masked au format numpy. Je vais vraiment l'essayer.
# fmri_masked se présente sous la forme d'un tableau numpy.
print(fmri_masked)
#La taille de ces données est le nombre de boxels et de points de séries chronologiques(Nombre de scans)Il se compose du nombre total de fichiers.
# shape[0]Est le nombre de scans, shape[1]Est le nombre de boxels
print(fmri_masked.shape)
Résultat de sortie
[[ 7.6757896e-01 2.3108697e+00 -2.0519443e-01 ... -1.0261143e+00
8.7993503e-02 2.0720518e+00]
[ 5.5640817e-01 1.6833434e+00 -2.4644937e-01 ... -7.0238107e-01
-3.4570047e-01 2.0341001e+00]
[ 7.6757896e-01 1.9186659e+00 1.0802225e-03 ... -9.9374104e-01
-2.7630943e-01 2.1479552e+00]
...
[-4.2905563e-01 -1.6896105e+00 -7.4150854e-01 ... -1.5440876e+00
1.8054217e+00 -1.6709718e-01]
[-1.4749455e-01 -1.8072717e+00 -2.4644937e-01 ... -1.7707009e+00
1.5452052e+00 7.8169477e-01]
[-2.1788482e-01 -1.4542881e+00 1.0802225e-03 ... -1.6412076e+00
1.2676411e+00 8.9554977e-01]]
(1452, 464)
Comme confirmé dans «fmri_masked.shape», 464 stocke le nombre de cellules de boîte dans la région masquée du cerveau, et 1452 stocke des informations de série chronologique (l'IRMf a acquis des données IRMf avec une résolution temporelle de TR). Normalement, il est réglé à des intervalles d'environ 1 à 2,5 secondes, mais le décodage analyse en utilisant des informations sur une région cérébrale spécifique et un neurofeedback qui mesure l'activité cérébrale en temps réel et entraîne le sujet. Dans la méthode appelée, je pense qu'un TR assez court était souvent fixé).
À titre de test, si vous retirez les informations de série chronologique des trois premiers boxels, ce sera comme suit.
#Affichez les trois premiers boxels sélectionnés dans la matrice.
import matplotlib.pyplot as plt
plt.plot(fmri_masked[5:150, :3])
plt.title('Voxel Time Series')
plt.xlabel('Scan number')
plt.ylabel('Normalized signal')
plt.tight_layout()
A propos, les gens avisés peuvent se demander pourquoi la partie fmri_masked [5: 150 ,: 3] n'est pas sous la forme de fmri_masked [: 150 ,: 3]. C'est ce qu'on appelle un ** scan factice **, car le premier scan qui commence à imager l'IRMf a une valeur de signal élevée et n'est souvent pas utilisé dans l'analyse réelle. La valeur du signal est élevée car l'effet T1 n'est pas stable (veuillez consulter Google sur d'autres sites pour plus de détails).
Capturez les étiquettes d'action pendant les tâches expérimentales
Maintenant que nous avons extrait les données IRMf du tractus visuel cortical ventral au format numpy, nous utiliserons ces données pour l'apprentissage automatique. ** Cette fois, nous allons apprendre avec un enseignant, nous avons donc besoin d'une étiquette en plus des données **. Par conséquent, nous allons travailler pour extraire les informations de la catégorie que le sujet recherchait au moment de l'expérience dans les données téléchargées sous forme d'étiquette.
import pandas as pd
#Lisez les informations sur le résultat de l'action. Il y a autant d'étiquettes que de scans
behavioral = pd.read_csv('~/Desktop/nilearn/subj2/labels.txt', delimiter=' ')
Résultat de sortie
labels chunks
0 rest 0
1 rest 0
2 rest 0
3 rest 0
4 rest 0
5 rest 0
6 scissors 0
7 scissors 0
8 scissors 0
9 scissors 0
10 scissors 0
11 scissors 0
12 scissors 0
13 scissors 0
14 scissors 0
... ... ...
1426 cat 11
1427 cat 11
1428 cat 11
1429 cat 11
1430 cat 11
1431 cat 11
1432 rest 11
1433 rest 11
1434 rest 11
1435 rest 11
1436 rest 11
1437 scissors 11
1438 scissors 11
1439 scissors 11
1440 scissors 11
1441 scissors 11
1442 scissors 11
1443 scissors 11
1444 scissors 11
1445 scissors 11
1446 rest 11
1447 rest 11
1448 rest 11
1449 rest 11
1450 rest 11
1451 rest 11
1452 rows × 2 columns
Puisque cette tâche est une tâche de reconnaissance visuelle comme mentionné ci-dessus, l'étiquette montre les conditions expérimentales (catégorie de l'objet montré au sujet). les morceaux représentent le nombre de sessions (car si vous laissez le sujet faire la tâche pour toujours, la charge sur le corps n'est pas étrange, généralement vous effectuez des tâches similaires plusieurs fois tout en vérifiant votre condition physique aux pauses de session. Je vais continuer).
conditions = behavioral['labels']
Résultat de sortie
0 rest
1 rest
2 rest
3 rest
4 rest
5 rest
6 scissors
7 scissors
8 scissors
9 scissors
10 scissors
11 scissors
12 scissors
13 scissors
14 scissors
15 rest
16 rest
17 rest
18 rest
19 rest
20 rest
21 face
22 face
...
1431 cat
1432 rest
1433 rest
1434 rest
1435 rest
1436 rest
1437 scissors
1438 scissors
1439 scissors
1440 scissors
1441 scissors
1442 scissors
1443 scissors
1444 scissors
1445 scissors
1446 rest
1447 rest
1448 rest
1449 rest
1450 rest
1451 rest
Name: labels, Length: 1452, dtype: object
Cette fois, à la suite du tutoriel, nous n'utiliserons que les conditions expérimentales pour les chats et les visages. Il existe de nombreux autres libellés, donc si vous êtes intéressé, essayez de vérifier dans d'autres catégories.
#Récupérez le boxel correspondant au nombre de scans dans les conditions expérimentales où le visage et le chat sont présentés à l'image cérébrale masquée.
condition_mask = conditions.isin(['face','cat'])
fmri_masked = fmri_masked[condition_mask]
print(fmri_masked.shape)
#Appliquer les mêmes conditions aux étiquettes
conditions = conditions[condition_mask]
print(conditions.shape)
Décodage avec SVM
Cette fois, nous utiliserons la boîte à outils d'apprentissage automatique de scikit-learn pour les données fmri_masked récupérées. J'utilise le noyau linéaire SVM comme décodeur. Comme il s'agit de scikit-learn, vous pouvez tout faire, de l'apprentissage à la prédiction en quelques lignes.
from sklearn.svm import SVC
svc = SVC(kernel='linear')
#Apprenez et faites des prédictions
svc.fit(fmri_masked, conditions)
prediction = svc.predict(fmri_masked)
Cependant, cela est totalement inutile pour l'analyse, nous allons donc effectuer une validation croisée (évidemment).
Effectuer une validation croisée pour calculer le score prévu
Effectuer la vérification d'intersection KFold
from sklearn.model_selection import KFold
import numpy as np
cv = KFold(n_splits=5)
scores = []
for train, test in cv.split(X=fmri_masked):
svc.fit(fmri_masked[train], conditions.values[train])
prediction = svc.predict(fmri_masked[test])
scores.append((prediction == conditions.values[test]).sum()
/ float(len(conditions.values[test])))
#Performance prédictive moyenne
print(np.mean(scores))
Résultat de sortie:
0.7628964059196617
Il est également possible d'utiliser la fonction cross_val_score pour effectuer un traitement distribué en utilisant plusieurs CPU sur la machine pour le traitement du calcul de l'évaluation pour chaque sous-ensemble. (Il peut être spécifié sous la forme n_jobs = n, et s'il est défini sur -1, le traitement parallèle peut être effectué en utilisant tous les processeurs disponibles sur la machine).
from sklearn.model_selection import cross_val_score
cv_score = cross_val_score(svc, fmri_masked, conditions)
print(cv_score)
Résultat de sortie:
[0.86363636 0.76744186 0.74418605 0.69767442 0.81395349]
** Cependant, cela est encore insuffisant pour la vérification. La raison en est que nous n'avons pas évalué les performances compte tenu de l'impact de chaque session. ** **
Évaluer la précision des prédictions pour chaque session
Les données IRMf sont acquises session par session et le bruit est auto-corrélé pour chaque session spécifique. Par conséquent, lorsque vous effectuez une validation croisée avec ces effets à l'esprit, il est nécessaire de faire des prédictions entre les sessions.
from sklearn.model_selection import LeaveOneGroupOut
cv = LeaveOneGroupOut()
cv_score = cross_val_score(svc,
fmri_masked,
conditions,
cv=cv,
groups=session_label,
)
#Obtenez la précision des prédictions dans les données de vérification pour chaque session
print(cv_score)
Résultat de sortie:
[0.55555556 1. 0.66666667 0.66666667 0.77777778 0.72222222
0.88888889 0.38888889 0.66666667 0.5 0.77777778 0.66666667]
Estimer le poids du modèle entraîné
Enfin, le poids du modèle est estimé et affiché. Pour ce faire, nous convertissons d'abord les poids en images NIfTY.
coef_ = svc.coef_
#Il se présente sous la forme d'un tableau numpy, avec un coefficient par boxel.(poids)avoir
print(coef_.shape)
Résultat de sortie:
(1, 464)
coef_img = masker.inverse_transform(coef_)
print(coef_img)
Résultat de sortie:
<class 'nibabel.nifti1.Nifti1Image'>
data shape (40, 64, 64, 1)
affine:
[[ -3.5 0. 0. 68.25 ]
[ 0. 3.75 0. -118.125]
[ 0. 0. 3.75 -118.125]
[ 0. 0. 0. 1. ]]
metadata:
<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr : 348
data_type : b''
db_name : b''
extents : 0
session_error : 0
regular : b''
dim_info : 0
dim : [ 4 40 64 64 1 1 1 1]
intent_p1 : 0.0
intent_p2 : 0.0
intent_p3 : 0.0
intent_code : none
datatype : float64
bitpix : 64
slice_start : 0
pixdim : [-1. 3.5 3.75 3.75 1. 1. 1. 1. ]
vox_offset : 0.0
scl_slope : nan
scl_inter : nan
slice_end : 0
slice_code : unknown
xyzt_units : 0
cal_max : 0.0
cal_min : 0.0
slice_duration : 0.0
toffset : 0.0
glmax : 0
glmin : 0
descrip : b''
aux_file : b''
qform_code : unknown
sform_code : aligned
quatern_b : 0.0
quatern_c : 1.0
quatern_d : 0.0
qoffset_x : 68.25
qoffset_y : -118.125
qoffset_z : -118.125
srow_x : [-3.5 0. 0. 68.25]
srow_y : [ 0. 3.75 0. -118.125]
srow_z : [ 0. 0. 3.75 -118.125]
intent_name : b''
magic : b'n+1'
Enregistrez l'image au format NIfTY créée au format nii.gz
# nii.Enregistrer au format de fichier gz
coef_img.to_filename('haxby_svc_weights.nii.gz')
Enfin, visualisons le coefficient de poids sous forme d'image cérébrale.
#Tracez en fait le poids de SVM
from nilearn.plotting import plot_stat_map, show
plot_stat_map(coef_img, bg_img=anat_filename,
title="SVM weights", display_mode='yx')
show()
Avec ce qui précède, nous avons effectué une série d'analyses de données IRMf à l'aide de l'apprentissage automatique.
Résumé
J'avais l'intention de publier cette analyse sur mon blog personnel, mais je n'ai trouvé personne qui analysait les données IRMf avec Qiita, alors je l'ai publiée ici (des blogs individuels sont publiés sur mon profil).
Image de données IRMf ... Cela peut sembler difficile juste de l'entendre, mais au final c'est juste plus de bruit qu'une image générale, une dimension de plus, et la suppression de la partie qui gère les formats spéciaux L'analyse d'image. L'un des objectifs est de pouvoir analyser les données maniaques de n'importe qui, et je les publie via mon blog personnel et Qiita, alors j'espère que vous serez intéressé par l'analyse des données que vous touchez rarement à travers cette activité.
Si cela ne vous dérange pas, LG ... Merci pour votre coopération! Lol
Une note nommée Addendum
Cette fois, nous analyserons en utilisant des données brutes de séries chronologiques IRMf, mais il est nécessaire de convoluer ce qu'on appelle à l'origine la fonction HRF et de la convertir de la présentation des tâches (activité nerveuse) à l'hémodynamique (réaction BOLD). De plus, afin d'obtenir une carte d'activité essai par essai (carte bêta), il est nécessaire d'ajuster au 1er niveau en utilisant un modèle linéaire général. Cela permet de faire des prédictions entre des conditions liées pour chaque essai. En fait, divers processus de traitement sont nécessaires pour l'analyse réelle des données IRMf, mais cette fois, il s'agit d'une introduction visant à appliquer l'apprentissage automatique aux données IRMf, j'ai donc décidé d'effectuer le traitement susmentionné (en fait) Si vous prenez autant de mesures, cela ne rentrera pas dans un seul article (transpiration).
Veuillez noter (je voudrais analyser quelque chose sérieusement s'il y a une opportunité)