SVILUPPO DI ALGORITMI PER LA
CLASSIFICAZIONE DI CROMOSOMI IN
IMMAGINI IN METAFASE
Enrico Fucigna
Relatore: Prof. Alfredo Ruggeri, Università di Padova
Correlatore: Enrico Grisan
Aprile 2006
2
INDICE
Titolo: Sviluppo di algoritmi per la classificazione di
cromosomi in immagini in metafase
SOMMARIO................................................................................................... 5
1
INTRODUZIONE .................................................................................... 6
1.1
2
Scopo della tesi..........................................................................................................6
TECNICHE DI ANALISI CROMOSOMICA............................................. 7
2.1
2.1.1
2.1.2
2.1.3
2.1.4
2.1.5
2.1.6
I cromosomi...............................................................................................................7
Morfologia .......................................................................................................................... 7
Cariogramma....................................................................................................................... 8
Banding pattern ................................................................................................................... 9
Ciclo cellulare ................................................................................................................... 11
Risoluzione........................................................................................................................ 11
Anomalie cromosomiche .................................................................................................. 11
2.2
2.2.1
International System for Human Cytogenetic Nomenclature ............................13
Nomenclatura .................................................................................................................... 13
2.3
2.3.1
2.3.2
2.3.3
2.3.4
2.3.5
2.3.6
Analisi computerizzata dei cromosomi .................................................................14
Analisi automatica............................................................................................................. 14
Pre-processing................................................................................................................... 15
Segmentazione .................................................................................................................. 15
Estrazione parametri.......................................................................................................... 16
Selezione dei parametri ..................................................................................................... 17
Classificazione .................................................................................................................. 17
3
SVILUPPO DI ALGORITMI PER LA CLASSIFICAZIONE................... 19
3.1
Materiali utilizzati ..................................................................................................19
3.2
3.2.1
3.2.2
3.2.3
Algoritmo di image enhancement..........................................................................21
Decomposizione e ricostruzione wavelet .......................................................................... 21
Multiscale pointwise product ............................................................................................ 22
Enhancement algorithm .................................................................................................... 23
3.3
3.3.1
3.3.2
3.3.3
Estrazione di parametri .........................................................................................24
Lunghezza ......................................................................................................................... 26
Density profile................................................................................................................... 26
Centromeric index............................................................................................................. 27
3.4
3.4.1
3.4.2
Calcolo dell’ asse cromosomico .............................................................................29
Algoritmo proposto ........................................................................................................... 31
Osservazioni...................................................................................................................... 33
3.5
Algoritmi che sfruttano il calcolo dell’ asse cromosomico ..................................34
3
3.5.1
3.5.2
3.5.3
Estrazione del banding pattern ...........................................................................................34
Calcolo della distanza del bordo cromosomico dall’ asse..................................................35
Calcolo del centromeric index ...........................................................................................36
4
RISULTATI............................................................................................39
5
CONCLUSIONE....................................................................................46
BIBLIOGRAFIA............................................................................................47
APPENDICE.................................................................................................49
Abbreviazioni ......................................................................................................................................50
4
Sommario
5
1 Introduzione
1.1 Scopo della tesi
Lo scopo di questo studio è la realizzazione di alcuni algoritmi di interesse nel
campo della classificazione di cromosomi.
6
2 Tecniche di analisi cromosomica
2.1 I cromosomi
I cromosomi sono i trasportatori fondamentali dell’ informazione genetica degli
individui. Il DNA umano è organizzato infatti in cromosomi attivi nelle cellule, le
quali si dividono, si differenziano e guidano lo sviluppo dell’ intero organismo.
Ci sono 46 cromosomi in una normale cellula umana, visivamente i cromosomi
possono essere catalogati in una di 24 classi (etichettate 1-22 X, Y) . Le ultime due
classi individuano i cromosomi sessuali, che si presentano entrambi nella classe X
per le femmine ed uno nella classe X, l’ altro nella classe Y per i maschi. Tutte le
altre classi negli individui normali contengono una coppia di cromosomi identici
detti omologhi. Ciascun cromosoma componente la coppia di omologhi viene
ereditato da ciascun genitore. I cromosomi non sessuali sono detti autosomi.
Dall’ osservazione al microscopio risulta che i cromosomi differiscono per
dimensione e morfologia. Ogni cromosoma, lungo la sua lunghezza spesso presenta
una costrizione chiamata centromero, che è importante per il comportamento dei
cromosomi dal punto di vista cellulare. La collocazione del centromero del
cromosoma in una delle quattro classiche posizioni è utile per la classificazione dei
cromosomi eucarioti. Un cromosoma metacentrico ha il centromero localizzato al
centro del cromosoma, che risulta così diviso in due bracci di lunghezza quasi
uguale. I cromosomi sub-metacentrici hanno un braccio più lungo dell’ altro, i
cromosomi acrocentrici mostrano un braccio con una breve appendice su cui spesso
si osserva un lieve rigonfiamento (denominato “satellite”) e i cromosomi telocentrici
hanno soltanto un braccio, poiché il centromero è situato ad una estremità (Figura
2.1).
Figura 2.1
Classificazione generale dei cromosomi in metacentrici,
submetacentrici, acrocentrici e telocentrici, basata sulla posizione del
centromero.
2.1.1 Morfologia
La lunghezza relativa dei cromosomi umani varia in modo considerevole, ed
unitamente alla posizione del centromero è di aiuto nell’ identificazione dei singoli
7
cromosomi. Gruppi di cromosomi umani con morfologia simile sono designati con
una lettera da A a G:
Gruppo A (1-3): Cromosomi larghi, metacentrici, facilmente distinguibili dagli
altri per le dimensioni e la posizione del centromero.
Gruppo B (4-5): Cromosomi larghi, submetacentrici.
Gruppo C (6-12,X): Cromosomi di taglia media, metacentrici o submetacentrici.
Gruppo D (13-15): Cromosomi di taglia media, acrocentrici, con satelliti.
Gruppo E (16-18): Cromosomi relativamente corti, metacentrici o submetacentrici.
Gruppo F (19-20): Cromosomi corti, metacentrici.
Gruppo G (21-22,Y): Cromosomi corti, acrocentrici, con satelliti. Y non possiede satellite.
Questa era l' unica classificazione possibile in assenza di informazioni oltre alla
lunghezza cromosomica ed alla posizione centromerica. Quando i diversi cromosomi
sono colorati in modo uniforme, è difficile infatti distinguerli in modo non ambiguo
basandosi soltanto su dimensione e forma.
2.1.2
Cariogramma
Un insieme completo di tutti i cromosomi metafasici di una cellula viene definito
cariogramma. Vengono utilizzati di consueto cromosomi metafasici per la
generazione di un cariogramma in quanto questo è il loro stato di maggior
condensazione nella mitosi (divisione cellulare). I cromosomi si presentano quindi in
questa fase nella loro forma più compatta e di più facile osservazione al microscopio
dopo appropriata colorazione. È scelta adottata disporre i cromosomi in ordine
decrescente secondo la dimensione e la posizione del centromero. In un cariotipo
umano i cromosomi sono numerati per una più facile identificazione. Per
convenzione il paio di cromosomi omologhi più lungo è designato come 1, il
successivo più lungo 2 e così via. (Figura 2.2). Benché in realtà il cromosoma 21 sia
più piccolo del cromosoma 22, è chiamato 21 per ragioni storiche.
8
(a)
(b)
Figura 2.2
Immagine di cromosomi estratti da una cellula in metafase di essere umano di sesso maschile,
colorati con Giemsa in (a) ed in (b) il corrispondente cariogramma.
2.1.3
Banding pattern
Per semplificare il riconoscimento dei cromosomi è stata sviluppata una serie di
tecniche che colorano alcune regioni cromosomiche o bande più intensamente di
altre. I bandeggi sono specifici per ogni cromosoma e questo consente di identificare
inequivocabilmente ogni cromosoma del cariotipo.
Una banda è definita come quella parte del cromosoma che è chiaramente
distinguibile dai suoi segmenti adiacenti perché appare più chiara o più scura
attraverso una o più tecniche di bandeggio. Bande che appaiono chiare con una
colorazione possono apparire scure utilizzando una tecnica diversa. I metodi di
colorazione pubblicati per primi utilizzano quinacrine dihydrocloride, un colorante
fluorescente che si lega preferenzialmente al DNA ricco di AT e osservati con
fluorescenza UV per produrre un banding pattern fluorescente e le bande risultanti
vengono definite bandeggio Q. Un pattern quasi identico viene ottenuto
sottoponendo i cromosomi a digestione controllata con tripsina e utilizzando in
seguito il colorante Giemsa che si lega anch’esso al DNA ; le bande scure
corrispondenti vengno nominate bandeggio G, le bande chiare si dicono G-negative;
le bande G scure individuano gli stessi segmenti cromosomici delle bande
fluorescenti Q. Le bande R sono essenzialmente il contrario delle bande G. I
cromosomi vengono denaturati al calore in soluzione salina prima di essere colorati
con Giemsa. Il trattamento al calore denatura il DNA ricco di AT e le bande R sono
Q-negative. Il metodo di bandeggio T identifica un gruppo di bande R concentrate in
particolare ai telomeri. Le bande T sono tra le bande R quelle che si colorano più
intensamente e vengono visualizzate sottoponendo i cromosomi a temperature
particolarmente elevate prima di colorarli con Giemsa, oppure utilizzando una
combinazione di coloranti e fluorocromi. Il bandeggio C si pensa sia in grado di
rivelare la presenza di eterocromatina costitutiva. Una nuova tecnica di colorazione
consiste nella Fluorescence In-Situ Hybridization (FISH) in cui il DNA viene
ibridato attraverso una sonda adatta marcata con fluorescenza consentendo di
9
accrescere significativamente la sensibilità e la risoluzione del’ osservazione (Figura
2.3). Recentemente è stata utilizzata un’ ulteriore tecnica nominata M-FISH in cui
viene generata un’ immagine composita di cinque differenti tipi di fluorescenza e a
ciascun cromosoma viene assegnato uno pseudocolore a seconda della sua
composizione di fluorescenze.
(a)
(b)
(c)
Figura 2.3
Tre esempi di immagini di cellule diverse in metafase analizzate con una colorazione differente:nella
cellula in (a) è stato messo in evidenza il bandeggio G, in (b) sono state ricavate le bande Q, in (c)
vengono evidenziati i centromeri con una colorazione FISH che utilizza la fluorescenza DAPI.
Cromosomi con bandeggi Q, G, R più lunghi mostrano un numero di bande maggiore
rispetto a quelli più corti. L’ aspetto dei cromosomi al microscopio dipende infatti
dallo stadio del ciclo di divisione cellulare in cui vengono visti. Il maggior numero di
bande presente in stadi di maggiore elongazione fornisce una descrizione con una
risoluzione più elevata della struttura del cromosoma, più vantaggiosa ai fini dell’
analisi. L` osservazione in questa fase viene tuttavia resa più difficoltosa a causa
della maggiore complessità del banding pattern (bandeggio) e del fatto che
cromosomi più lunghi tendono a toccarsi e sovrapporsi molto di più rispetto a
cromosomi più corti.
10
2.1.4
Ciclo cellulare
Per la maggior parte del ciclo cellulare (interfase), i singoli cromosomi non possono
essere distinti. La mitosi, ossia il processo di duplicazione cellulare viene solitamente
distino in quattro fasi: profase, metafase, anafase e telofase. All’ inizio della profase,
i cromosomi sono molto allungato, simili a sottili stringhe longilinee, che in seguito
cominciano ad avvolgersi, per cui appaiono al microscopio più corti e spessi. Nella
profase avanzata si osserva che ciascun cromosoma è costituito da due subunità
longitudinali definite cromatidi fratelli. La metafase inizia quando la membrana
nucleare si è dissolta completamente, i cromosomi si scindono originando due copie
esatte, unite al centromero, giusto prima che la divisione cellulare abbia luogo. L’
anafase ha inizio quando i centromeri uniti dei due cromatidi fratelli si separano
dando origine a due cromosomi figli. La citocinesi (divisione cellulare generalmente
inizia negli stadi tardivi dell’ anafase). Durante la telofase si compie il processo di
divisione nucleare ed al termine di questo stadio la cellula possiede due nuclei.
Benché non significativo dal punto di vista biologico è comune riferirsi ad uno stadio
intermedio di contrazione sito tra profase e metafase indicato con il nome di
prometafase.
2.1.5
Risoluzione
Operativamente questi stadi del ciclo di divisione cellulare vengono definiti dal
numero di bande visibili nella cellula. Le cellule in metafase possiedono
approssimativamente 450 bande e si definiscono in prometafase quando hanno 550
bande e oltre (vedi Figura 2.4). Una cellula può essere definita in profase soltanto se
almeno 850 bande sono visibili. Il numero complessivo di bande visibili in un’
immagine cellulare viene definito risoluzione del bandeggio.
2.1.6
Anomalie cromosomiche
I cromosomi possono presentare delle anomalie ossia cambiamenti che producono un
un’ alterazione visibile dei cromosomi, quali di queste siano visibili dipende dalla
tecnica utilizzata. La maggior parte delle aberrazioni cromosomiche è prodotta da
errori nella riparazione di rotture cromosomiche, da errata ricombinazione o da errata
segregazione dei cromosomi durante la mitosi. Le anomalie cromosomiche rientrano
per lo più in due categorie: la anomalie numeriche e quelle strutturali. Tra le
anomalie numeriche sono note triploidia (presenza di cromosomi in eccesso),
trisomia (presenza di autosomi in eccesso es. tre cromosomi 21 indicano la sindrome
di Down), monosomia (mancanza di un cromosoma ), mosaicismo, tra quelle
strutturali si conoscono delezione (perdita di un pezzo di cromosoma), duplicazione
inserzione (un frammento staccato si attacca in un altro punto del cromosoma),
cromosoma ad anello, marcatore e traslocazione (scambio di frammenti tra
cromosomi). Il cariogramma si presenta quindi con tutti i cromosomi con i bandeggi
standard soltanto in individui sani.
11
(a)
(b)
(c)
Figura 2.4
Tre immagini in banda Q di cariogrammi raccolti da cellule di essere umano di sesso maschile a
differenti livelli di risoluzione: l’ immagine in (a) rappresenta un cariogramma estratto da una cellula
in metafase iniziale. I cromosomi mostrano approssimativamente 500 bande e presentano una forma
allungata e sottile. L’ immagine in (b) contiene un cariogramma in cui i cromosomi sono in metafase
intermedia e mostrano una struttura accorciata e compatta, in cui si osservano circa 300 bande. Nell’
immagine in (c) i cromosomi si trovano in metafase avanzata, i due cromatidi fratelli sono in via di
separazione ma rimangono uniti all’ altezza del centromero. Non sono visibili in questo caso più di 70
bande.
12
2.2 International System for Human Cytogenetic Nomenclature
2.2.1
Nomenclatura
Il Sistema Internazionale per la Nomenclatura Umana (ISCN [1]) è stabilito dallo
Standing Committee of Human Cytogenetic Nomenclature. La terminologia
fondamentale per i cromosomi bandeggiati fu decisa in un convegno a Parigi nel
1971 ed è spesso indicata come “nomenclatura di Parigi”. Le localizzazioni sul
braccio corto sono indicate da una p (petite), quelle sul braccio lungo con una q
(queue). Si considera che i bracci di ciascun cromosoma siano suddivisi in differenti
regioni che vengono indicate come p1, p2, p3 etc. e q1, q2, q3 etc., contando a partire
dal centromero e andando verso i telomeri. Le varie regioni sono delimitate da punti
di riferimento specifici, con caratteristiche morfologiche costanti e distinte, che
comprendono le estremità dei bracci cromosomici, il centromero e alcune bande . Le
regioni vengono divise in bande indicate come p11, p12, p13, etc., ed in sottobande
p11.1, p11.2, etc., con una numerazione che va sempre dal centromero verso i
telomeri. Il centromero è indicato come “cen” ed il telomero con “ter”. Xq
prossimale significa il segmento del braccio lungo dell’ X che è più vicino al
centromero, mentre 2p distale significa la porzione del braccio corto del cromosoma
2 che è più distante dal centromero e, quindi, più vicina al telomero.
La posizione del centromero viene definita attraverso il centromeric index ossia il
rapporto tra il braccio più corto e l’ intera lunghezza del cromosoma.
Figura 2.5
Cromosoma in banda G con indicazioni relative al centromero e ai
bracci p e q. per definizione il braccio più corto viene designato p.
la freccia indica la posizione del centromero
Ideogrammi. L’ ISCN ha definito inoltre un insieme di modelli chiamati
ideogrammi volti a fornire una descrizione visuale dei bandeggi cromosomici. Per
convenzione questi disegni rappresentano le bande G .
13
Figura 2.6
Ideogrammi del bandeggio G per cromosomi umani normali a 3 differenti livelli di
risoluzione. Da sinistra, i cromosomi in ogni gruppo rappresentano un cariogramma di
rispettivamente 400, 550 ed 850 linee di bandeggio. Le bande G scure corrispondono a
bande Q fluorescenti, con l’ eccezione delle regioni variabili (in grigio). La numerazione
delle bande R è esattamente la stessa con la condizione di complementare bande chiare e
scure. Negli ideogrammi viene rappresentata anche la posizione del centromero
2.3 Analisi computerizzata dei cromosomi
Ai fini dello sviluppo di un cariogramma una cellula viene fotografata al microscopio
ottico durante lo stadio di metafase, dall’ immagine ottenuta vengono ritagliati i
cromosomi e disposti nella locazione appropriata, in accordo con la classificazione
visiva operata dal citologo
Il cariogramma è uno strumento molto utile per riconoscere deviazioni dalla
normale struttura di una cellula ma tuttora l’ analisi dei cromosomi e la relativa
classificazione vengono eseguite manualmente nella maggior parte dei laboratori di
citogenetica, con un procedimento ripetitivo, dispendioso in termini temporali e
molto costoso. Il Canadian Workload Measurement System [9] impiega 465 min per
la creazione del cariogramma degli amniociti, una della operazioni più comuni in
citogenetica. La gran parte del tempo viene dedicata al microscopio, una fase
tediante, stressante per l’ occhio e che richiede una meticolosa attenzione ai dettagli,
oltre a necessitare di personale molto qualificato.
2.3.1
Analisi automatica
L’ analisi automatica dei cromosomi è stata proposta per la prima volta nel 1960. Nei
tre decenni passati sono stati compiuti molti sforzi al fine di risolvere il problema
della classificazione completamente automatizzata dei cromosomi. Considerando
tuttavia gli innumerevoli possibili aspetti con cui si presentano le immagini
cromosomiche, bisogna prendere atto che la complessità di questo problema ascende
ad un livello elevatissimo.
14
I software di creazione automatica del cariogramma disponibili attualmente in
commercio si rivelano essere di grande aiuto ma ancora inadeguati. Alcune delle
ragioni che spiegano queste modeste prestazioni vanno ricercate nel carente utilizzo
della conoscenza e della competenza degli esperti, nell’ insufficiente capacità di
eseguire confronti o eliminazioni tra cromosomi nella stessa metafase. In aggiunta i
sistemi attuali richiedono sempre l’ interazione dell’ operatore per separare
cromosomi che si toccano o sovrappongono e per verificare il risultato della
classificazione.
L’ analisi al computer dei cromosomi consiste in primo luogo di pre-processing dell’
immagine ottenuta al microscopio, segmentazione,(eventualmente seguita da un
processing intermedio), estrazione dei parametri, selezione degli stessi ed infine
classificazione.
2.3.2
Pre-processing
La fase di pre-processing ha come obiettivo il miglioramento della qualità dell’
immagine cellulare attraverso tecniche di riduzione del rumore, ottimizzazione dei
bordi del cromosoma, aumento del contrasto. Procedure di questo tipo si possono
riscontrare ad esempio nei lavori di Moradi e Setarehdan [10], in cui, prima di
effettuare il calcolo di qualsiasi parametro dipendente dalle condizioni di
illuminazione, viene eseguita una normalizzazione dell’ intensità dell’ immagine
attraverso la histogram stretching technique. Anche nello studio di Shunren et al.
[19], a precedere il successivo stadio di segmentazione viene modificato l’
istogramma dell’ immagine eseguendo uno stretching dei livelli di grigio. Uno dei
lavori migliori nell’ambito della realizzazione di un software destinato al image
enhancement di immagini cromosomiche è quello compiuto da Wang e Wu [2],[3],
in cui viene esaminato quantitativamente l’ effetto del miglioramento dell’ immagine
sulle prestazioni del pattern recognition per immagini cromosomiche e si arriva a
verificare come l’ image enhancement comporti un significativo incremento nella
precisione dell’identificazione. L’ algoritmo che introducono si basa su filtri
differenziali che facilitano l’ estrazione delle caratteristiche geometriche multiscala
dell’ immagine cromosomica. Considerando che caratteristiche geometriche come i
bordi sono correlate a scale differenti nella rappresentazione in wavelets, viene
utilizzato un multiscale point-wise product (MPP) per caratterizzare la correlazione
degli aspetti peculiari dell’ immagine a diverse scale.
2.3.3
Segmentazione
Il proposito della fase di segmentazione verte invece sull’ isolamento dei cromosomi
in metafase dallo sfondo, da nuclei cellulari indivisi, da materiale biologico
irrilevante presente all’ interno dell’ immagine e altro. In materia è stato scritto molto
in letteratura e generalmente si conviene che il primo passo della segmentazione
consti di un processo a soglia per la separazione dello sfondo dagli oggetti di
interesse. Karvelis et al. nel loro studio [20] ricorrono al processo di thresholding di
Otsu così come Guimaraes et al. [21] e Piper e Granum [13]; questo consente di
ricondursi ad un’ immagine binaria da cui si ricavano facilmente i cromosomi singoli
ma si incontrano ingenti difficoltà nella separazione di cromosomi che si toccano o
sovrappongono.
15
Tra le ricerche compiute nell’ intento di risolvere questo problema critico sono di
notevole interesse quelle effettuate da Karvelis et al. [20] (relativamente tuttavia a
gruppi di cromosomi che si toccano soltanto), ove si esegue una iniziale
segmentazione utilizzando la watershed transform ed in seguito questa
trasformazione viene applicata ricorsivamente ad ogni area segmentata finché non
vengono generate ulteriori aree. Infine vengono prodotte le aree segmentate
definitive calcolando tutti i percorsi che partono da punti ad alta concavità.
Agam e Dinstein [22] propongono un criterio che tratta il problema della
separazione di clusters di cromosomi sovrapposti come un problema di
identificazione, Shunren et al. [19] introducono ben due metodi di incisione dei
cromosomi sovrapposti: il metodo della caratteristica del contorno, che sfrutta
elementi del contorno, ad esempio la presenza di punti concavi sul bordo di un
cluster costituito da cromosomi sovrapposti, ed il metodo della trasformata di
Fourier, che impiega un filtro passa basso di Fourier per elaborare l’ immagine
cromosomica e creare una nuova immagine in cui il tono di grigio degli oggetti viene
ridistribuito con il gradiente. Da questa utilizzando a seguire un processo a soglia l’
interno del cromosoma viene rivelato. Ulteriori lavori rilevanti in questo senso sono
da riferire a Lerner et al. [7] [6], che introduce il metodo CPOOS (Classification
Driven Partially Occluded Object Segmentation) che viene applicato, in seguito alla
classificazione degli oggetti isolati, sui clusters che non è stato possibile separare. Il
metodo prende i punti più concavi lungo la curvatura dell’ immagine e
congiungendoli ricava la possibili linee di separazione dei cromosomi; infine le linee
di taglio più probabili vengono scelte come effettive.
Charters e Graham, [15], cercano di svincolare le sovrapposizioni attraverso
modelli del bandeggio cromosomico che possono essere sottoposti ad un training. I
modelli consistono di templates dei profili del banding pattern di lunghezza inferiore
a quella cromosomica. I segmenti cromosomici candidati vengono classificati in
accordo con il loro confronto con l’ intero set di templates, e uniti in base alla
classificazione. In uno studio successivo [15], unitamente ai modelli parziali di
bandeggio cromosomico vengono considerati modelli addestrabili della forma
introducendo una parametrizzazione definita Chord Distribution Model (CDM). La
percentuale complessiva di cromosomi identificati correttamente migliora
sensibilmente.
In seguito alla segmentazione può venire introdotta una ulteriore fase di processing
dell’ immagine ottenuta in precedenza, con l’ obiettivo di perfezionare quegli aspetti
dell’ immagine che facilitano il riconoscimento dei cromosomi, quali ad esempio il
banding pattern. Da noverare in questo senso l’ esempio offerto da Shunren et al.
[19] nel loro scritto, ove alla fase di segmentazione segue un’ operazione
matematica di opening morfologico, atto ad assottigliare gli oggetti e rimuovere il
rumore ancora presente.
2.3.4
Estrazione parametri
Quanto segue le fasi descritte in precedenza viene comunemente ricondotto allo
stadio di estrazione delle features. Questo punto richiede di effettuare la misura dei
parametri cromosomici per estrarne le caratteristiche che in maggior misura si
prestano a discriminare le differenti classi cromosomiche. Il risultato poi della
selezione delle caratteristiche, se applicata, consente di ricondurre il cromosoma ad
una rappresentazione più sintetica che tuttavia contenga la maggior parte dell’
16
informazione cromosomica. A partire da quest’ ultima descrizione del cromosoma
viene eseguita la classificazione vera e propria, usualmente valendosi di un metodo
statistico.
L’ utilizzo di parametri piuttosto che dell’ immagine stessa rende il processo di
classificazione più semplice, veloce ed accurato. Sono state utilizzate numerose
caratteristiche differenti per descrivere un cromosoma, ad esempio nella
pubblicazione di Piper e Granum del 1989 [13] vengono vagliati diversi parametri,
taluni misurabili direttamente dall’ immagine cromosomica (area, densità relativa,
perimetro del convex hull) altri invece ricavati in seguito al calcolo dell’ asse
(lunghezza, banding pattern) ed ulteriori che richiedono, oltre al computo dell’ asse,
l’ informazione circa la polarità del cromosoma, che può essere ottenuta
determinando anche la locazione del centromero (Centromeric Index, ossia il
rapporto tra il braccio più corto del cromosoma e l’ intera lunghezza cromosomica).
Sweeney at al. [23] calcolano la lunghezza unitamente ai coefficienti wavelet e di
Fourier derivati dal profilo densitometrico per formare un vettore di parametri per
ciascun cromosoma, Delshadpour [24] aggiunge a lunghezza e banding pattern la
decomposizione dei livelli di grigio o componenti della trasformata di Fourier
bidimensionale, Guimaraes [21] considera invece il rilevamento della forma del
cromosoma, ne ricava il contorno e quindi la relativa decomposizione attraverso la
wavelet packet transform. Nonostante molte possibili analisi effettuate su immagini
cromosomiche siano state prese in considerazione, quelle maggiormente
discriminanti ai fini della classificazione sembrano essere tuttavia la lunghezza ed il
centromeric index (ossia due caratteristiche geometriche), ed il density profile
(integrale o densità media lungo sezioni perpendicolari al medial axis del
cromosoma), come riportato negli studi di Lerner et al. [7] [8], in quelli di Moradi et
al. [10], [11], [12], ed in quelli di Cho [17].
2.3.5
Selezione dei parametri
Una volta estratti quelle features del cromosoma che si ritengono più utili alla
classificazione è possibile cercare di ridurne il numero mantenendo intatto o
diminuendo lievemente il contenuto di informazione attraverso la feature selection,
che può essere intesa come la ricerca, tra tutte le possibili trasformazioni dello spazio
dei parametri, di quella che preserva la separabilità tra le classi in uno spazio con il
minor numero possibile di dimensioni. Lerner [7] propone lo knock-out algorithm
che valuta, attraverso una matrice che pesa lo scatter within-class e between-class, l’
efficacia di vettori di parametri costituiti di sottoinsiemi del vettore iniziale,
eliminando quindi le features meno discriminanti.
2.3.6
Classificazione
Il problema che si pone al termine della convenzionale procedura di analisi dei
cromosomi è quello della classificazione vera e propria, l’ argomento che è stato
maggiormente esaminato nell’ analisi cromosomica. In quest’ ambito il classificatore
più popolare è una rete neurale con multi-layer perceptron, usualmente addestrata
con un algoritmo di backpropagation come in Delshadpour [24], Lerner [8], Cho
[17], Moradi [10]; altri classificatori sono stati studiati, come gli inferred Markov
network models di Granum e Thomason [14], statistici come in Piper and Granum
17
[13] (maximum likelihood classifier) o ancora fuzzy o nearest neighbor. La maggior
parte di questi classificatori ha due inconvenienti: innanzitutto offre prestazioni di
livello inferiore rispetto ad un citologo esperto (70-80% paragonato a 99.7%) ed in
secondo luogo richiede l’ intervento di un operatore per risolvere casi di erronea
classificazione di cromosomi. Alcune delle ragioni alla base di queste limitazioni
vanno ricercate nell’ uso inadeguato della metodologia di un esperto (ad esempio il
confronto tra cromosomi o l’ eliminazione), l’ utilizzo di parametri limitati nella
qualità o nel numero se raffrontati con il potentissimo meccanismo di sintesi del
cervello umano, l’ impiego delle stesse caratteristiche specifiche per descrivere tutti i
cromosomi e l’ ingenuo tentativo di rappresentare un problema molto complesso, che
comporta il riconoscimento di schemi molto sensibili al contesto ed esibisce una
naturale varietà, attraverso le restrizioni di una rigida struttura matematica.
In ogni caso un completo ed efficiente sistema di classificazione automatica dei
cromosomi è un obiettivo molto ambizioso. Molti dei lavori precedentemente
riportati sono in effetti da ritenere studi di fattibilità, in cui gli autori si collocano in
un contesto limitato (es. analisi di un numero esiguo di classi cromosomiche, o
utilizzo di cromosomi privi di difetti genetici, non sovrapposti, non eccessivamente
curvati ) ed inoltre i cromosomi considerati vengono prelevati da data sets di utilizzo
pubblico, dove i cromosomi vengono colorati con Giemsa il metodo più diffuso per
evidenziare il banding pattern.
18
3 Sviluppo di algoritmi per la classificazione
3.1 Materiali utilizzati
Le immagini che sono state utilizzate sono di tre tipologie:
- immagini in banda Q (cromosomi bianchi su sfondo nero) con risoluzione PAL
(768x576)
(a)
(b)
Figura 3.1
Immagine di cellula in metafase in figura (a) e relativo cariogramma rappresentato in
(b).
(a)
(b)
(c)
Figura 3.2
Immagini cromosomiche in dettaglio: cromosoma preso da immagine di
una cellula in metafase in (a), proveninti da immagini di cariogramma in (b)
ed in (c). In (b) si nota la presenza dello sfondo da cui l’ immagine è stata
tratta, in (c) la presenza di di cromosomi sovrapposti.
Oltre alle immagini metafasiche in banda Q, sono state fornite immagini dei
cariogrammi relativi, ottenuti da classificazione eseguita da esperti citogenetisti. Le
immagini in metafase con bandeggio Q sono quelle maggiormente utilizzate dai
laboratori di Padova in quanto più economiche e semplici da ottenere. Tuttavia
19
presentano spesso condizioni di illuminazione non ottimali. Le immagini di
cariogramma inoltre sono state scontornate manualmente in modo a volte
approssimativo, rendendo l’ operazione di estrazione dei parametri piuttosto
difficoltosa. Un problema presentatosi, che rende critico l’ ottenimento del
bandeggio, consiste nella presenza di overlaps (sovrapposizioni) anche tra
cromosomi del cariogramma.
- Immagini in banda G acquisite da telecamera digitale (1024x1280)
(a)
(b)
Figura 3.3
Cellula in metafase in banda G in (a), relativo cariogramma in (b).
Queste immagini sono state utilizzate in misura minore in quanto, pur essendo il
bandeggio G il più studiato in letteratura, le richieste dei laboratori propendono per
lo studio del bandeggio Q.
- Immagini in banda Q (768x576) con cariogrammi passati attraverso un programma
di cariotipizzazione automatica.
(a)
(b)
Figura 3.4
Cellula in metafase avanzata con banda G (a) e cariogramma filtrato corrispondente (b).
Queste immagini sono state le meno utilizzate in quanto già filtrate (quindi con un
banding pattern modificato) e con cellule in metafase avanzata (bassa risoluzione del
bandeggio).
20
Il software utilizzato per l’ implementazione degli algoritmi è MATLAB release 12,
che fornisce una serie di utili toolboxes per semplificare numerose operazioni di uso
comune (Wavelet Toolbox, Neural Network Toolbox etc.).
3.2 Algoritmo di image enhancement
Il primo degli algoritmi sviluppati i questo studio è un tipo di algoritmo di processing
che cerca di soddisfare la necessità di ottenere una trasformazione dell’ immagine in
metafase che faciliti l’ estrazione del banding pattern.
Questa caratteristica geometrica si presenta come una serie di bande che
attraversano il cromosoma con direzione perpendicolare al medial axis. Tra le varie
tecniche utilizzate, gli operatori differenziali sono stati usati per molto tempo in
quanto semplificano l’ estrazione di importanti features geometriche come i
contorni. Questo tipo di operatori può essere considerato come un filtro passa-alto,
quindi una tecnica di questo genere raffina l’ immagine estraendo l’ informazione ad
alta frequenza. L’ algoritmo implementato prende le mosse dai lavori di Wang e Wu
et al. [2], [3], [4], in cui presentano una struttura formale per la rappresentazione di
immagini utilizzando operatori differenziali multiscala. Sono stati disegnati molti tipi
di wavelets per la rappresentazione multiscala di immagini. Un’ immagine può
essere decomposta e ricostruita a partire dalle sue componenti a scale differenti.
3.2.1 Decomposizione e ricostruzione wavelet
Avvalendosi delle proprietà delle splines si è ottenuto un rapido algoritmo e le
relative rappresentazioni differiscono dai consueti modelli wavelet in quanto sono
completi ed invarianti rispetto alla traslazione. In generale queste rappresentazioni
sono più flessibili rispetto a quelle ortogonali e si prestano maggiormente al
trattamento di immagine in base alla geometria, permettendo di eliminare spiacevoli
artefatti introdotti da modifiche dei coefficienti di espansione delle wavelets
ortogonali (vedi Appendice).
Sfruttando le proprietà delle splines si può ottenere un algoritmo ricorsivo per il
calcolo di queste tre derivate parziali locali lungo scale diadiche {2j}j ∈ Ζ .Il costrutto
teorico che sta alla base dell’ algoritmo può essere ritrovato in [2], [3], [4].
Si ottiene la seguente formula di decomposizione discreta:
S 2 j f = S 2 j −1 f ⊗ (h, h) ↑ 2 j −1
W 1j f = S 2 j −1 f ⊗ ( g ( 2 ) , d ) ↑ 2 j −1
2
(3.1)
W 2j f = S 2 j −1 f ⊗ (d , g ( 2 ) ) ↑ 2 j −1
2
W 3j f = S 2 j −1 f ⊗ ( g (1) , g (1) ) ↑2 j −1
2
dove I ⊗ (h, g ) ↑ 2 j −1 rappresenta la convoluzione separabile delle righe e delle
colonne dell’ immagine con i filtri 1-D [h]↑2 j −1 e [ g ]↑ 2 j −1 rispettivamente. Il simbolo d
21
denota il filtro di Dirac, il cui impulso vale uno nell’ origine e zero altrove. L’
operazione di upsampling di una sequenza {h(n)} di un fattore intero m è definita
come [h]↑ m .
L’ immagine può essere ricostruita seguendo la formula di ricostruzione discreta
S 2 j −1 f = W21j f ⊗ ( g~ ( 2 ) , u ) ↑ 2 j −1 + W22j f ⊗ (u, g~ ( 2 ) ) ↑2 j −1 +
+ W23j f ⊗ ( g~ (1) , g~ (1) ) ↑2 j −1 + S 2 j f ⊗ (h, h) ↑ 2 j −1
⎛ 2n + 2 ⎞
⎟⎟ ,
u ( j ) = (1 /(2 2 n + 2 ))⎜⎜
⎝ j ⎠
impulsiva finita (FIR) descritti sotto:
(3.2)
0 ≤ j ≤ 2n+2 sono i filtri a risposta
Dove
(a)
(b)
Figura 3.5
Filtri FIR di decomposizione e ricostruzione: (a) basato sulle
linear splines e (b) basato sulle cubic splines
3.2.2
Multiscale pointwise product
Seguendo l’ analisi effettuata da Wang-wu et al. nel loro studio [2] viene osservato
come bordi e patterns permangono nella decomposizione attraverso scale differenti e
sono correlate spazialmente. In questa tesi, come nel loro lavoro viene pertanto
adottato il Multiscale Pointwise Product (MPP) per misurare la correlazione crossscale. Il MPP è definito come
K
p K ( n ) = ∏ W2 j f ( n )
(3.3)
j =1
22
Questo criterio è stato usato per il rilevamento e la localizzazione, per eliminare il
rumore e per il filtraggio di immagini di risonanza magnetica. In realtà, anche prima
dell’ utilizzo della wavelet transform, il MPP era stato usato per rafforzare i picchi
multiscala dei segnali e contemporaneamente ridurre il rumore, sfruttando la
correlazione multiscala dei segnali desiderati. Dato che i massimi di W2j f (n) dovuti
ai contorni nel segnale f (n) tendono a propagarsi attraverso scale differenti
contrariamente ai massimi dovuti al rumore, pK rafforza la risposta del segnale
piuttosto che del rumore. Il multiscale product ha quindi la capacità intrinseca di
sopprimere impulsi stretti ed isolati e al contempo preservare la risposta dei contorni
a scale differenti.
3.2.3
Enhancement algorithm
L’ algoritmo di miglioramento dell’ immagine può essere riassunto nei tre passi
seguenti:
1. Eseguire delle decomposizioni multiscala dell’ immagine usando le spline
wavelets ottenendo così una serie di decomposizioni lungo le direzioni
orizzontale, verticale e diagonale. L’ utilizzatore può decidere di modificare
queste componenti direzionali in modo selettivo attraverso dei coefficienti
passati in ingresso, al fine di evidenziare alcuni dettagli di informazione
piuttosto che altri.
2. Ordinare i valori di p Kd (n) d = 1, 2, 3 ad ogni pixel in ordine crescente. In
questo modo si nota quali componenti wavelet contribuiscono in misura
maggiore all’ immagine. L’ utilizzatore imposta una soglia μ per i diversi
gradi di enhancement. La fase di thresholding usando la formula seguente:
~
W2dj f (n) = λ( j )W2dj f (n)
= 0
if
if
p Kd (n) ≥ μ
p Kd (n) < μ
(3.4)
3. Dopo aver modificato le componenti wavelet nei passi precedenti viene
operata la trasformazione wavelet inversa usando la formula di ricostruzione
S 2 j −1 f = W21j f ⊗ ( g~ ( 2 ) , u ) ↑ 2 j −1 + W22j f ⊗ (u, g~ ( 2 ) ) ↑2 j −1 +
+ W23j f ⊗ ( g~ (1) , g~ (1) ) ↑2 j −1 + S 2 j f ⊗ (h, h) ↑ 2 j −1
(3.2)
4. Si esegue un histogram stretching dell’ immagine ottenuta.
Nella formula, la scelta di μ dipende dal livello del rumore nell’ immagine. Poiché le
immagini cromosomiche utilizzate in questo studio hanno poco rumore, il valore di μ
viene impostato molto piccolo. L’ utilizzatore ha l’ opzione di impostare differenti
valori di λ(j) per ottenere diversi gradi di miglioramento. Generalmente aumentando il
valore di λ(j) si ottiene un aumento più evidente del contrasto. Poiché le componenti
ad alta frequenza vengono di consueto evidenziate a risoluzioni a bassa scala,
conviene impostare un valore più alto di λ(j) a bassa scala.
23
(a)
(b)
Immagine di cellula in metafase in banda G, in (a) prima
ed in (b) dopo l’ algoritmo di image enhancement
proposto
3.3 Estrazione di parametri
L’ obiettivo finale di questa fase di un software di riconoscimento cromosomico è la
determinazione di un vettore di parametri, eguale in numero per ciascun cromosoma
che sia in grado di fornire una univoca descrizione di ciascuna classe. Tale vettore,
qualora passato in ingresso ad un classificatore opportuno dovrebbe essere in grado
di determinare la classe di appartenenza del cromosoma.
Si possono distinguere parametri cromosomici di quattro livelli: quelli di livello 1
vengono ottenuti da una misura diretta dell’ immagine cromosomica e sono ad
esempio l’ area cromosomica, la densità relativa (rapporto tra densità ottica totale e
area), e il perimetro del convex hull. Le caratteristiche di livello 2 richiedono l’ asse
come ad esempio la lunghezza cromosomica, quelle di livello 3 richiedono sia l’
informazione circa l’ asse , sia quella relativa alla polarità del cromosoma ed un
esempio è il density profile, infine i parametri di livello 4 necessitano ulteriormente
della conoscenza della posizione del centromero, ad esempio il banding pattern.
24
(a)
(b)
(c)
Figura 3.6
Parametri cromosomici di differenti livelli: in (a) un cromosoma, il suo convex hull
ed il rettangolo circoscritto di larghezza minima derivato dal convex hull, che viene
utilizzato per ruotare il cromosoma in posizione verticale (parametro di livello 1).In
(b) l’ asse mediano del cromosoma utile per trovare la lunghezza cromosomica
(livello2). In (c) viene evidenziato il banding pattern di un cromosoma ossia la media
dell’ intensità del cromosoma lungo linee perpendicolari al medial axis (livello 3).
Tra i lavori di letteratura utili alla scelta delle caratteristiche cromosomiche che si
prestano in maggior misura all’ identificazione, vi sono gli studi di Lerner et al. [6],
[7], [8], in cui viene chiaramente indicato come, benché le caratteristiche globali del
cromosoma come l’ area possano essere calcolate molto facilmente e siano adatte per
l’ implementazione, una combinazione di caratteristiche geometriche (lunghezza e
centromeric index) e di intensità (density profile) consenta di ottenere un miglior
discriminazione nella classificazione dei cinque tipi cromosomici considerati nel data
set dello studio. Il grafico riportato in Figura 3.7, mostra come un vettore di
parametri cromosomici costituito da lunghezza, CI e density profile consenta di
ottenere una percentuale di riconoscimenti corretti ottimale (compresa tra il 90% ed
il 100%), con l’ aggiunta di soltanto 5-10 parametri del density profile, mentre vettori
privi della lunghezza ma contenenti il centromeric index raggiungono percentuali
simili con l’ aggiunta di 18-20 caratteristiche del bandeggio ulteriori. Vettori infine
contenenti soltanto parametri del density profile necessitano di almeno 28-30
25
caratteristiche per raggiungere percentuali superiori al 90% (ricordando che vengono
considerati soltanto 5 classi cromosomiche) .
Figura 3.7
Studio effettuato da Lerner et al. [7] [8] per valutare l’ importanza del
numero di elementi (in ascissa) del vettore di parametri con cui viene
descritto un cromosoma per ottenere una classificazione corretta
(valore in percentuale in ordinata). La linea continua rappresenta il
vettore costituito dal solo density profile, la linea con” x” rappresenta il
vettore con centromeric index e density profile, la linea con “o” indica
l’ utilizzo di lunghezza, CI e DP.
Per questa ragione in questo studio è stato concentrato il lavoro sull’ ottenimento di
lunghezza, posizione del centromero (quindi centromeric index) e density profile.
3.3.1
Lunghezza
Questo parametro cromosomico può essere ottenuto semplicemente come una misura
in pixel della lunghezza dell’ asse cromosomico. Nel caso l’ asse venga calcolato
attraverso il computo della Medial Axis Transform (vedi par. 3.4) la lunghezza può
essere considerata quella del medial axis finché questo non si dirama, e quella della
mediana del triangolo costituito dalle due biforcazioni e dal bordo del cromosoma
dove iniziano le ramificazioni.
3.3.2 Density profile
Ogni campione del density profile viene definito da Moradi et al. [10] come la media
o l’ integrale dell’ intensità dei pixels che giacciono su di una linea perpendicolare
all’ approssimazione dell’ asse del cromosoma. A questo proposito,in corrispondenza
di ogni pixel che appartiene all’ approssimazione della curva centrale del cromosoma
26
viene considerata una retta e vengono mediati i livelli di grigio dei pixel che vi
appartengono, generando così un campione del DP-signal ossia un numero compreso
nel range 0-255. Nel caso di cromosomi a bande G (come quelli nel data set di
Moradi ) poiché usualmente è desiderabile attribuire i picchi del DP-signal alle
regioni più scure del cromosoma, ogni campione del DP-signal viene rimpiazzato dal
suo complemento a 255.
Figura 3.8
Density profile signal calcolato per i due cromosomi 17 e 18 appartenetnti al gruppo E.
3.3.3
Centromeric index
Nel lavoro di Moradi et al. [11] si considera che da un punto di vista morfologico il
centromero è il punto lungo l’ asse del cromosoma che possiede la minor larghezza.
In corrispondenza del centromero e distinguibile un restringimento lungo i bordi del
cromosoma. Per questo motivo è sufficiente ottenere anche in questo caso una
versione binaria dell’ immagine attraverso un processo di thresholding, con una
soglia che differisce da un’ immagine all’ altra. Una volta ottenuta la binarizzazione
vengono calcolati i vettori proiezione nelle direzioni orizzontale e verticale
sommando i valori dei pixel lungo le linee di proiezione. Quando sono stati ottenuti i
valori delle proiezioni orizzontali, per l’ intera lunghezza del cromosoma, il minimo
globale che si presenta nella parte centrale del cromosoma può venire considerato
come come locazione del centromero. Calcolando il medial axis, il punto di questo
che corrisponde alla posizione del centromero sul piano proiettivo può essere
approssimato come locazione del centromero sull’ asse e da questo è possibile
ottenere i vettori proiettivi. Questo metodo di calcolo del centromero necessita di
cromosomi che si trovi in posizione verticale e presenta tuttavia l’ inconveniente di
non riuscire a determinare una locazione corretta se vengono analizzati cromosomi
troppo piegati.
27
Figura 3.9
Il metodo di calcolo del centromero basato sul computo dei
vettori proiettivi.
La metodologia proposta da Lerner et al. [7] prevede di localizzare la posizione
centromerica ricercando la coppia più vicina di punti opposti del contorno lungo
linee perpendicolari all’ asse. Questo metodo da risultati migliori rispetto al
precedente tuttavia il problema critico di questo approccio si trova nel fatto che la
“testa” e la “coda” del cromosoma debbono essere tagliate, sarebbero infatti dei punti
di minimo assoluto, e la posizione del taglio è di difficile determinazione (non si
conosce a priori infatti se il cromosoma è metacentrico o acrocentrico).
Un tentativo di ovviare a questo inconveniente è stato proposto da Piper e Granum;la
soluzione adottata in questo lavoro [13] è quella di ottenere il profilo del cromosoma
inteso come distanza del bordo del cromosoma dall’ asse, e rifletterlo vicino alle
estremità. Se il profilo del cromosoma viene rappresentato come una funzione S di
una variabile intera che ha valori non nulli nel range [0-l] ed è zero altrove, dove l è
la lunghezza del profilo, allora il profilo è una funzione S’(i) = S(i) + S(k-i) + S(2l-ki) dove k < l è la lunghezza della sovrapposizione tra il profilo originale e le sue
riflessioni. S’ è non nullo nel range [k-l, 2l-k] . L’ effetto di questa trasformazione e
indurre un minimo locale artificiale nel profilo si un cromosoma acrocentrico,
distinto da uno metacentrico. Questo giace usualmente a k/2 dall’ estremo del profilo
originale e quindi k viene preso tale che k/2l è la lunghezza media CI per cromosomi
acrocentrici normali. Nel caso di cromosomi metacentrici il minimo globale del
profilo trasformato è quasi sempre ancora nel centromero, mentre nel caso di
cromosoma acrocentrico il minimo locale indotto nel profilo trasformato costituisce
anche il minimo globale (vedi ). Anche questo metodo tuttavia necessita di
cromosomi con bordi sufficientemente regolari, che consentano con una semplice
funzione di minimo, di determinare la posizione centromerica, condizione non
verificata in molte immagini di cromosomi in particolare in banda Q.
Figura 3.10
Trasformazione del shape profile di cromosomi con bandeggio G
attraverso la riflessione alle due estremità:Cromosoma metacentrico. La
posizione del centromero viene indicata dalla linea trasversale all’ asse.
28
Figura 3.11
Trasformazione del shape profile di cromosomi con
bandeggio G attraverso la riflessione alle due estremità:
caso di cromosoma acrocentrico.
Tutte le caratteristiche considerate maggiormente discriminanti ai fini della
classificazione richiedono quindi il computo dell’ asse del cromosoma. Ecco perché
il primo passo è la determinazione di un algoritmo che riesca ad ottenere una buona
curva centrale anche per cromosomi morfologicamente complessi.
3.4 Calcolo dell’ asse cromosomico
Gli algoritmi utilizzati per ottenere il medial axis possono essere suddivisi in due
grandi categorie dipendentemente dal tipo di informazione preservata: algoritmi di
tipo “skeleton” e algoritmi di “thinning”. Negli algoritmi skeleton ogni punto viene
specificato attraverso la sua distanza dal punto più vicino del contorno. Lo skeleton
viene a seguire definito come l’ insieme di punti la cui distanza dal più vicino punto
del contorno è localmente massima. Negli algoritmi di thinning invece i punti del
bordo vengono gradualmente eliminati se la loro rimozione non modifica la
topologia dell’ oggetto. I più utilizzati algoritmi di calcolo dell’ asse cromosomico si
fondano sul computo della Medial Axis Transform (MAT), una trasformazione molto
conveniente ai fini della rappresentazione di oggetti di forma allungata, dove la
larghezza degli oggetti contiene poca (se non nessuna) informazione utile, come ad
esempio nel riconoscimento di caratteri o proprio nell’ identificazione di pattern
cromosomici. La MAT di un oggetto preserva le proprietà topologiche dell’ oggetto e
riduce le esigenze di memoria e tempo di computazione.Dal punto di vista
geometrico il medial axis di un oggetto solido può essere definito come il luogo dei
centri dei dischi massimali inscritti nel dominio dato.
Nello studio di Moradi et al. [10] il computo del medial axis fornisce una curva che
passa attraverso il centro del cromosoma lungo la sua direzione longitudinale e si
divide in due parti nei pressi delle due estremità. Per ottenere il medial axis si
necessita di una versione binaria dell’ immagine, in cui l’ oggetto e lo sfondo sono
separati. La versione binaria può essere ottenuta attraverso un opportuno processo a
soglia. In seguito, nell’ algoritmo considerato nel lavoro di Moradi viene ricavata una
versione primaria dell’ asse basata sulla Euclidean distance transform ed in seguito
un processo di raffinamento consente di ottenere il medial axis. La lunghezza del
cromosoma viene quindi definita come la lunghezza della curva centrale che coincide
con il medial axis per la maggior parte, eccetto ai due estremi del cromosoma. In
questi punti il medial axis è diviso in parti e la curva centrale viene approssimata con
la mediana del triangolo formato da questi due rami e dal bordo del cromosoma.
29
Figura 3.12
Medial axis transform ottenuta secondo il procedimento descritto da Moradi: in(a) un
tipico cromosoma 17 dove è stata identificata la curva centrale; in (b) l’ immagine in
versione binaria dello stesso cromosoma; in (c) la distance transform dell’ immagine
binaria ed in (d) il medial axis del cromosoma dove sono ancora presenti le due
diramazioni.
L’ algoritmo determinato da Moradi in realtà non garantisce un asse cromosomico
connesso e della larghezza di un pixel. Altri metodi sono stati studiati per la
determinazione dell’ asse. Nello studio di Piper e Granum [13] vengono utilizzati due
differenti algoritmi per il calcolo dell’ asse cromosomico. Il primo, definito “poor
man’s skeleton”, è più rapido ma necessita di cromosomi non troppo piegati essendo
fondato sul calcolo del punto medio di ogni linea orizzontale attraverso il cromosoma
orientato verticalmente. Il secondo asse viene invece ricavato a partire dal
convenzionale algoritmo di skeletonization (operazione morfologica di dilatazione
seguita da una di erosione), seguito da un’ estensione delle estremità e da un’
operazione di smoothing (vedi Figura 3.14).
Figura 3.13
Stadi del computo dell’ asse a partire dall skeleton secondo Piper e
Granum; in (a) viene rappresentata l’ operazione di
skeletonization, in (b) la rappresentazione poligonale di (a); in (c)
l’ operazione di estensione delle estremità ed in (d) l’ asse dopo lo
smoothing realizzato attraverso una convoluzione con un filtro
passa-basso.
30
3.4.1 Algoritmo proposto
In questo studio l’ algoritmo utilizzato per il calcolo dell’ asse è stato ottenuto
prendendo le mosse da uno simile utilizzato già in ambito biomedico per effettuare il
tracking dei vasi sanguigni dell’ occhio umano. La computazione dell’ asse
cromosomico è basata sulla sua rappresentazione come immagine binaria. L’
immagine binaria viene ottenuta attraverso l’ applicazione di un processo a soglia all’
immagine in livelli di grigio. Le variazioni dei livelli di grigio tra tipi diversi di
cromosomi e le diverse condizioni di illuminazione impongono la scelta di una soglia
adattativa. MATLAB dispone di istruzioni che consentono, a partire da un’
immagine in toni di grigio, il calcolo di una soglia secondo il metodo di Otsu
(graythresh). Questa è la soglia che è stata utilizzata per ottenere la binarizzazione.
(a)
(b)
Figura 3.14
Operazione di thresholding con soglia di Otsu: in (a) l’ immagine da segmentare ed in (b) l’
effetto della binarizzazione a soglia.
Il software di calcolo dell’ asse implementato, cui in ingresso viene quindi passata la
versione binaria dell’ immagine, si compone di quattro step principali:
1) Determinazione del punto di partenza da cui in seguito viene generato
l’ asse completo. Il pixel da cui ha inizio l’ algoritmo può essere
scelto in modo abbastanza arbitrario, è sufficiente che si trovi in una
posizione sufficientemente “centrale” dell’ immagine cromosomica e
non eccessivamente in prossimità del bordo. Nell’ algoritmo
implementato viene scelto come iniziale il primo punto di una
versione dell’ immagine generata attraverso una operazione
morfologica di thinning, seguita da una di pulizia (pruning).
2) Il punto ottenuto viene usato come centro per la generazione di un
fascio di rette che vengono intersecate con l’ immagine binaria d’
ingresso ottenendone dei segmenti. Tra tutti i segmenti trovati
vengono messe in memoria le coordinate cartesiane degli estremi del
segmento più corto. Il punto medio del segmento minimo ottenuto
viene considerato il primo punto del medial axis del cromosoma.
Vengono determinati due punti successivi dell’ asse cromosomico
31
come punti medi di due segmenti paralleli al primo ma traslati di un
passo “sopra” e “sotto”.
3) I due nuovi punti trovati identificano col primo due vettori. Per
ottenere i punti ulteriori, questi vettori vengono prolungati di un passo
ottenendo due punti prova e i punti medi dei segmenti passanti per il
punto prova e perpendicolari ai vettori precedenti costituiscono due
nuovi punti dell’ asse. I vettori vengono aggiornati con gli ultimi punti
trovati. Questo passo viene ripetuto finché un punto non esce dal
bordo del cromosoma.
4) Vengono ricalcolati gli estremi dell’ asse trovato in precedenza e
considerati come prolungamento dei due precedenti.
Figura 3.15
Immagine cromosomica binaria con i punti della primaria
versione dell’ asse calcolato indicati con x.
5) Si esegue un’ operazione di smoothing dell’ asse con un’
interpolazione spline dei punti.
Figura 3.16
Immagine cromosomica con la versione finale dell’ asse.
32
3.4.2
Osservazioni
Osservazione 1. Al di là di una generica conformazione longitudinale, la struttura
naturale dei cromosomi è costituita da molte possibili configurazioni e forme, ecco
perché si rende necessario l’ ausilio di un efficace sistema di controllo, per evitare di
incorrere in anomale deviazioni dell’ asse dal suo regolare percorso. Nel caso in
esame il controllo considera distanza di due successivi punti dell’ asse ed angolo del
vettore che essi formano e, nel caso superino una soglia pre-impostata sulla base di
numerose osservazioni su immagini cromosomiche, termina la chiamata all’
algoritmo di calcolo dell’ asse precedente, l’ ultimo punto viene eliminato e si
riprende dal punto precedente con il calcolo del fascio di rette. Questo computo ha lo
svantaggio di contenere molti cicli annidati e quindi avere un peso notevole in
termini di velocità computazionale, per questo motivo si è cercato di accelerare il
processo attraverso una riduzione del numero di step adibito ad ogni ciclo, e l’
utilizzo dell’ algoritmo che procede per vettori successivi.
Osservazione 2. Nelle vicinanze di testa e coda del cromosoma si è reso necessario il
ricalcolo dei punti finali in quanto spesso il cromosoma ha in questi punti critici un
forma complessa (cromatidi in separazione, frammenti mancanti etc.).
Osservazione 3. L’ operazione terminale di approssimazione con una spline consente
anche di utilizzare un passo qualsiasi per i punti dell’ asse per ciascun cromosoma,
possibilità utile nel caso si desideri variare la dimensione del vettore con cui
descrivere un cromosoma.
Osservazione 4. Una volta ottenuto l’ asse del cromosoma è stato immediato
giungere alla lunghezza cromosomica come dimensione in pixel dell’ asse.
Osservazione 5. L’ asse così costruito, pur essendo specifico per forme
cromosomiche è risultato di ausilio anche nella realizzazione di un software di
segmentazione dove vengono considerati frammenti cromosomici di forme qualsiasi.
Figura 3.17
Immagine binaria di cromosoma con due touchings. L’ asse
del cromosoma principale non ha subito anomale
deviazioni.
33
3.5 Algoritmi che sfruttano il calcolo dell’ asse cromosomico
Una volta ottenuto l’ asse l’ estrazione degli ulteriori parametri di caratterizzazione
di un cromosoma si fonda su di esso. Gli algoritmi successivi cercano di estrarre le
altre due features di maggiore importanza, oltre la lunghezza cromosomica: banding
pattern e centromeric index.
3.5.1
Estrazione del banding pattern
Il banding pattern è stato ottenuto anche in questa tesi come la media delle intensità
dei pixel lungo linee perpendicolari al medial axis del cromosoma. Il metodo qui
utilizzato prevede i seguenti passi:
1) Ricavare l’ asse sfruttando l’ algoritmo precedente.
2) Per ogni punto dell’ asse determinare la tangente all’ interpolazione spline
della curva, calcolarne il vettore perpendicolare e determinare tutti i punti che
giacciono su questa retta compresi entro il cromosoma.
3) Mediare le intensità dei pixel determinati.
(a)
Density profile signal
180
160
140
DP-signal
120
100
80
60
40
20
0
0
10
20
30
40
Punti del medial axis
50
60
70
(b)
Figura 3.18
Density profile signal: in (a) è rappresentato il complemento di un’
immagine cromosomica in banda Q. In (b) il grafico del density profile, in
ascissa il numero di punti dell’ asse cromosomico, in ordinata la media dell’
intensità dei pixel lungo linee perpendicolari all’ asse.
34
Osservazione. Per ciascuno dei cromosomi presi in esame è stato trovato un vettore
di valori del Density profile signal, poi normalizzato in dimensione ed in valore, in
accordo con quanto svolto in molti dei lavori correlati visionati.
3.5.2
Calcolo della distanza del bordo cromosomico dall’ asse
L’ algoritmo di calcolo della distanza del contorno del cromosoma dalla sua curva
centrale è risultato essere utile non soltanto ai fini della determinazione della
posizione del centromero, ma anche in sede di segmentazione dell’ immagine.
Il calcolo di questa distanza consente infatti di ricostruire in un grafico una
rappresentazione rettificata della forma del cromosoma, utile a verificare la
somiglianza con un modello di cromosoma. Anche in questo algoritmo viene
utilizzata soltanto la versione binaria dell’ immagine cromosomica.
L’ algoritmo è stato costruito in questo modo:
1) Estrazione dell’ asse attraverso l’ algoritmo trattato in precedenza.
2) Dall’ interpolazione dell’ asse determinazione delle rette a essa
perpendicolari, computo delle intersezioni di queste con l’ immagine binaria
del cromosoma.
3) Calcolo della distanza tra gli estremi ottenuti in precedenza e l’ asse
cromosomico, positiva o negativa a seconda che i punti si trovino da una
parte o dall’ altra rispetto all’ asse cromosomico.
Chromosome border
15
Distanze del bordo
dall' asse in pixel
10
5
0
-5
-10
-15
0
10
20
30
40
50
Punti dell' asse cromosomico
60
Figura 3.19
Distanza dell’ asse cromosomicocalcolata per il cromosoma di
Figura 3.18 (a): l’ asse divide il cromosoma in due parti. Per
ciascuna viene calcolata la distanza del bordo dall’ asse. Il
grafico ottenuto mostra la sagoma del cromosoma rettificata.
35
70
Contour function
25
Funzione contorno
20
15
10
5
0
0
10
20
30
40
50
Punti dell' asse cromosomico
60
70
Figura 3.20
Funzione contorno: è ottenuta dalla somma tra le due
distanze dai bordi calcolate nel grafico precedente.
3.5.3 Calcolo del centromeric index
La posizione del centromero nei cromosomi delle immagini utilizzate è di difficile
localizzazione. Le condizioni di illuminazione delle immagini in possesso sono
scadenti e i bordi dei cromosomi sono poco definiti. Si è cercato di determinare il
centromeric index prendendo le mosse da alcuni algoritmi presenti in letteratura.
Come il precedente anche questo algoritmo necessita soltanto dell’ immagine in
versione binaria; i passi:
1) Computo di asse e distanza del bordo dallo stesso.
2) Determinazione della funzione contorno, ossia per ogni punto della curva
centrale, la somma tra le distanze del bordo cromosomico da una parte e dall’
altra dell’ asse.
3) Calcolo del minimo della funzione contorno dopo aver eliminato il 10% dei
punti dell’ asse agli estremi del cromosoma.
La posizione centromerica così ottenuta trascura tuttavia la possibilità che il
centromero sia situato alle estremità del cromosoma, evento che si verifica nel caso
di cromosomi acrocentrici e telocentrici
36
Contour function and centromere position
25
Contour function
Centromero
20
F u n z io n e c o n to rn o
15
10
5
0
0
10
20
30
Punti dell' asse
40
50
60
Figura 3.21
Determinazione della posizione del centromero come minimo della funzione
contorno nei punti “centrali” del cromosoma.
37
38
4 Risultati
4.1 Risultati dell’ image enhancement
Il miglioramento del contrasto dell’ immagine si può riscontrare visivamente
operando un confronto tra immagini prima e dopo l’ applicazione dell’ algoritmo. L’
immagine in fig mostra il risultato dell’ algoritmo su di un’ immagine di cellula
metafasica in banda Q. Inizialmente l’ immagine viene complementata per rendere il
bandeggio visivamente più chiaro (questa operazione non aggiunge informazione).
Viene poi applicato il metodo di enhancement, dove nella wavelet decomposition si è
giunti fino al livello 3 e i coefficienti λ(j) utilizzati sono 5, 2, 2 per la risoluzione in
bassa, media ed alta scala rispettivamente, mentre per μ è stato scelto 0.01.
(a)
(b)
39
Dare una quantificazione dei risultati del miglioramento di un’ immagine è
generalmente difficoltoso. Il contrasto viene spesso definito come la differenza di
luminosità media tra gli oggetti e lo sfondo. L’ effetto del miglioramento dell’
immagine alla fine viene giudicato sulla base dell’ accuratezza nella classificazione
dei cromosomi, per determinare se il miglioramento conseguito ha come effetto una
migliore cariotipizzazione e quindi una diagnosi semplificata. A questo proposito da
Wang et al. [2] viene utilizzato un classificatore bayesiano per effettuare la
classificazione dei cromosomi. Estraendo sempre le stesse caratteristiche utilizzando
la stessa procedura (che include lunghezza normalizzata, densità, centromeric index e
distribuzione di densità pesata), viene effettuato un confronto tra il metodo di
enhancement proposto e altri metodi, da cui risulta che il primo ottiene il più basso
error rate.
Figura 4.1
Curve dell’ error rate di classificazione come funzione del numero delle features
utilizzato. I metodi testati sono: adaptive contrast stretch (ACE), adaptive
contrast enhancement (ACE) e contrast gain transform (CGT). Il metodo
proposto è il multiscale contrast enhancement (MCE)
40
4.2 Risultati del calcolo dell’asse e degli algoritmi correlati
Risultati del DP signal
Risultati funzione contorno
41
Precisione determinazione del centromero
4.3 Determinazione distanza tra classi
Per avere una valutazione circa la capacità degli algoritmi di estrazione delle
caratteristiche implementati di determinare l’ appartenenza di un cromosoma ad una
determinata classe, si è ottenuto un vettore di parametri cromosomici di 65 elementi
per ciascun cromosoma. Il primo parametro è costituito dalla lunghezza del
cromosoma, gli ulteriori 64 sono valori del density profile. Le immagini da cui le
features sono state estratte sono 7 cariogrammi di cellula in metafase cui è stata
effettuata un’ operazione manuale di scontornatura per togliere materiale non
cromosomico presente a causa di una approssimativa segmentazione.
42
Utilizzando il vettore 65-dimensionale, innanzitutto è stato realizzato un plot delle
lunghezze dei cromosomi in relazione al loro tipo. Si vede chiaramente come sia
possibile differenziare abbastanza cromosomi, o restringere la possibile appartenenza
di un cromosoma con una determinata lunghezza ad 1-3 tipi.
Lunghezza normalizzata del cromosoma
Lunghezza cromosomica
Intervallo di confidenza
Lunghezza cromosomica normalizzata
Media
1
0.8
0.6
0.4
0.2
0
2
4
6
8
10
12
14
N° della classe
16
18
20
22
Figura 4.2
Sono stati poi calcolati il vettore medio per ogni tipo di cromosoma delle 64 features,
ottenendo un ‘prototipo’ per ogni classe. I seguito si è ricavata la correlazione di ogni
cromosoma coi prototipi, sperando che la correlazione massima si ottenesse col
prototipi della classe di appartenenza.
Pattern correlation
Dati ordinati per tipo di cromosoma
0.8
50
0.6
100
0.4
0.2
150
0
200
-0.2
-0.4
250
-0.6
5
10
15
Tipo di cromosoma
43
20
Figura 4.3
Il passo successivo è stato il calcolo di un vettore di features ottenuto prendendo il
segno della derivata del vettore di caratteristiche precedente. Il prototipo per ogni
classe è ottenuto calcolando il segno della derivata sui prototipi calcolati al punto
precedente
Correlazione della derivata
0.8
Dati per tipo di cromosoma
50
0.6
100
0.4
150
0.2
0
200
-0.2
250
5
10
15
Tipi di cromosoma
20
Figura 4.4
Si sono moltiplicate le correlazioni ai due passi precedenti ed in effetti sembra che i
massimi di questa misura si trovino in corrispondenza delle classi corrette.
Correlazione pattern X correlazione derivata
Dati ordinati per tipo di cromosoma
0.7
50
0.6
0.5
100
0.4
0.3
150
0.2
200
0.1
0
250
-0.1
5
10
15
Tipo di cromosoma
Figura 4.5
44
20
Osservazione 1: il cromosoma 22 ha correlazioni bassissime, ma probabilmente si
può classificare per esclusione.
45
Conclusione
46
Bibliografia
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
ISCN: an international system for human cytogenetic nomenclature, in
Cytogenetics and Cell Genetics
Y.-P. Wang, Q. Wu, K. R. Castleman and Z. Xiong, “Chromosome image
enhancement using multiscale differential operators” IEEE Trans. on
Medical Imaging, vol. 22, n.5 May 2003
Q. Wu, Y.-P. Wang, Z. Liu, T. Chen and K. Castleman, “The effect of
image enhancement on image pattern recognition”, Proc. Of the Second
Joint EMBS/BMES Conference, Houston, TX, USA, October 23-26, 2002.
Y.-P Wang, “Image representation using multiscale differential operators”,
IEEE Trans on Image Processing, vol. 8, n. 12, December 1999.
Q. Wu, Z. Liu, T. Chen, Z. Xiong, K. R. Castleman, “Subspace-based
prototyping and classification of chromosome images” IEEE Trans on
Image Processing, vol. 14, n. 9, September 2005.
B. Lerner, H. Guterman and I. Dinstein, “A classification-driven partially
occluded object segmentation (CPOOS) method with application to
chromosome analysis” IEEE Trans. on Signal Processing, vol. 46, n.10,
October 1998.
B. Lerner, H. Guterman, I. Dinstein and Y. Romem, “Medial axis
transform-based features and a neural network for human chromosome
classification”. Pattern Recognition, vol. 28, n.11, pp. 1673-1683 Great
Britain, 1995.
B. Lerner “Toward a completely automatic neural-network-based human
chromosome analysis”, IEEE Trans. on Systems, Man and Cybernetics –
Part B: Cybernetics, vol. 28, n. 4, August 1998.
Canadian Workload Measurement System, Minister of Supply and
Services Canada, Toronto, 88 (1986).
M. Moradi, S. K. Setarehdan, “New features for automatic classification of
human chromosomes: a feasibility study” Pattern Recognition Letters 27
(2006) 19-28.
M. Moradi, S. K. Setarehdan and S. R. Ghaffari, “Automatic landmark
detection on chromosome images for feature extraction purposes” Proc. Of
the 3rd International Simposium on Image and Signal Processing and
Analysis (2003).
M. Moradi, S. K. Setarehdan and S. R. Ghaffari, “Automatic locating the
centromere on human chromosome pictures” Proc. of the 16th IEEE
Symposium on Computer-Based Medical Systems (CBMS’03).
J. Piper and E. Granum, “On fully automatic feature measurement for
banded chromosome classification” .
E. Granum and M. Thomason “Inferred Markov network models for
classification of chromosomal band pattern structures” Cytometry 11 pp.
26-39 1990.
G. C. Charters and J. Graham, “Disentangling chromosome overlaps by
combining trainable shape models with classification evidence” IEEE
Trans. on signal Processing, vol. 50, n.8, August 2002.
G. C. Charters and J. Graham, “Trainable gray-level models for
disentangling overlapping chromosomes” Pattern Recognition 32 13351349 (1999).
47
[17] J. M. Cho, “Chromosome classification using backpropagation neural
networks” IEEE Engineering in Medicine and Biology January/February
2000.
[18] M. R. Ortiz-Posadas, “A mathematical model for classical chromosome
identification using the logical combinatory approach” IEEE EMBC pp.
1342-1345 (2003) .
[19] X. Shunren, X. Weidong and S. Yutang “Two intelligent algorithms
applied to automatic chromosome incision” IEEE ICASSP vol. 3, pp. 697700 (2003).
[20] P. S. Karvelis , D.I. Fotiadis, M. Syrrou and I. Georgiu, “Segmentation of
chromosome images based on a recursive watershed transform” IEEE The
3rd European Medical and Biological Engineering Conference, Prague
Czech Republic, November 20-25, 2005.
[21] L. V. Guimaraes, A. Schuck Jr and A. Elbern, “Chromosome classification
for kariotype composing applying shape representation on wavelet packet
transform” European Medical and Biological Engineering Conference.
[22] G. Agam and I. Dinstein “Geometric separation of partially overlapping
non-rigid objects applied to automatic chromosome classification” IEEE
Transactions on Pattern Analysis and machine intelligence, vol. 19, n. 11,
November 1997.
[23] N. Sweeney, R. L. Becker and B. Sweeney “A comparison for wavelet and
Fourier descriptors for a neural network classifier” IEEE/EMBS Oct 30Nov 2, 1997 Chicago, Illinois, USA.
[24] S. Delshadpour “ Reduced size multi-level neural network for human
chromosome classification” IEEE/EMBC 2003
48
Appendice
Nel caso unidimensionale (1-D) una wavelet differenziale che assomiglia al mexican
hat o Lapalcian of Gaussian (LoG) è definita come la derivata seconda di B-splines
di ordini differenti:
1 d 2 n+ 2
ψ =
β 2−1 ( x) = 2( β n (2 x + 1) − 2β n (2 x) + β n (2 x − 1))
2
4 dx
n
La wavelet transform è definita come:
W2 j f = f ⊗ψ 2nj
j= 1, 2, .... , J
Se introduciamo l’ operazione di smoothing
S 2 j f = f ⊗ β 2nj
j=1, 2, … , J
Si può dimostrare che le trasformazioni sopra sono le derivate seconde delle
operazioni di smoothing.
49
Una rappresentazione multiscala di un’ immagine può essere ottenuta dalla
convoluzione dell’ immagine con una funzione di smoothing bidimensionale (2-D)
definita come
Formula 3
a diverse scale.
Per l’ implementazione discreta, si approssimano le tre componenti delle derivate
direzionali, o decomposizione wavelet come:
formula
Abbreviazioni
DNA
AT
deossiribonucleic acid
adenosine trifosfate
50