Bayes' theorem applied to chromosome classification

COMPUTERS

AND

Bayes’

BIOMEDICAL

Theorem

RESEARCH

4, 561-570 (1971)

Applied

DAVID A. Low, WILLIAM

to Chromosome

Classification

D. SELLES,AND PETERW. NEURATH

Physics Division, Department of Therapeutic Radiology, Tufts-New England Medical Center Hospitals, Boston, Mass. 02111

ReceivedFebruary4,197l As part of aninteractive,computer-aided imageprocessing systemdeveloped at TuftsNewEnglandMedicalCenter,a computerprogramto classifyhumanchromosomes into ten major groupshasrecently beendevised.The programappliesBayes’theoremto rankingsto determinethe mostprobablegroup into whicha chromosome will fall. A methodof calculatingthe requireddistributionsfrom a relativelysmallsampleis given. The programthenproceedswith additionalconstraintsto constructa karyotype. Four parametersderivedfrom the measurements are used:normalizedlength,normalized area,andcentromericindicesby lengthandarea.On a sampleof 13800chromosomes from 300“normal” spreads, therewere1.3% discrepancies between manuallydetermined andcomputedgroupassignments. Two hundredandnineteenout of the 300spreads had no discrepancies. Pair assignments within the ten major groupswerealsocomputed. Agreementwith manualassignments wasrelatively low but muchbetter than chance, asshownby quantitativeestimates. The computedmethodappearsto he similarto humanperformance,quantitatively justifying karyotyping into the ten groupsbut pointingout the needto developnew methodsfor identifyingpairsotherthan 1,2,3, and16. The adaptationsneededto makeBayes’theoremapplicableto this problemmay extendits usefulness to otherbiomedicalareas. I. OBJECTIVE One objective of any system for automatic chromosome analysis is to make largescale screening both feasible and economical. This involves, in addition to the optimization of measurementhardware and software, a computer program for classifying individual chromosome spreads.A “good” classification schemewill minimize the number of samplesrequired from an individual and will also aid in the measurement procedure, e.g., by presenting “doubtful” chromosomes for remeasurement. This paper describessuch a classification method currently in use at the Tufts-New England Medical Center. The method is based on Bayes’ theorem, which basically provides a straightforward approach to the classification problem. To make it work we used the rankings of four computer measuredparameters as variables and devised a method to obtain a large designset from which to calculate the a priori distribution coefficients. Results on 13 800 chromosomes from normal spreads show discrepancies at an 561

562

LOW,SELLES,

AND NEURATH

average rate of less than one chromosome per spread, indicating the superiority of this method over that of stepwise discriminant analysis (I). II. MATERIALS

AND METHODS OF MEASUREMENT

All spreads were prepared by a modification of the Moorehead procedure (4. Thirty-five millimeter photographs of nonoverlapping, unbroken, complete human chromosome spreads were subjected to analysis using an on-line programmable flying spot scanner and interactive display monitor (3). Chromosomes were located and separated under program control (4, 5) and on each a number of parameters were obtained with methods utilizing boundary string analysis (6). The measured parameters represent total chromosome length (L) and area (A), and centromeric indices (ratio of small arm measurements to those of the whole chromosome) by length (CZ,) and area (CZ,). A preliminary classification scheme by a method described previously (4) was used recursively with this measurement program in about 20% of the spreads reported here. This method categorized chromosomes on the basis of rankings on various combinations of L and Cl, or A and CI,. Those chromosomes whose “L” and “A” classifications did not agree were subjected to immediate remeasurement, without any indication to the operator as to the nature of the disagreement. This served as a check on chromosomes which might have been improperly measured the first time and helped avoid the collection of poor data. The remaining 80% of the spreads were measured with a program utilizing the present classification scheme recursively in a similar manner. Finally, the measurements on all spreads were subjected to the present classification scheme. 111. GROUP CLASSIFICATKON METHODS A general method applicable to both normal and abnormal spreads has been devised. The preliminary method, while yielding low discrepancy rates between manual and automatic classification was not sufficiently general and could only handle normal spreads. A. Useof Ranking

In developing this general method, it was decided to retain rankings as a classification criterion, in emulation of the cytogeneticist who ranks when determining his karyotype. A disadvantage of rankings is that they can disguise differences in an unevenly spaced distribution. Despite this, and the fact that available statistical methods are relatively primitive, rankings have several advantages for our purpose : (1) They provide a discrete event space which speeds computation. (2) Adjustment of normalization factors is not necessary for spreads of different sex or chromosome number.

BAYES’

THEOREM

AND

563

CHROMOSOMES

(3) Group separations are as good as those obtained from normalized measurements. Although not strictly applicable to nonnormal distributions such as rankings, Mahalonobis D* separations (7) can be used to qualitatively compare ranking vectors with normalized measurements. Table 1 lists such separations for several groups which require differentiation. In no case is the loss of information due to ranking vs. measurements serious; and for the smaller chromosomes, which are more difficult to classify, the ranking separations compare favorably. TABLE

MAHALONOBISSEPARATIONOF

Groups Differentiated

----__. N

Al-A2 Al-A3 B-C(X) E6-F E-D E-G E&E

14.8 6.0 14.6 5.0 17.8 1.4

L

.I

R .9 11.0

4.6 15.4 6.0 22.2 1.7

-

1

MANUALLY

~._ N

AWGNED GROUPS A ___--_ N

cl, R

22.2 .3 3.0 1.3 10.5

14.3 .2 3.1 1.4 11.9

.l

.l

1.7

9.3

.2 12.3 6.0 13.8 4.2 24.1 1.5

R .3 9.0 4.1 15.3 5.1 22.1 1.7

____---. N 53.5 2.2 3.1 1.4 18.3 1.7 11.7

crA

R 25.0 2.0 3.2 1.4 11.0 1.6

15.1

“100 Spreads: 50 male, 50 Female. N = Direct Measurements. R = Rankings of Measurements.

B. Bayesian Method

With the Bayesian approach (8) we adopted, we computed the probability P( M lr,, . . ., rr) of a chromosome belonging to a group M on the basis of its rankings rl,. . ., rt. The incidence ratios used to compute these probabilities are derived from a design set of N spreads, involving the number of times C1, (1 < i < t) a chromosome of rank ri belonged to group M. In particular, if the mean complement of group m is S,. g**ti P(MIr,,...,rt)=

numerator +

46-S,,, 46

&‘” N-Ci ’ * iG N(46 - S,)

(1)

This approach discriminates only between two groups: “M” and “Not M” and because our parameters are not independent, its use is only justified by the results. For any particular chromosome, the sum of its probabilities for belonging to any group need not be one, and in fact can be zero. Thus a chromosome need not fit well into any group, allowing us to classify a chromosome as “abnormal.”

564

LOW,

SELLES,

AND

NEURATH

C. Most Probable ClassiJication In an experimental spread about to be classified, the P’s are computed for each chromosome j, for each group m, to yield Pj,,,. Thus classification of 46 chromosomes into ten groups requires calculation of 460 probabilities. For chromosome j, that group M for which Pjm is a maximum, is called its “most probable” classification. A karyotype for any one spread made up of chromosomes so classified has had no constraints applied as to maximum group size (S,). D. Additional Non-Bayesian Constraints Such constraints can be used to produce an alternative karyotype which we have called the computed karyotype. For each group m, up to S, chromosomes are chosen by ranking the Pjm for that group. The number chosen will be less than S,,, if all remaining P,im are zero before S, is reached. The groups are filled in a predetermined manner, no chromosome being considered if it has already been placed in another group. If, at the conclusion of the process, a chromosome has not been assigned to a group, it is given the designation “?“. This can result from all Pjm for that chromosome being zero, or from all Pj,, not equal to zero ranking less than S,,,among unclassified chromosomes. A sufficiently unusual chromosome may thus remain unclassified. To put these procedures into use, a way must be found to determine the coefficients Ci in Eq. (1) and to determine the optimal order of redistribution according to the secondary constraints. E. Creation of Data Base In order for Eq. (1) to approximate population probabilities, N must be very large. Since a sufficiently large amount of data was not available, a “design” set of 45 spreads (3 females: 15 spreads each) was used to create a much larger “artificial” design set of 1000 spreads. Later in the experiment, 150 measured spreads (75 male, 75 female) were treated in the same manner to yield an artificial design set of 5000 spreads. It was found that this did not change or improve the classification. Indeed, more than 90% of the chromosomes misclassified by the 45-spread base were also misclassified by the 150-spread base. Therefore, the original data base was retained. Each chromosome in the original design set was classified manually into one of 23 “pairs.” Then 1000 artificial spreads (500 male, 500 female) were created by randomly selecting 46 chromosomes for each: two from each of the 23 sets of 45 pairs, except that for males only one was chosen from pair 23, the “X,” and three from pair 21. These artificial spreads were then ranked on each of the four parameters, and the results used to determine the Ci in the Eq. (1). As an example of the results obtained, a chromosome ranking first in length has probabilities of 38.6 %, 30.0 %, 18.1 %, 13.09 %, .27 7; of belonging to pairs 1, 2, 3,4, or 5, and the C group, respectively, and zero for all others.

BAYES'THEOREMAND

CHROMOSOMES

565

F. Choice qf’ Groups and Order of Assignments

Ten groups provided the base for classification: Al, A2, A3, E, C(6-12, X), D, E16, E( 17-18) F, G(21-22, Y). In order to decrease parameter dependence within the C group, it was divided into two groups consisting of chromosome pairs 6, 7, X, 8, and 9, 10, 1 1, 12. Both groups were filled separately, the size restrictions being only that their sum must be less than or equal to the complement for the entire C category. They were then merged again at the end of the procedure. The order in which the groups are filled is Al, A2, A3, B, F, E16, G, D, E, C. G. Sex Determination

The patient phenotype was used to determine whether to expect a Y(G) or X(C) chromosome. If this information was incorrect, it was readily evident in the resulting classifications. H. Rankingsfor Non-Modal Spreads

Rankings

for spreads with chromosome count (CC) # 46 were normalized to off the r * to the nearest integer provides duplicate or missing rankings around 23 depending on whether there are more or fewer than 46 chromosomes, respectively. Since rank 23 is usually occupied by a member of the C group, resultant problems are few; specifically that is the case because the C group consists of all those left over after all other groups have been chosen, unless they have no probability of being a C. r * = r *(46/CC). Rounding

I. Computational Features

The classification subroutine is written in IBM 360 assembly language. It takes about 5 seconds on an IBM 360/30 to compute the two classification vectors for 46 chromosomes from the original Ci. The subroutine and Ci matrix occupy approximately 12 000 bytes of core storage. Listings are obtainable from the authors. IV. RESULTS OF GROUP CLASSIFICATION A. Accuracy

qf Analysis for Normal Spreads

Three hundred normal spreads, including 45 “design” set spreads whose discrepancy rates did not differ from the remainder of the set, were classified with the “computed” (M) and “most probable” (P) methods and compared with previous manual classifications. Discrepancy rates are 1.3 % for the A4 and 1.8% for the P classification. The discrepancy matrix for the M classification is shown in Table 2. The M classification resulted in 219 (73%) zero discrepancy spreads; the P classification resulted in 150 (50 ‘A). A stepwise linear discriminant analysis program developed by Neurath and Enslein (1) was applied to a representative sample of the same 300 spreads. It used various linear and nonlinear combinations of the same parameters and resulted in a discrepancy rate for most probable assignment of 4.1 “/L

566

LOW, SELLES,AND NEURATH TABLE 2 300 NORMAL

Manual assignments Al A2 A3 B c D E6 E F G

SPREADS

(150

150 FEMALE)

MALE,

M

CLASSIFICATION

DISCREPANCIES

Computer assignments Al

A2

A3

Ill 4 4

1 /ii 1

4

B

C

D

Eb

E

F

G

3 I 2 4

/ii I

llll 14

?

9 llil //ii 2 9

//ii 31 9

zo //ii 4 4

10 ,;,i I

5 l/l/

:‘:,

8 5 7

1.3 .8 1.2

13

1.1

5 8

26 9 40

.6 .5 6.7

I :

43 I3 8

4.1 1.2 .7 1.3

6

5

Total

179

on the design set. This was better than on the data originally reported, reflecting the better measurements of the current data, but still much inferior to the method here reported. A detailed examination of those chromosomes which produced discrepancies between the manual and M classifications showed that, in general, the M classification reflected the recorded measurements, which either did not quite agree with what the eye “saw” in a photograph of the spread or else gave rise to ambiguous chromosomes, which the program quite logically placed differently than the technician in a number of cases. Changing the logic only produced more errors than were corrected. B. Abnormal Spreads Because of the different criteria used to determine the M and P classifications, discrepancies between the two classifications can be used to spot and identify abnormalities. One hundred fifteen abnormal spreads consisting of incorrect sexing, deletions, XXY, X0, D trisomies, E trisomies and G trisomies have been subjected to analysis. Although the overall success rate is high (70x), for several kinds of abnormalities the classifications err frequently and additional work is required to determine whether more reliable estimates of classifying missing or extra chromosomes are possible. V. COMPUTED PAIRASSIGNMENTS A. Method Having met with considerable success in correctly assigning chromosomes into ten groups by means of a program, we extended the procedure to assign pairs, at least

BAYES’

THEOREM

AND

CHROMOSOMES

567

for “normal” karyotypes. Qualitatively, such pair assignments are known to be highly arbitrary, even when made by a trained technician; but we lacked any quantitative estimate of how good, or bad, the agreement between assignments made manually and by the computer would be. In addition, the karyogram display produced as part of the computer output requires us to order the chromosomes within each group in some way: we preferred to do it in a manner that would approximate the technician’s procedure more closely than by simply sorting according to decreasing size. These motivations led to a procedure which is applied to all groups containing the nominal number of chromosomes. The decision rules are based on the distribution of mean values for manually assigned pairs in a two-dimensional space defined by total arm length (L) and centromeric index by area (Cf,, referred to below as CX). The criteria were arrived at in collaboration with a cytogeneticist and are consistent with the Denver classification. Group A. Pairs are determined by the existing group classification program. Group B. The two longest are assigned to Pair 4; the two shortest to Pair 5. Group C. Of the five (males) or six (females) longest, one (or two) with the maximum CXjL are assigned to Pair X, of the four remaining long C’s, the two with greatest CX.L are assigned to Pair 6 and the others to Pair 7. Of the remaining shorter five pairs of the C group, the two with the largest CX (excluding the longest pair) are assigned to Pair 11. Of the four pairs remaining, the two with the smallest CX are assigned to Pair 10. Of the remaining three pairs, the two smallest are assigned to Pair 12. Of the remaining two pairs, the two having greatest CX/L are assigned to Pair 9, and the last two to Pair 8. Group D. The two longest are assigned to Pair 13, the two next longest to Pair 14, and the two smallest to Pair 15. Group E. Pair 16 is assigned by the existing group classification program. Of the four chromosomes remaining, the two with greatest CX. L are assigned to Pair 17: the other two to Pair 18. Group F. The two longest are assigned to Pair 19 and the two shortest to Pair 20. Group G. For males, the longest is assigned to Y. Of the four remaining, the two having the greatest CX.L are assigned to Pair 21; the other two to Pair 22. B. Results Agreement with our manual pair assignments was evaluated as follows: Each computed group in 300 spreads was considered separately. In order to allow a valid comparison between programmed and manual pair assignments without introducing discrepancies in groups, those groups which did not contain the nominal number of chromosomes, or which contained chromosomes from other groups (according to their manual assignments) were excluded from the analysis. The agreement in assigning chromosomes to pairs within each group by the above procedure (excluding the single pair groups Al, A2, A3 and El6) was calculated and compared with the agreement expected when pairs are assigned randomly

568

LOW,

SELLES,

AN5

NEURATH

within the group. These data are shown in Table 3. The highest ratio of improvement over chance was 4.9 for the Cl I pair, the next highest 4.7 for the Cl0 pair, with other TABLE

AGREEMENT RATES

Pair 84 B5 C6 Cl C8 c9 Cl0

Cl1 Cl2 x 013

014 D15 El7

El8 F19

F20 G21 G22 Y

FOR PAIR

No. of subjects

No. of chromosomes

288 288 285 285 285 285 285 285 285 285 290 290 290 256 256 283 283 295 295 295

576 576 570 570 570 570 570 570 570 423 580 580 580 512 512 566 566 590 590 145

-

ASTIGNMENTS

‘A agreement 74.0 74.0 42.8 28.1 38.8 32.5 61.2 63.0 43.6 32.4 65.6 47.9 65.7 85.7 85.7 60.1 60.1 53.4 55.3 44.9

3 MADE

WITHIN

KNOWN

GROUPS

% agreement expected with random assignment 50.0 50.0 12.9 12.9 12.9 12.9 12.9 12.9 12.9 9.5 33.3 33.3 33.3 50.0 50.0 50.0 50.0 45.1 45.1 20.0

Ratio improvement

of

over chance 1.5 1.5 3.3 2.2 3.0 2.5 4.7 4.9 3.4 3.4 2.0 1.4 2.0 1.7 1.7 1.2 1.2 1.2 1.2 2.2

improvement rates ranging downward to a low of 1.2 for the G21 pair. The higher rates obtain for chromosomes which are most easily distinguished from others in their group by human inspection; thus the Cl0 is well known to be the most acrocentric of the C’s, the Cl1 the most metacentric, and the P19 and F20 to be practically indistinguishable from each other, as are the G21 and G22. Pair assignment agreement rates were also compared by groups, as shown in Table 4. Among those groups having more than one pair of chromosomes, the assignment procedure was most successful in group E (El7 and El 8), where agreement with the cytogeneticist was obtained 86% of the time. Perhaps the most interesting group, because of its size, is the C group. Although its agreement rate of 43 % was the lowest, its improvement ratio over chance assignment, at 3.4 was the highest of any of the groups. An attempt was made to assign the chromosomes of group Cinto two classes rather than eight pairs with 6,7, and Xin one class and 8-12 in the other. Assignment to these classes resulted in agreement rates of 81.6% and 89.7%, respectively, representing improvement ratios of 2.3 and 1.4. Though all of these rates of pair prediction

BAYES'THEOREM

PAIR ASSIGNMENT

Group

No. of subjects

Al A? A3 B c D El6 E F G

300 300 300 288 285 290 300 256 283 295

No. of chromosomes 600 600 600 1152 4413 1740 600 1024 1132 1325

569

AND CHROMOSOMES

TABLE 4 AGREEMENT RATES BY GROUP

‘A agreement 99.0 98.2 97.7 74.0 43.1 59.7 92.7 85.7 60.1 53.3

% Agreement expected with random assignment

50.0 12.6 33.3 50.0 50.0 42.4

Ratio of improvement over chance 1.5 3.4 I .8 1.7 1.2 1.3

within groups appear significant, they do not equal our predictive ability for pairs A 1, A2, A3 and E16, which are considered groups by the first part of the program. As shown in Table 5, the agreement rates for these pairs are 99 “,/,, 98 %, 98 % and 93 %, respectively. Having thus established reliability estimates for our computed pair assignments independently of their group assignments, we have included the pair assignment procedure in our analytical programs. Chromosomes having probable or undetermined group assignments are placed into unfilled groups for display in the karyogram, and when this cannot be done reasonably, no attempt is made to assign pairs in the affected groups. VI. CONCLUSIONS Our conclusion at this point is that the computer classification works well, perhaps nearly as well as the manual methods. Its very good agreement for groups with the manual method confirms our objective basis for the latter. It also demonstrates the usefulness of a method for obtaining distributions for the application of Bayes’ theorem from a limited measured design set and the advantage of Bayes’ over stepwise discriminant analysis. On the other hand, its partial agreement for within group pairing may simply point out in a quantitative way that pair assignments from the usual karyogram have only very limited meaning except for pairs 1, 2, 3, and 16 which are single pair “groups” in our terms. Because the identification of a chromosome pair is the genetically important goal, these methods and results point to a need to use better pair identification methods such as those possible with fluorescent staining procedures (9). In spite of their technical difficulties they should be performed to make detailed chromosome analysis more meaningful rather than ritualistic.

570

LOW, SELLES, AND NEURATH REFERENCES

1. ENSLEIN, K. AND NEURATH, P. W. Augmented stepwise discriminant analysis applied to two classification problems in the biomedical field. Comput. Biomed. Res. 2, 568-581 (1969). 2. MOOREHEAD, P. S., NOWELL, P. C., MELLMAN, W. J., BATTIPS, D. M., AND HUNGERFORD, D. A. Chromosome preparations of Ieucocytes cultured from human peripheral blood. Exper. Cell Res. 20,613-616 (1960). 3. NEURATH, P. W., BRAND, D. H., AND SCHREINER, E. D. Man-machine interaction for image processing. Ann. N. Y. Acad. Sci. 157,324338 (1969). 4. NEURATH, P. W., AMPOLA, M., Low, D., AND SELLES, W. Combined interactive computer measurement and automatic classification of human chromosomes. Cytogenetics 9, 424-435 (1970). 5. NEURATH, 6. 7.

8. 9.

P. W., GALLUS, G., AND SELLES, W. Interactive computer-aided chromosome analysis. J. Assoc. Adv. Med. Instr. 4, 6-14 (1970). GALLUS, G. AND NEURATH, P. W. Improved computer chromosome analysis incorporating preprocessing and boundary analysis. Phys. Med. Biol. 15,435-445 (1970). RUTOVITZ, D. Instrumentation and organization for chromosome measurement and karyotype analysis. In “Human Population Cytogenetics,” p. 292. University of Edinburgh Press, through Williams & Wilkins, Baltimore (1970). For a derivation of the generalized theorem, see H. E. Kyburg, “Probability Theory,” p. 216. Prentice Hall, Englewood Cliffs, N.J., 1969. CASPERSSON, T., ZECH, L., JOHANSSON, C., AND MODEST, E. J. Identification of human chromosomes by DNA-binding fluorescent agents. Chromosoma 30,215-227 (1970).