Chromosome analysis by minicomputer

cmlPUTEH GllAPIIICS A~D nlAGE PHOCESSI~G (1973) 2, 402-416

Chromosome Analysis by Minicomputer MOHlo ONOE, Muao TAKAGI At~D KENICIII YUKIMATSU Institl/te of Industrial Science, Unit;ersity of Tokyo, Roppongi, Tokyo, japan

Communicated by A. Rosenfeld Receit;cdjuly 23,1973 The development of a small-scale chromosome analysis system is described. The system includes a minicomputer and an on-line microscope. It will pennit automatic stage control, focusing, and metaphase spread screening, and interactive processing of spreads. 1. INTRODUCTION

Chromosome analysis by computer has been one of the most interesting problems in biomedical image processing. Much excellent work [1-8] have been reported in this area utilizing relatively large image processing systems. These researches have shown that easy chromosomes can be analyzed automatically. In cases of more difficult chromosomes such as touching and overlapping chromosomes, they can be analyzed automatically to some extent, but interactive processing is necessary to check the results. Taking the str~>ng and rapidly increasing demand for automated medical inspection by image processing into consideration, it seems that the time has come to develop inexpensive biomedical image processing systems for practical applications. From this point of view, a small-size chromosome analysis system has been under development. The main aim of this project is not to develop a fully automated chromosome analysis system but to build up an inexpensive image processing system for saving human labor. The development of a system with the following features is now under consideration. A minicomputer should be the processor in this system so as to realize an inexpensive image processing system. This minicomputer will be used not only for chromosome analysis but also for other processing, such as control of an on-line microscope, acquisition of image data, and screening of metaphase spreads. An on-line microscope should be attached to this system to read image data directly without taking any pictures. This will save time and labor of technicians. Interactive processing is essential to tllis system, because it is extremely difficult for tile computer to analyze chromosomes automatically, especially for a minicomputer with severe hardware limitations. Assuming an operator to be available, it will be much better to fully utilize the excellent human pattern COPlTighl © 1973 by Aoademio Press. Ino. All rights of reproduction in any fonn res("l"\'ed.

402

CHRO~lOSOME

ANALYSIS DY

MINICO~lPUTEn

403

recognition capability through interactive processing. And this will lead to savings in cost. Automatic screening of good metaphase spreads that are worth analyzing is necessary [9]. It is time-consuming labor for a technician to find good spreads and to take their pictures. Even if automated chromosome analysis can be realFORTRAN 0 , except for image input and output operations. ized, it will only be a partial saving of human labor; much will still remain. Automatic focusing and automatic stage control of an on-line microscope are essential capabilities for perfonning the above-mentioned functions. An image processing system with these capabilities has been under consideration. As the first step, the feasibility of chromosome analysis by minicomputer has been checked, and it has been shown that in case of easy chromosomes, chromosome analysis is possible with this system. But this is only a small first step, and there still remain difficult problems in realizing the above described system. 2. THE SYSTEM

In this system, we have tried to develop a chromosome analysis system with a minicomputer and simple input and output equipment from the economical point of view and to supplement the insufficiency of hardware with software as far as possible. A mechanical scanner developed in our laboratory [10] is used as input equipment, and a typewriter and a storage-type CRT display (Tektronix 611) are used as output equipment. The minicomputer (HITAC-lO, 16 bits/word) has a core memory with 8 kilowords. At present, only paper tape is available as an external memory medium, so it takes a considerably long time to input and output data, because the processing is divided into several steps and the intennediate results are output on paper tapes. This disadvantage will be much reduced by using a magnetic disc and magnetic tapes which are now attached to the minicomputer. In this system, the programming language is FORTRAN 0 , except for image input and output operations. The block diagram of the system is shown in Fig. 1. It is extremely difficult to develop a fully automated chromosome analysis system with a minicomputer, but in this system the difficulty is solved by interactive processing with a storage-type CRT display and typewriter. STORAGE CRT

FIG.

1. Hardware system.

404

ONOE, TAKAGI, AND

YUKI~lATSU

3. FLOW OF PROCESSING

The processing is divided into five steps, and its flow is shown in Fig. 2. In· teractive processing is done in Steps II and IV, if necessary. The final result is displayed on the storage-type CRT display.

,

~

II

.: t

Scanning

Jl

/

Original

,

'It ~

~: \x

Segment

~) /

I

Axis

Variance

/

Rotation

L±l Profile

~

/Mm,eo"h ..

)

I'f. ,

~~ U:/!J~

Arm ratio

!Ii ~R ~W, .~'

Karyotype FIG. 2. Flow of processing.

The details of processing at each step are as follows:

Ste], 1. Digitization of the Image A picture of the chromosome spread is digitized into 350 X 350 picture elements with 128 gray levels by the scanncr, and the density of each picture element is convcrted into binary (black and white) data. The density slice level is detcnnined by the density histogram obtained by the scanning of one appropriate line as typical data. A string of picture elements on a scan line with gray levels above the ap-

CHHOMOSO~lE

ANALYSIS BY

MINICO~IPUTEH

405

SCANNING LINE

BUCK

I

WHITE SEG~ENT

7-7"" 'j"""\O ~.....-..-.

DATA POIlIlAT

lOS

114

0

121

126

o .•••

1

FIG.3. Digitization of image.

propriate threshold is named a "segment." Only the coordinates of the end points of cach segment are output on paper tape as the data (Fig. 3).

Step II. Separation of Chromosomes and Calculation of Center of Gravity and Principal Axis By detecting the overlap of each segment on each line, each connected object is separated individually, and it is determined to which object each segment belongs. The center of gravity and principal axis of each object are calculated. A certain threshold is applied to the area of each object to discriminate chromosomes from other objccts. Separation of objects. A segment number is assigned to each segment in the following manner, and the number shows to which object a segment belongs. When a new segment which has no overlapping segment in the previous line appears, a new segment number is assigned to the segment. If a segment has an ovcrlapping segment in the previous line, the same segment number as that on the previous line is assigned to that segment. A group of segments which have the same segment number is called a "segment group," and an object is defined as a group of several segment groups. The relation between an object and a group of segments is identified by an identification number. If more than two segments are merged into one segment as shown in Fig. 4, the segment number of the segment group on the left side is assigned as the identification number to all the merged segment groups. In the case of Fig. 4, the identification number assigned to the segment groups 1 and 3 is 1, and the identification number assigned to the segment groups 2 and 4 is 4.

406

mWE, TAKAGI, AND 11-

22--

1-- ) -

1-)-

1- )1)1)1- )-1--1--

1---

1--1-

1-1-1-

YUKI~IATSU

2--

224-24-- 2-

4-

4--4--4-44-4-

111-

1-

1-

1-

1-

1- 11 - 11-

4~­

5-

55-5--

4-44-44 'I 4--4-4-

FIG. 4. Segment numbers.

At the same time that this process is done, the summation of the coordinates of each picture element and the summation of their squared values are calculated for each object, and from these data, the location of the center of gravity and the direction of the principal axis of each object are calculated. Center of gravity. Let the coordinates of a sampled point with black level belonging to an object i be expressed as P(Xu,Yu) and the total number of the sampled points in the object i be expressed as PI. Then the coordinates of the center of gravity of the object i (XOhYO/) are obtained as X Oi =

{~Xu}/p, = E(Xu),

YOI =

{~ Yu}/P, = E(X u)· j=1

Direction of principal axis. As an object has a two-dimensional distribution of sampled points, the shape of the object is approximated by an ellipse. This ellipse is called as an eigenellipse and is calculated by the following statistical method: Let A = E(X U)2 - {E(Xu)}2, B=E(YiJ )2-E{(Yu)}2, . C = E(XijYij) - E(Xu)E(Yij).

A covariance matrix H is determined by A, B, and C: H=

(~ ~).

Let the eigenvalues of this matrix be AI> A2 (AI ~ A2) and the eigenvectors ell> (12. Then the ratio of the long and short axes of the eigenellipse corresponding to this object is AI/ 2: A~/2, and the direction of the eigenvector (11

CHI\O~IOSO~IE

FIG.

ANALYSIS BY

~IINICO~JPUTEn

407

5. Contour and principal axis of chromosome.

corresponding to Al coincides with the direction of the long axis (principal axis). The principal axes and contours of each object in the original picture are displayed on the screen of the storage-type CRT display as shown in Fig. 5. If the direction of the principal axis is not calculated correctly, this can be modified interactively from the keyboard of the teletype while looking at the screen of display. The principal axis can be rotated by an arbitrary angle in integer multiple steps of 15°.

Stel) III. Rotation of the Object (Iud Profile Calculation Each object is rotated so as to make its principal axis vertical, based on the data about the center of gravity and the direction of the principal axis acquired in Step II. The number of the sampled points on each line perpendicular to the principal axis and the average and variance of the coordinates of the sampled points perpendicular to the principal axis are calculated, and the pro-

408

ONOE, TAKAGI, AND YUKIMATSU y

---~~'-TI'-IIr--i'.-------

X

FIG. 6. Rotation of chromosome.

jected profile along the principal axis of the variance of the coordinates perpendicular to the principal axis is calculated. Rotation of the object. The object is rotated about its center of gravity so as to make its principal axis vertical. The relative coordinates (X', Y'), in which the center of gravity is the origin, of the object after rotation are calculated as X'

=

(X - Xo) cos () - (Y - Yo)sin (),

Y' = (X - Xo) sin () - (Y - Yo) cos (), where (X,Y) are the coordinates before rotation, and (Xo,Yo) are the coordinates of the center of gravity (Fig. 6). The calculated coordinates are real, but the value of the Y coordinate is made integer, counting as 1 fractions more than 0.5 inclusive and cutting away the rest, to make later processing easier. Projected profile ofvariance. Let the relative coordinates of a sampled point after rotation be (Xj,YJ. Also, let the average value of the X coordinates of the points (Xj,Yj ) on the line perpendicular to the principal axis (Y axis), where i = 1, 2, ... ni> and Yj is constant, be XOj (Fig. 7). The variance defined here is calculated as

All U"j for all j (-25 .s; Yj along the Y axis is made.

.s;

25) are calculated, and a profile of these values

Step IV. Determination of Centromere and Calculation of Arm Ratio and Arm Length From the profile of each object the arm ratios and arm lengths are calculated.

CHHOMOSO},lE ANALYSIS BY MINICO},IPUTEH Y

409

(Principal a.is)

00 00 00000 a

===o:J;[Q:O-=:=
Y=YI

000 0 00 00

x

o FIG. 7. Sampled points on a line perpendicular to the principal axis.

Determination of centromere. Several methods to determine the centromere, e.g., methods using contour coding, using the profile of the density distribution along the principal axis, have been proposed. In this system, a method using a profile is applied; the contour tracing method cannot be applied due to the shortage of core memory. The location of the centromere cannot be clearly determined by only using the profile of the sampled points as shown in Fig. 9(A), because in this system, the sampled points have only binary values (black and white) and have no weighting of density to help determine the centromere. Therefore, a new method has been introduced in which, using the profile of the variance described above, the location of the centromere is determined by finding the minimum point of this profile after smoothing of the histogram to remove noise. Figure 9(B) shows the profile of variance for the chromosome shown in Fig. 8(A). This method is based on the fact that a chromosome has nearly an "X" shape, and the location of the centromere can be determined quite well with this method. But in the case of acrocentric autosomes as in the D and G groups, the minimum point of the profile cannot be found. In such cases, the location of the centromere is determined using the derivative of the profile curve (curvature). With this method, the location of the centromere can be determined correctly for almost all chromosomes, but the location is sometimes determined

FIC. 8. Chromosomes.

410

ONOE, TAKAGI, AND YUKnlATSU

(A) POINT

(B)

VARIANCE

FIG. 9. Profiles.

L

--±-

--=L _~

.

~L

__ ~L

--.L

~_.

__1__

~

-lL_

~

2_ ---:L L _L 3-.-- --l ~

---±-_

~

-_1- _-=L _L

-L~_~_~_L~

,~ ~ ~ , .. 1-

~

.1..1_ --L- -~

-±L -L

l~ -..-'.-

~

~-L

FIG. 10. Profile and centromere.

--±---±.... 1/__

~_

.. I~_.

--=:±.:..:-

CIlHO~lOSO~lE

ANALYSIS BY

~IINICO~IPUTEn

411

erroneously. Such cases occur due to the following reasons. This system first looks for the dip in the pro6.le, and if no dip can he found, it then looks for the maximum change of the curvature of this profile. Therefore, a dip in the histogram due to noise is assumed to be the location of the centromere in the case of acrocentric autosomes, which have no real dip in the histogram. To remove such errors, the profiles and the locations of the centromere determined by the computer are displayed on the storage-type CRT display and if any error in the location of the centromere is found, this error can be corrected interactively. In Figure 10, vertical lines show the locations of centromeres, and more than two lines on the profile show modifications that were made in centromere location. Calculation of arm ratio and arm length. Once the location of the centromere is determined, the arm ratio and arm length can be calculated from the coordinates of the end points of the profile [7]. But in this method, much error will occur if a chromosome bends considerably. Therefore, a simple method, taking bends of the chromosome into consideration, has been found. Since the number of the sampled points, the averaged value of the coordinates of the sampled points, and the histogram of variance have already been calculated, the following method is applied to the calculation of the arm ratio and arm length. The locations of arm ends are determined at the points shifted horizontally by the averaged value of the X coordinates from the principal axis at the points of maximum and minimum Y coordinates, as shown in Fig. 11. Also, the arm length is not calculated by the length of a stmight line between an arm end and centromere, but is calculated by the lengths of two lines as shown in Fig. 12. Shifting the midpoints (Pa, P4 ) by an amount proportional to the standard deviation of the X coordinates, to points (P~, P~), the arm lengths are calculated as the length from PI to Pc via P~ and the length from P2 to Pc via P~. By this method, the arm length clm be measured approximately, while also taking the shape of chromosome into account. Table I shows a comparison of the measured values of total length and arm ratio on the chromosomes shown in Fig. 8 by several methods. In the case of the method described here, the proportional constant may vary depending on

I Yj

'000. ~

ooo~

ooq

~ogg

\ 00 0\ \0000\

\

I

\

\

\

"

AIII~~ x"" 00t?

r~ I. 1

I

\.

,..., /0000/

I

" \J/ / \ "" I ~Bwr~Bi "I

"

I

6

I

FIG.

1

\

11. Ann end.

/

/

//

412

ONOE, TAKAGI, AND YUKIMATSU

FIG. 12. Mcasuremcnt of ann Icngth.

TABLE I CO~IPARISO="l OF AIL\! LEl\'GTIl r-.IEASURE~IE="IT ~IET"ODS

A

Length Method used by JPL [7J Approximation by lines 0-1.5: proportional constant

59.0

o

59.4

0.3 0.7 4.0 1.5

60.0 62.4

65.3 71.8

Human measurement

59.0

or

A Proportional constant Method used hy JPL [7] Approximation by lines

Human measurement

Length 59.0

o

0.3 0.7 1.0 1.5

59.4

60.0 62.4

65.3 71.8

59.0

TABLE II CLASSIFICATIO:"< BY MI;o.;ICO~IPUTER

No.

Length

Ratio

Group

38 40

53.2 50.3

0.513 0.500

A A

3 19

48.8 46.7

0.608 0.582

A A

35 31

44.5 39.4

0.523 0.56.3

A A

30 6 9 5

40.4 39.7 38.9 36.7

0.721 0.691 0.756 0.658

B B B B

13 8 2 42

39.5 35.0 3·1.6 33.6

0.614 0.635 0.614 0.571

27 24 21 10 1 39 45 28

33.1 30.5 30.5 30.4 30.1 29.5 28.6 27.7

0.5·19 0.661 0.597 0.629 0.674 0.586 0.606 0.621

C C C

34 41 4 32

27.6 27.3 25.7 23.8

0.623 0.591 0.557 0.610

C C C C

16 20 26 18 25 43

23.5 22.5 22.4 21.6 21.5 20.4

0.777 0.813 0.856 0.711 0.803 0.84·1

D D D D D D

12 7

20.6 20.1

0.552 0.582

E E

14 15 37 44

18.9 18.2 17.9 16.6

0.762 0.690 0.622 0.805

E E E E

22 11 29 17

17.8 17.7 16.9 15.8

0.542 0.527 0.560 0.585

F F F F

23 33 36 46

13.6 12.5 12.4 11.5

0.608 0.737 0.741 0.809

G G G G

C C C C C C C C C

414

OXOE, TAKAGI, AND YUKnlATSU

•

.1

Zf" -F22

I,.

FIG.

... $6

2J

"IJ

-6-

" II

X

t

13. Karyotyping of chromosomes.

the shapes of the chromosomes and the original image quality, but a range between 0.3 and 0.4 wiII be preferable. Step V. Karyotyping Using the infomlation about the amI ratio and ann length of each chromosome obtained in Step IV, karyotyping is perfonned in Step V. First, the chromosomes are arranged in order of total length and then partially rearranged taking ann ratio into account and classified into seven groups. At present, chromosomes are classified into 12 groups (3 groups in A, 3 groups in C, 2 groups in E, and 1 group in each ofB, D, F, and G). Each group has its own threshold value of ann ratio. Figure 13 shows the original chromosome picture and its manual karyotyping. The results of classification of the same picture by this system is shown in Table II and Fig. 14. Figure 14 shows the results of Table II

CHRO~IOSO~IE

FIG.

ANALYSIS BY

MINICO~IPUTER

415

14. Karyotyping by minicomputer.

displayed on the storage-type CRT display. The results coincide well with these obtained by manual karyotyping. 4. CONCLUSIONS

It has been shown that in the case of images without overlapping chromosomes and with comparatively good quality, the chromosomes can be classified correctly into seven groups. At present, this system has only small merit, because it needs to input and output data by paper tape several times, and this takes a long time. But with a magnetic disc and magnetic tapes, processing time will be much shortened. Actually, it is rather difficult to find images such as those processed in this system without overlapping chromosomes and with good quality; most chromosome images have touching and overlapping chromosomes. Therefore, interactive processing as the next step should be developed in this system to separate overlapping chromosomes from each other. An on-line microscope now under development will be attached to this system. The system should be developed to be an integrated, highly economical, and practical system.

416

O;,\,OE, TAKAGI, AND YUKIMATSU ACKNOWLEDGMENT

The authors express their gratitude for the encouragement and advice received from Prof. A. Koizumi and Dr. S. Kaihara of the Department of Medicine, University of Tokyo. REFERENCES 1.

2. 3. 4. 5. 6. 7. 8. 9. 10.

J. W. BUTLER, ~r. K. BUTLER, A:"D A.

STROUD, Automatic classification of chromosomes-III, in Data Acquisition and Processing in Biology and Medicine Vol. 5, (E. Enslein, Ed.), Pergamon Press, New York, 1968, pp. 21-38. ]. HILDITCH A:"D D. RUTOVITZ, Chromosome recognition, Ann. Nr Acad. Sci. 157, 1969, 339-36·1. G. GALLUS AND P. W. NEUHATH, Improved computer chromosome analysis incorporating preprocessing and boundary analysis, Phys. Aled. Bioi. 15, No.3, 1970,435-445. F. RUDDLE, S. S~lITH, n. S. LEDLEY, AND ~1. BELSO:", Replication.precision study of manual and automatic chromosome analysis, Ann. Nr Acad. Sci. 157, 1969,401-423. II. S. FHEY, An interactive computer program for chromosome analysis, Comput. Biomed. Res. 2, 1969,247-290. M. L. ME:-;DELSOIIX, D. A. IIU:"CERFOHD, B. II. MAYALL, B. PERRY, T. CO:"WAY, AND J. M. PnEWIIT, Computer-oriented analysis of human chromosomes - II, Integrated optical density as a singl~ parameter for karyotype analysis, ArlII. NY Acad. Sci. 157, 1969, 377-391. K. R. CASTLE~IA.-':, R. ]. WALL, II. J. FRIEDE:", J.]. RuSSO, A.'o;D J. T. E:-;CLlSII, Method and means for chromosome karyolyping, JPL New Technology Report, Invention Report No. 30-11852, Case No. 2500, 1969. S. P. STO:-;E, J. L. LIITLEPAGE, AlI:D B. R. CLEGG, Second Report on tllC Chromosome Scanning Program at the Lawrence Radiation Laboratory, UCRL-71493, 1969. N. WALD A/I;D K. PHESTO="', JR., Automatic screening of metaphase spreads for chromosome analysis, in Image Processing in Biological Science (D. M. Ramsey. Ed.), University of California Press. 1968, pp. 9-34. M. O:"OE, ~f. TAKAGI, T. MASU~IOTO, AND N. HA~IA:-;O, Simple input/output equipment for digital processing of pictorial infornlation, Seisankenkyu 24, No.4, 1972, 127-134.