E Freyssingeas

M.-J Dalbe

J.-C Géminard

Flowers in flour: avalanches in cohesive granular matter

We report on the intermittent dynamics of the free surface of a cohesive granular material during a silo discharge. In absence of cohesion, one observes the formation and the growth of a conical crater whose angle is well-defined and constant in time. When the cohesion is involved, the free surface exhibits a complex dynamics and the crater, resulting from a series of individual avalanches, is no longer axi-symmetric. However, in spite of the intermittent behavior of the free surface, the flow-rate is observed to remain constant throughout the discharge.

INTRODUCTION

The flow of granular media through an orifice, because of its practical importance, has been intensely studied in the last 50 years [1][2][3][4][5][START_REF] Duran | Sands, Powders and Grains[END_REF][START_REF] Ristow | Pattern Formation in Granular Materials[END_REF][START_REF] Nedderman | [END_REF][9][10][11]. The phenomenon also presents some interesting fundamental features. For instance, what mainly differentiates the discharge of a container filled with granular matter from one filled with a viscous liquid is that the mass flow rate does not depend on the height of material above the outlet. The explanation most frequently used for this independence is based on the Janssen effect: the distribution of the weight of the material onto the silo walls, due to the friction forces, leads to a saturation of the pressure at the bottom, which results in a constant flow rate [12]. As we shall see, even if the conditions for the Janssen effect to be involved are not fulfilled, the flow rate remains constant because it is independent of the pressure in the outlet region [13]. In general, the discharge of a silo through an orifice can present three regimes: a continuous flow, an intermittent flow, or a complete blockage of the flow due to arching depending on the ratio of the typical grain diameter d to the aperture diameter D [14][15][16]. However, for a large outlet, in the continuous regime, the mass flow rate Q is satisfactorily accounted for by the empirical Beverloo law which, for instance, predicts that Q ∝ D 5/2 [1,[START_REF] Brown | Principles of Powder Mechanics[END_REF].

In addition, adding even minute amounts of liquid can change dramatically the mechanical properties of the granular matter. For instance, during the building of sand castles, one observes a transformation from a fluidlike to a sticky and deformable material with increasing water content. Indeed, at very low water content, the formation of partially developed capillary bridges leads to a fast increase of tensile strength whereas, for large enough fluid content, tensile strength is nearly constant [START_REF] Fournier | [END_REF]. Cluster formation was identfied as the main mechanism responsible of such a behavior [19]. Depending on the granular system, the cohesion can be due to electrostatic interactions and van der Waals forces, magnetic interactions [20], and even solid bonds [21,22]. The cohesion has been lately proven to give rise to puzzling surface patterning when the material is deformed [23].

Suprizingly, among the broad litterature about the silo discharge, only a few studies were dedicated to the case of cohesive materials and the dynamics of the free surface has not yet been reported. In the present article, we report a study of the surface dynamics during the continuous discharge of a weakly cohesive material and show that it exhibits a complex and interesting behavior.

EXPERIMENTAL TECHNIQUES

Experimental principle

We aim at assessing the deformation of the free surface of a granular material during a silo discharge. To do so, we consider the flow a granular material through an orifice pierced at the base of a cylindrical container. The flow-rate is measured and the free surface is observed so as to determine its geometry.

First, it is a common observation that the discharge leads to a significant deformation of the free surface only when the height of the granular column compares with the diameter of the container. In our experimental conditions, the diameter of the container is always larger than the depth of the granular bed. We note that, in this limit, the granular column is not subjected to any dynamical version of the usual static Janssen effect [3] and the pressure at the outlet thus continuously decreases during the discharge [24]. Second, it has long been identified that fluctuations of the outflow and, even, blockages are observed if the diameter of the outlet is too small, typically less than ten times the grains diameter [25]. In order to avoid system blockages, the diameter of the aperture is at least ten times the diameter of the grains in our experimental conditions. At last, we consider the case of slightly cohesive granular matter, meaning that the angle of avalanche is less than 90 degrees. The cohesion is tuned by considering various grain sizes, whereas the typical atmospheric conditions remain constant. The cohesion is accounted for by the measurement of the angle of avalanche [26,27].

Thus, in brief, the experiment consists in studying the discharge of a wet granular material through a hole pierced at the bottom of a bucket.

Experimental setup

The experimental setup (Fig. 1) consists of a cylindrical container with a hole at the bottom. The upper surface is observed thanks to two CCD cameras whereas the flow-rate is measured by means of a scale placed underneath.

FIG. 1. Experimental setup -The experimental setup consists of a cylindrical container with a hole at the bottom. The mass of the falling grains is measured by means of a scale placed underneath whereas the deformations of the free surface are observed thanks to two CCD cameras. A ring of LEDs is used to obtain a good contrast in the images from above. A LASER diode is used to produce a light sheet which makes possible to report the surface profile in a vertical, radial, plane. Two thin lamina, located at the base of the container are used to tune the size of the aperture and to open the container (see text for details).

The container is made of a Plexiglas cylinder (inner diameter 15 cm, height 15 cm) glued on a Plexiglas plate (Fig. 1) supported by four vertical legs. A large orifice (diameter 1 cm) is drilled in the plate at the center. Two additional lamina sliding in a slot (width 2 cm) are used to tune (Fig. 1, aperture) the diameter D of outlet (2, 4, 6 or 8 mm) and to suddenly initiate the flow (Fig. 1, open).

An electronic scale (KERN 440-51N) placed under the container measures the mass M of the deposited grains with an accuracy to within 0.1 g. The temporal resolution of the measurements is achieved by means of an RS-232 connection to a PC, which makes it possible to record the mass as a function of time at 2 Hz.

The free surface is imaged by means of 2 CCD cameras (JAI CB-080 GE, monitored by CVB 10.0 software from STEMMER IMAGING S.A.S.). The first camera (Fig. 1, camera 1) images the free surface from the top, along the vertical (Fig. 2a). We are interested in observing the vertical deformation of the free surface: A good contrast is achieved by the use of a ring of LEDs, placed slightly above the level of the free surface, which produces an homogeneous grazing light. The second camera (Fig. 1, camera 2) is used to obtain the profile of the free surface in one radial direction. To do so, the latter camera makes an angle of about 30 • with respect to the vertical whereas a vertical light-sheet from a LASER diode and a glass cylinder is casted onto the granular surface. In this latter configuration, the profile of the free surface can reconstructed from the deformation of the bright line visible in the image (Fig. 2b). The images from both cameras are recorded during the container discharge and subsequently analyzed with ImageJ [28].

The granular material

The granular material consists of sodosilicate glass beads (Matrasur, sodosilicate glass, density 2200 kg/m 3 ). In order to probe different situations, experiments are carried out on seven different granular batches (Table I). For all batches, the measured density of the granular ma-terial is ρ = (1510 ± 20) kg/m 3 , corresponding to a compacity of about 0.68±0.02. Using these different batches we expect to observe different behaviors. Indeed, on the one hand, for batches 4 to 7, the granular material flows easily and does not exhibit cohesion. In other words, the beads do not stick significantly to each other. On the other hand, for batches 1 and 2, the granular material is strongly cohesive and any attempt to make the material flow leads to the appearance of fractures. For these batches, the adhesion force between the grains is significant. Finally, batch 3 exhibits an intermediate state (a weak cohesion between beads). Several reasons can explain the cohesion in batches made of small grains. It may originate from capillary bridges between beads due to room humidity [26,27], or from Van der Walls interactions as well [29,30]. In our experimental conditions, both reasons are likely to be involved in the cohesion. Whatever the physical origin of the cohesion, it will be accounted for by the measurement of the angle of avalanche, θ a , which indeed depends on the adhesion force between the grains [31].

Experimental procedures

The experiment is prepared by filling, without any specific care, the container with the granular material up to a height of about 15 cm and by subsequently levelling the free surface in order to start from a flat surface and a well-defined granular height. Acquisitions of the images and of the deposited weight are started a few seconds before the aperture is opened. The experiment finishes once the grains stop flowing. Acquisitions are stopped a few seconds later.

EXPERIMENTAL RESULTS

Cohesionless material

In the case of cohesionless material, the discharge process in entirely reproducible involving a discharge at continuous flow-rate and a regular shape of the crater. We thus first report experimental results obtained for batches 4 to 7, therefore for grain diameters larger than 150 µm typically.

Discharge rate

For large grains, one observes a linear increase of the deposited mass M as a function of time t (Fig 3), and thus a constant flow-rate Q, thourough the entire discharge: We mention that, on the one hand, the stationary regime is reach immediately after the outlet is opened. On the other hand, a short transitory regime is observed just before the flow stops, at end of the reservoir drainage. Its duration, a few seconds, much less than the duration of the entire discharge, depends on the flow-rate.

In the final state an annular pile remains around the outlet. In the whole experimental range of the aperture diameter D, the experiments are reproducible and the flow rate Q only depends on the diameters d and D of the beads and outlet, respectively. It is particularly interesting to note that the flow-rate Q is constant even if the granular height above the outlet is smaller than the container diameter: even if the Janssen effect cannot insure that the pressure at the outlet is constant, the flow-rate remains constant [13]. Considering various choices of d and D (Fig. 4), we recover that the mass flow-rate Q is indeed given by the Beverloo law [1,[START_REF] Brown | Principles of Powder Mechanics[END_REF]:

Q = A ρ g 1 2 (D -k d) 5 2 (1)
where ρ = 1510 kg/m 3 denotes the density of the granular material and g the acceleration due to the gravity.

From the interpolation of the experimental data we determine the prefactor A = 0.63 ± 0.06. As expected from simple geometrical arguments, the size of the aperture is effectively reduced by a constant length of the order of the grain diameter d and we estimate k = 1.0 ± 0.1.

Q * ≡ Q/[ρ √ g(1 -d/D) 2.5
] as a function of the aperture diameter D (The symbol size is proportional to the grain diameter d). We observe that, for cohensionless materials (larger grains, batches 4 to 7), the flow-rate is satisfactorily given by Beverloo law (Eq. 1) with k = 1.0 ± 0.1 and A = 0.63 ± 0.06 (continuous line). For cohesive grains, the flow-rate varies a lot from one experimental run to another (typically up to 30% in the case of batch 2 and up to 50% for batch 1). In average, the flow-rate is smaller, given by Beverloo law with A = 0.63 ± 0.06 and larger values of k that are reported in the inset (error bars of about 10%, barely visible in the log-lin scale used to point out that k is about 1 for the cohesionless materials).

Dynamics of the free surface

For cohesionless materials, soon after the reservoir is opened, a crater forms at the free surface at the vertical of the outlet. The crater is axi-symmetric and the contour observed from above is a circle of radius R(t). For all experimental conditions, we observe that the slopes make a constant angle θ c with the horizontal. The grains avalanching continuously along the slope, we deduce that the experimental angle θ c equals, by definition, the dynamical angle of avalanche θ a . We measure to be θ a = (21.5 ± 1) deg (Fig. 5). Finally, after the discharge, when the lateral piles are at rest, we measure that θ c is slightly smaller, corresponding then, by definition, to the angle of repose θ r = (19.5 ± 1) deg.

One can easily account for the temporal evolution of the radius of the crater, R(t), by considering that the volume of the crater corresponds to the volume of the deposited mass M that increases linearly with the time t. Taking into account the flow-rate Q, the density of the material ρ and the dynamical angle of avalanche θ d , one obtains after simple algebra:

R 3 (t) = 3 Q πρ tan(θ c ) t, (2)
in good agreement with the experimental data (Fig. 6). In the next section, we will consider the behavior of the flow-rate as well as that of the free surface when the cohesion is involved in the discharge process. We will make use of the results previously reported for the cohensionless material to analyze the phenomenon in the same framework.

Cohesive material

For small grains (batches 1 and 2), typically smaller than 100 µm, when cohesion plays a role, one observes that the discharge leads to the formation of an irregular crater which results from a sequence of independent avalanches. In this case, the experiments are not, strictly speaking, reproducible. From one experiment to another, the discharge process may vary a lot. In particular, the quantity of granular material that fell through the orifice during the reservoir drainage can be very different for two experiments performed in the same experimental conditions (same grains and same outlet diameter D). In practice, we note variations of the average flow rate Q of about 30% in the case of batch 2 (d = 45 -90 µm) and of about 50% for the smallest grains (batch 1, d = 0 -45 µm), from one experimental run to another. This may come from the difficulty of preparing identical initial states. Indeed, the way the beads are poured in the tank significantly affects the initial compaction, and thus the discharge process, of cohesive granular materials. Nevertheless, the typical behavior of the material during the discharge exhibits several common features that we report in the present section.

Discharge rate

Although the dynamics is intermittent at the free surface, the mass of granular material flowing out through the hole still exhibits a linear increases, with little fluctuations, as a function of time thorough the entire discharge (Fig. 7). Therefore, as for cohesionless material, Q is constant. In this case again, the flow rate, which depends on the diameters d and D of the grains and of the outlet respectively, is independent of the amount of material above. It is particularly interesting to point out that, surprisingly, the flow rate is constant in spite of the intermittent dynamics of the avalanche process above. The explanation of such behavior is found by direct observation of the crater. Indeed, one easily observes that the drainage involves two different processes. First, the grains that detach from the lateral slopes fall down in the crater. This fluidized material flows continuously through the outlet. Second, the decrease of the height of fluidized material in the crater leads to a destabilization of the lateral walls of the crater. The resulting avalanches feed, in turn, the fluidized region above the outlet. The flow rate which does not depend on the height of grains above the outlet remains constant.

For cohesive grains, the experimental flow-rate Q is observed to be significantly smaller than that expected with the values A = 0.63 and k = 1 (cohesionless samples) in Beverloo law (Fig. 4). However, the experimental data from the cohesive samples are correctly accounted for by Eq. 1 with A = 0.63 ± 0.06 (∀d), but with k increasing with the cohesive nature of the material (small grain size, Fig. 4,inset). This might be because the cohesion is only partially destroyed when the grains flow. Due to the remaining cohesion, grain aggregates, rather than individual grains, flow out through the aperture. In this case, the effective diameter d eff of the falling grains is larger than the individual grain diameter d, leading to a smaller effective outlet diameter [ie to a larger value of k in Berverloo law (1)].

It is particularly interesting to comment here the measurements of the flow rate, even if they were not the scope of the present article. Indeed, the fact that the Q is constant throughout a silo discharge is not an obvious result. In the litterature, the constancy of Q has been often implicitely assumed to be a result of the so-called Janssen effect, a saturation of the pressure at the base of a static silo. Indeed, in this case, the pressure in the outlet region would be constant, independant of the height of material above. Even if the authors generally did not explicitely state that there was a relation between the constancy of the flow-rate and the Janssen effect, they made erroneously think that the flow-rate was constant because the pressure in the outlet region was constant. Note first that the Janssen effect is a static effect and that there is no reason for the pressure to be constant when the material is flowing. We now know that, indeed, the pressure is not constant during the discharge, as proven by direct measurements of the forces exerted by the material on the container [24,32]. Moreover, in Ref. [32], the authors show that the flow rate is not altered by vibrations of the container that avoid static arching. Independently, it has been proven in a horizontal system (a conveyor belt) that the flow-rate does not depend on the pressure and obeys Beverloo law, the typical velocity at the outlet to consider in this case being the belt velocity instead of the free fall velocity √ gD [13]. Our experiment interestingly provides an additional experimental proof of the same assertion, the flow-rate being constant even if the height of the granular column is never larger than the radius of the container. The pressure surely varies dur-ing the discharge but Q is constant and accounted for by Beverloo law. Furthermore, we observe that the cohesion mainly alters, k, thus the typical size of the particles to be considered; the size of grain aggregates rather than the diameter of individual grains.

The physical process leading to the reduction of the flow rate being out of the scope of the present article, we do not discuss further this point and focus now on the dynamics of the crater.

Dynamics of the free surface

For cohesive material, the crater geometry and dynamics are noticeably different than those reported for the cohesionless case.

First, we observe a significant delay between the beginning of the discharge and the appearance of a crater at the free surface. Because of the cohesion, the material can sustain a cavity at the vertical of the outlet and the deformation of the free surface result from the granular flow at the bottom only when a significant amount of grains has flown out. Once the crater appears at the free surface, the angle that the lateral slopes make with the horizontal is large (Fig. 8). Note that, as the grains at the lateral slopes are not flowing, the crater angle θ c does not correspond to the dynamical angle of avalanche anymore. We also point out that the crater angle does not depend on the diameter D of the outlet. However, contrary to what is observed for non-cohesive grains, the angle depends on the grain size d, i.e. on the cohesion (Table I): From the vertical profiles of the crater (Fig. 8), we obtain, for instance, θ c = (88 ± 2) deg for batch 1 (0 -45 µm) and θ c = (56 ± 1) deg for batch 2 (45 -90 µm). In agreement with the expectation of a larger effect of the cohesion force on the smaller grains, the angle of avalanche significantly increases when the grain size is decreased [31]. From now on, we shall consider the crater angle θ c as the quantity accounting for the cohesion. As mentioned above, the crater is no more circular and its growth is due to successive, partial, collapses of the lateral walls. This behavior is clearly revealed by the dynamics of the contours obtained from the images from top (Fig. 9, inset). The superposition of the contours at several equally-spaced times reveals that the intermittent nature of the discharge. In such an image, in which the superposition of the contours draws a "flower" at the surface, one directly observes that the size of the avalanches and the delay (lines are darker away from the center) between successive events increase with time. In order to report quantitative results, we define the average radius R of the crater, R ≡ A/π, where A is the surface area of the crater in the horizontal plane. Reporting R as a function of time (Fig. 9), we observe that R increases by jumps, each corresponding to a partial collapse of the lateral wall. Again, the results point out that both the typical size of the avalanches and the delay between them increase thoughout the discharge process. In addition, one observes, in average, the scaling R ∝ t 1/3 , which is compatible with a constant flow-rate and constant angle θ c of the lateral walls.

Geometrical characterization of the "flower"

To characterize further the discharge process, we now characterize the geometry of the crater seen from above.

In particular, we observe that, during the discharge process, the typical size (arc length) of the avalanches increases with the distance to the center so that one can consider the angular size of the avalanche, α (the angle that the straight lines joining the centre to the two ends of the avalanche make between them, Fig. 10). Such an analysis is possible only when the material is sufficiently cohesive so that one can distinguish individual avalanches, thus for the smallest grains. We thus report results obtained for batches 1, 2 and 3, only. We thus consider the distribution of α measured throughout several discharges (typically, 200 to 300 avalanches). We observe that the distributions of α do not depend significantly on the aperture diameter D and report results obtained for D = 6 cm (Fig. 11). The distributions are broad and asymmetric but we note that the typical angular size is significantly larger for the smaller grains. For instance, the experimental results lead to α = (73±23) deg for d = 0-45 µm and α = (54±17) deg for d = 45 -90 µm (Fig. 11). For the cohesionless samples such an estimate cannot be obtained whereas for the intermediated case d = 75 -150 µm the individual avalanches can barely be isolated and we can only get an estimate α = 17 deg with an error of about 100 % (Fig. 12). We observe that α → 0 for vanishing small cohesion (θ c = 22 deg), which can be understood by considering that the continuous collapse of the walls in this case consists in an infinite number of infinitely small avalanches.

Additional characterization of the discharge dynamics

In order to characterize further the intermittent nature of the discharge, we define an estimate, M (t), of the mass of material that fell down the slopes at time t. To do so, we use the apparent radius R(t) and the angle θ c of the crater and write, considering the volume of the corresponding cone, M = ρ πR 3 3 tan(θ c ).

(3)

Then, from the estimated mass M (t) and the actual mass m(t) that is directly measured with the balance, we define:

Φ = 1 T m 2 ∞ T 0 [m(t) -M (t)] 2 dt ( 4
)
where T stands for the total duration of the discharge and m ∞ for the total mass of material measured after the discharge. The parameter Φ accounts for the intermittent nature of the discharge. Indeed, M (t) = m(t) and thus Φ = 0 for a cohesionless material. For cohesive materials, due to the avalanches at the free surface, the temporal fluctuations in R (and thus in M ) are enhanced whereas the flow rate is constant and m increases linearly with time. Thus, Φ is expected to increase with the intermittent nature of the discharge process and, thus, to account for the cohesion. In the figure 11, we report the experimental values of Φ as a function of the crater angle θ c . We note that, even if it, indeed, increases significantly for the cohesive sample, Φ gently vanishes for vanishing cohesion: the occurrence of small avalanches observed in the case of batch 3 does not produce any important oscillation of the estimated mass M (t) and the parameter Φ is not very appropriate to reveal the cohesive nature of the material. The measure of the angular size α is much more efficient for this purpose as it significantly depends on the crater angle θ c even for almost cohesionless materials.

DISCUSSION AND CONCLUSION

The experimental results can be summarized as follows: for large enough aperture, the discharge is continuous in the sense that the flow rate is constant throughout the discharge. The results holds true in absence or presence of cohesion. However, when the material is (even weakly) cohesive, the dynamics of the free surface is intermittent. Instead of an axi-symmetric conical crater, one observes an irregular crater which results from a series of individual avalanches. The typical size of the avalanches increases with the crater size, which permits to define, for a given material, an angular size which increases with the cohesion. These results suggest that there is a geometrybased relation between the angular size of the avalanche and the angle that the crater walls make with the horizontal. However, we have not been able to establish this relation theoretically.

In the effort to understand, from the fundamental point of view, the mechanisms governing the discharge of cohesive materials, the present work shall be extended to measurements of the forces exerted on the base plane and on the wall while tuning the cohesion.

FIG. 2 .

 2 FIG. 2. Images of the free surface -(a) Top view -A good contrast is obtained by the use of an annular light source providing an homogeneous grazing light. The bright line is from the vertical light sheet casted onto the surface. (b) Side view -The deformation of the bright line, observed from the side thanks to the second camera makes it possible to reconstruct the surface profile along one radial direction (d = 45 -90 µm, D = 6 mm) .

  45-90 75-150 150-200 180-300 250-425 425-600

FIG. 3 .

 3 FIG.3. Deposited mass M vs. time t -For cohesionless materials, the deposited mass M increases linearly with time t thourough the entire discharge. The permanent regime is reached immediately after the outlet is opened and the flowrate Q can be accurately determined from the slope. The final transitory regime, during typically 10 seconds in this example, corresponds to the fall of the last grains along the granular slopes that remain around the outlet in the final state (d = 250 -425 µm, D = 6 mm).

FIG. 4 .

 4 FIG. 4. Reduced flow-rateQ * vs. D -We report the reduced flow-rate Q * ≡ Q/[ρ √ g(1 -d/D)2.5 ] as a function

FIG. 5 .

 5 FIG.5. Vertical profile at various times t -We report the profile of the crater in a vertical plane of symmetry every 5 seconds. We observe that the angle θc that the slopes make with the horizontal is constant throughout the entire discharge (Thick lines: grey, profile of the initial surface; black, t = 25 s. Thin lines: from light grey to black, t = 5 s to 20 s . d = 425 -600 µm, D = 6 mm).

FIG. 6 .

 6 FIG. 6. Radius R vs. time t -In agreement with Q constant and simple geometrical arguments, R scales like t 1/3 [Black dots, experimental data; grey line, Eq. (2)]. Inset: contours of the crater in the horizontal plane at different, equally-spaced (2.5 sec), times t (The scale bar is worth 2 cm. d = 150 -250 µm, D = 6 mm).

FIG. 7 .

 7 FIG. 7. Deposited mass M vs. time t -For cohesive materials, the deposited mass M increases linearly with time t thorough the entire discharge in spite of the intermittent avalanche process occurring in the container (d = 45 -90 µm, D = 6 mm).

FIG. 8 .

 8 FIG.8. Vertical profile at various times t -We report the profile of the crater in a vertical plane every 5 seconds. We observe that the slopes on the side are much steeper than in the cohesionless case (Thick lines: grey, profile of the initial surface; black, t = 20 s. Thin lines: from light grey to black, t = 5 s to 15 s . d = 45 -90 µm, D = 6 mm).

FIG. 9 .

 9 FIG. 9. Effective radius R of the crater vs time t -The average radius of the crater increases by independant jumps. However, in average, R scales like t 1/3 . The vertical dashes mark the avalanche events observed directly in the images [Black dots, experimental data; grey line, Eq. 2 with θc = 56 deg, ρ = 1510 kg/m 3 and the measured dimensionless flowrate Q * = 1.75]. Inset: superposition of the crater contours in the horizontal plane at different, equally-spaced, times t. The darker the line is, the longer the contour remained at the given position (The scale bar is worth 2 cm. d = 45 -90 µm, D = 6 mm).

FIG. 10 .

 10 FIG. 10. Geometrical characterization of the "flower" -We determine the angular size α of each of the individual avalanches (d = 45 -90 µm, D = 6 mm).

FIG. 11 .

 11 FIG. 11. Angular size α and parameter Φ vs crater angle θc -The typical angular size of the avalanches, α and the parameter Φ (Eq. 4) increase significantly with the typical angle θc that the lateral walls make with the horizontal, thus with the cohesion (The cross indicates θc = 22 deg for cohesionless batches and the continuous lines are only a guide to the eyes). Inset: Distribution p(α) for d = 45 -90 µm (D = 6 mm).

TABLE I .

 I Granular

	•	•	•/•	•	•	•	•
	88	56	34	22	22	22	22

batches -Seven different batches, consisting of glass spheres, were used in the experiments. Small grains clearly exhibit cohesion whereas large grains do not exhibit any apparent cohesive behavior. Batch 3 is an intermediate case (Top: batch number. Center : diameter range in µm. Center : • cohesion; • no apparent cohesion. Bottom : lateral slopes angle θc.)


Keywords: Composite Likelihood, Association Mapping, Breast Cancer, Imputed Genotypes, FGFR2 gene

We describe composite likelihood-based analysis of a genome-wide breast cancer casecontrol sample from the Cancer Genetic Markers of Susceptibility project. We determine 14,380 genome regions of fixed size on a linkage disequilibrium map which delimit comparable levels of linkage disequilibrium. Although the numbers of SNPs are highly variable each region contains an average of ~35 SNPs and an average of ~69 after imputation of missing genotypes. Composite likelihood association mapping yields a single P-value for each region, established by a permutation test, along with a maximum likelihood disease location, standard error and information weight. For single SNP analysis the nominal P-value for the most significant SNP (msSNP) requires substantial correction given the number of SNPs in the region. Therefore imputing genotypes may not always be advantageous for the msSNP test, in contrast to composite likelihood. For the region containing FGFR2 (a known breast cancer gene) the largest chi-square is obtained under composite likelihood with imputed genotypes ( 2 2 χ increases from 20.6 to 22.7), and compares to a single-SNP based 2 2 χ of 19.9 after correction. Imputation of additional genotypes in this region reduces the size of the 95% confidence interval for location of the disease gene by ~40%. Amongst the highest ranked regions, SNPs in the NTSR1 gene would be worthy of examination in additional samples. Meta-analysis, which combines weighted evidence from composite likelihood in different samples, and refines putative disease locations, is facilitated through defining fixed regions on an underlying linkage disequilibrium map.

INTRODUCTION

Genome-wide association mapping studies based on large case-control samples [START_REF] Easton | Genome-wide association study identifies novel breast cancer susceptibility loci[END_REF][START_REF] Hunter | A genome-wide association study identifies alleles in FGFR2 associated with risk of sporadic postmenopausal breast cancer[END_REF] have identified common genetic variants associated with increased risk of breast cancer. Most analyses of genome-wide case-control data sets employ tests based on individual single nucleotide polymorphisms (SNPs) [START_REF]WTCCC): Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls[END_REF] . Meta-analysis (combining evidence across samples) is facilitated by imputation of 'missing' SNP genotypes, using the HapMap samples (http://www.hapmap.org/ ) as a reference population [START_REF] Zeggini | Meta-analysis of genome-wide association data and largescale replication identifies additional susceptibility loci for type 2 diabetes[END_REF] . An alternative approach to single SNP tests [START_REF] Gibson | A multimetric approach to analysis of genome-wide association by single markers and composite likelihood[END_REF][START_REF] Morton | Genome scanning by composite likelihood[END_REF] undertakes composite likelihood analysis of multiple SNPs in a region and determines a location for a putative disease influencing variant on an underlying linkage disequilibrium unit (LDU) map [START_REF] Maniatis | The first linkage disequilibrium (LD) maps: Delineation of hot and cold blocks by diplotype analysis[END_REF] . When plotted against physical (kb) locations the LDU map describes the underlying pattern of linkage disequilibrium (LD) as a series of plateaus (strong LD) and steps (where LD is breaking down, such as at the location of recombination hot spots). The LDU map provides a framework for characterising small chromosome regions which may differ substantially in physical size but share comparable levels of linkage disequilibrium. Modelling the pattern of association with disease at multiple markers in a region generates a single P-value for disease association, a disease location, standard error and corresponding information weight. As there is just one statistical test in a region there is a reduced Bonferroni correction relative to single SNPbased tests which require consideration of the number of tests made at nearby SNPs.

Gibson et al [START_REF] Gibson | A multimetric approach to analysis of genome-wide association by single markers and composite likelihood[END_REF] evaluated the composite likelihood approach using relatively low density genotype data (~200K SNPs) in a relatively small sample (403 cases and 395 controls) with an undisclosed disease phenotype. Larger and more comprehensively genotyped samples are now available. The genome-wide breast cancer association analysis by Hunter et al [START_REF] Hunter | A genome-wide association study identifies alleles in FGFR2 associated with risk of sporadic postmenopausal breast cancer[END_REF] utilised samples from 1,145 postmenopausal women of European ancestry with invasive breast cancer contrasted with 1,142 controls analysed with 528,173 SNPs. These data are made available through the Cancer Genetic Markers of Susceptibility (CGEMS) project data portal (http://cgems.cancer.gov/). The data are valuable for comparing composite likelihood and single marker analyses and for developing strategies for meta-analysis. These data present significant evidence for a now well established breast cancer gene, FGFR2, which has been verified in several studies 1 . We describe the application of a composite likelihood modelling approach to this higher density SNP sample, evaluate relative power for composite likelihood and single-SNP based tests and test the impact of increasing marker coverage through genotype imputation. The chromosome-region based approach used in composite likelihood, with regions defined on the underlying LDU map, is highly suited to meta-analysis which is essential to increase the sample size for the identification of novel causal variants.

MATERIALS & METHODS

Data Preparation & Quality Control

Following successful application for permissions, data comprising 1,145 cases and 1,142 controls and genotypes for 555,148 SNPs were downloaded from the CGEMS data portal.

Data files were converted into PLINK format [START_REF] Purcell | PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses[END_REF] and quality control (QC) procedures undertaken. Samples rejected through the QC employed by Hunter et al 2 had already been excluded in the downloaded dataset from an original set of 1,183 cases and 1,185 controls.

The QC we applied resulted in the removal of 93 SNPs with inconsistent or ambiguous kilobase locations, 8,648 SNPs with a high proportion (>10%) missing genotypes, 53,615 SNPs with minor allele frequencies lower than 0.05 and a further 4,308 SNPs with large deviations from Hardy-Weinberg (HW, χ 2 ≥ 10) in the controls (Supplementary Table 1). In addition, one individual with >10% missing genotypes was excluded at the QC stage. To minimise biases created by population stratification we identified individuals with possible non-Caucasian ancestry through multidimensional scaling cluster analysis 8 (Supplementary Figures 1 and2) using 73,560 "LD-independent" SNPs from CGEMs and HapMap. A total of 12,907 of these SNPs showed strand mismatches and were flipped accordingly. No A/T or G/C SNPs were genotyped in the CGEMS data because of the chemistry of the genotyping beadchip (Infinium II). This cluster analysis identified four individuals which were judged to be outside the CEU cluster, suggesting admixture, and were excluded from further analysis at this point. Following QC we analysed a total of 498,786 SNPs in 1,143 cases and 1,139 controls.

Genotype Imputation

After flipping strands for 94,489 SNPs, to ensure strand concordance of the two SNP data sets, a combined CGEMS and HapMap (CEU, phase 3) dataset was produced for genomewide genotype imputation using the PLINK software. 3).

Composite Likelihood Tests

The program CHROMSCAN 9 develops the model described by Maniatis et al 10 utilising data from SNPs in a chromosome region to compute a maximum likelihood location, S, for a causal variant, a standard error, a 95% confidence interval and a permutation-based Pvalue. The underlying LD structure is incorporated into the model through LDU maps, which represent the association mapping analogue of the linkage map [START_REF] Collins | Mapping genes for common diseases: the case for genetic (LD) maps[END_REF] . Disease mapping on the underlying LDU scale has been shown to increase fine mapping resolution and power relative to the physical (kilobase) map [START_REF] Collins | CHROMSCAN: genome-wide association using a linkage disequilibrium map[END_REF] . We constructed LDU maps from the CEU sample from HapMap Phase II based on physical locations from build 36 of the human genome sequence (University of California, Santa Clara, March 2006). Genome-wide LDU maps [START_REF] Lau | Exploiting large scale computing to construct high resolution linkage disequilibrium maps of the human genome[END_REF] are available on request from the authors.

The association test reduces the 3 2 × table of SNP genotype counts by disease affection status at the i th SNP to the corresponding 2 2 × table of allele counts by affection status, with cell totals a, b, c, d, giving n haplotypes from 2 n diplotypes. Association of disease phenotype with SNPs in a region is modeled using a composite likelihood approach.

Observed association with disease at the i th SNP is:

( ) ( ) ˆ( ) ( ) i a d b c z a b b d × -× = + × + , with information ( ) ( ) ( ) ( ) i z a b b d K n a c c d + × + = × + × + .
Expected association, i z , is modeled using the Malecot equation 10 :

( )

1 i S S i z L Me L ε - - = - +
, where i S is the location of the i th SNP in LDU and the S parameter represents the LDU map location showing maximal association with disease. The ε parameter describes the decline of association with map distance and has a value ~1 for LDU maps [START_REF] Collins | CHROMSCAN: genome-wide association using a linkage disequilibrium map[END_REF] , M is the intercept, and L is the asymptote, representing association not due to linkage which is estimated ( L ) or predicted ( p L ). The predicted asymptote is taken as the mean absolute value of a standard normal deviate, weighted by information z K [START_REF] Collins | CHROMSCAN: genome-wide association using a linkage disequilibrium map[END_REF] . Composite likelihood is defined as:

2 lk e -Λ =
, where ( ) The association test statistic for each region is the difference

2 î z i i i K z z Λ = - ∑ . CHROMSCAN
A D X = Λ -Λ . The difference
is computed for the real data ( 1 H ) and a large number of replicates ( 0 H ), as j X for the j th replicate, in which the disease phenotype is randomised (shuffled). The distribution of P-values under 0 H is obtained from fractional ranks in a large sample of replicates. From each of the replicate P-values the corresponding 2 3 χ for the contrast between models 'A' and 'D' is obtained from the GNU Scientific Library (GSL) function gsl_cdf_chisq_Pinv (http://www.gnu.org/software/gsl/ ), and hence the variance for the j th replicate is:

2 3 X V χ = .
Variances for replicates, j V , , are used to predict, by regression, variance V ( 1 H ) and hence 2 3

χ ( 1 H ). The computation of V ( 1 H ) requires a sorted sub-set of replicates which are centered on the value X ( 1 H ), and the model: lnV j = A + B lnX j , with X centered between the 20 closest replicates with X j ≤ X and the corresponding 20 with X j ≥ X; if X is an outlier, the 20 closest values are taken. From this model V ( 1 H ) is estimated as exp(A + B ln X), and

2 3 χ ( 1 H ) = X/V.
Simultaneous estimates of M , S and L give an information matrix which is inverted to provide the nominal variance (U) for location S. Using V ( 1 H ), the information weight, W , about disease gene location, S , is computed as:

1 3 U W V =
and the standard error ( SE ) of

S is: 1 SE W = .
We revised CHROMSCAN to increment the number of replicates adaptively to ensure that the P-value ( 1 H ) predicted from the replicates is accurately determined, with a minimum of 50 replicates and maximum of 20000 per region (or more for refining evidence in a significant region of interest). Gibson et al [START_REF] Gibson | A multimetric approach to analysis of genome-wide association by single markers and composite likelihood[END_REF] , in their analysis of a relatively low density SNP data, used non-overlapping regions spanning at least 10 LDUs and containing a minimum of 30 SNPs. More recent high-density panels enable analysis in smaller regions and higher resolution with reduced possibility of confounding adjacent independent signals.

We used regions of fixed LDU size which facilitates combination of evidence in metaanalysis. Regions of four LDUs contain an average of over 30 SNPs in a ~550,000 SNP scan, assuming ~60,000 LDUs in the CEU genome [START_REF] Lau | Exploiting large scale computing to construct high resolution linkage disequilibrium maps of the human genome[END_REF] . However, there is wide variation in the number of SNPs per region, although coverage is increased with genotype imputation (Table 1).

Single SNP Tests

For single SNP tests we identify the most significant (ms) SNP in a region, from the nominal [START_REF] Hunter | A genome-wide association study identifies alleles in FGFR2 associated with risk of sporadic postmenopausal breast cancer[END_REF] 1 χ (from the 2 2 × table between affection and the two SNP alleles). Selecting the msSNP from a large number of SNPs in a region biases the nominal P-value ( n P ), computed on the null hypothesis. To correct for the number of SNPs we firstly grouped four-LDU regions into ranges which show relatively limited diversity in the number of SNPs they contain (Table 1). The ranges (SNP range, Table 1) were defined to include approximately similar numbers of four-LDU regions, with the exception of regions containing more than 250 SNPs. This enabled the relationship between T (the observed mean number of SNPs in the range) and R (the effective mean number of SNPs in the range) to be characterized in the tail of the distribution. We determined the distribution of numbers of SNPs in each of 28,750 four-LDU regions (original and imputation inclusive data sets combined) and computed the weighted mean number of SNPs, T , in a range (for each range

i i i f m T f = ∑ ∑
, where f is the number of four-LDU regions containing m SNPs, with summation over i = 1, N regions, Table 1). Under the null hypothesis P-values for random SNPs have expectation 2 2 2 ln P χ = , with an expected variance of four and a mean of two [START_REF] Gibson | A multimetric approach to analysis of genome-wide association by single markers and composite likelihood[END_REF] . For each range we computed, R, the effective number of independent SNPs (Table 1) by regula falsi. Bonferroni correction assumes a corrected P-value χ from the msSNPs (μ = 2.5 and 2.2 for original and imputation datasets respectively) and the 2 2 χ , from permutation-based P- values in composite likelihood analyses (μ = 1.9 for both original and imputation inclusive datasets), were multiplied by 2/μ to correct the deviation from the expected mean of 2.

RESULTS

CHROMSCAN analysis yields 14,370 four-LDU regions containing at least one SNP from the original genotype data and 14,380 from the data containing imputed genotypes. The distribution of SNPs in each region is very variable (Table 1, Figure 1). Many regions contain 20 or fewer SNPs and, since the LDU map describes regions with comparable levels of linkage disequilibrium, this suggests that a substantial proportion of the genome may be poorly screened by this set of genotypes. Coverage is increased by imputation of missing genotypes, with the mean number of SNPs per region increasing from 34.7 to 68.6

(Figure 1). However, following imputation ~15% of the regions still have 20 or fewer SNPs and may be poorly represented by both single SNP and composite likelihood tests. SNP panels with more uniform coverage of markers on the LDU, rather than kilobase scale, would reduce the possibility of overlooking regions associated with disease. In higher SNP density panels the magnitude of the Bonferroni correction required for single SNP analysis will be greater. In contrast, more comprehensive genotyping may increase power for composite likelihood tests because one permutation-based P-value is obtained for every region.

The distribution of nominal single SNPs ( 2 1 χ ) in the FGFR2 gene region (Figure 2) show a cluster of SNPs localized in a region with extensive linkage disequilibrium, represented as a plateau on the underlying LDU map. Composite likelihood mapping in this region (Table 2) indicates that, after imputation adds ~64% more SNPs, there is an increase in 2 2 χ from 20.6 to 22.7. The 95% confidence interval for the location of the causal variant decreases by 40% from 1.5 to 0.9 LDUs using the more densely genotyped imputation data set. This reduction in the confidence interval, which spans intron 2 of FGFR2, is reflected in the composite likelihood surface (Figure 3) which shows the difference in

A D X = Λ -Λ
between the A (null) and D (causal variant location) models for the original and imputation inclusive data sets.

DISCUSSION

Comparison of composite likelihood and single SNP tests suggest higher power of the former for the FGFR2 association, which is well established as breast cancer risk gene.

Power is further increased with imputation of missing genotypes (Table 2). None of the other genes identified in Table 3 contain well established breast cancer risk variants although it is notable that the NTSR1 (neurotensin receptor 1) gene ranks highly in both composite likelihood and single SNP tests. NTSR1 is a candidate risk factor involved in ductal breast cancer progression [START_REF] Dupouy | The Neurotensin Receptor-1 Pathway Contributes to Human Ductal Breast Cancer Progression[END_REF] . The authors note that in breast cancer cells functionally expressed NT1 receptor coordinates transforming functions including cellular migration and invasion. High expression of NTSR1 is associated with the tumour grade, size and number of metastatic lymph nodes. Given that the well established breast cancer genes only account for a small proportion of the familial genetic risk, regions that fail to achieve genome-wide significance, but rank highly, are worthy of examination in larger samples. A worthwhile focus of future analyses includes screening highly ranked variants in breast cancer phenotypic sub-types, including those which describe tumor characteristics [START_REF] Tapper | The influence of genetic variation in 30 selected genes on the clinical characteristics of early onset breast cancer[END_REF] .

Hunter et al 2 describe the original analysis of these data and the identification of SNPs in the FGFR2 gene as highly associated with sporadic postmenopausal breast cancer. These findings were confirmed by the authors in a second sample. Although strong evidence for the involvement of FGFR2 is a feature of our analysis, comparison with the results presented by Hunter et al is difficult. Differences in the quality control procedures employed (Supplementary Table 1), their use of additional phenotypic data (details of age and hormone replacement therapy use) and differences in analytical methods employed, including their use of logistic regression models, underlie the difficulty of comparison.

In the FGFR2 gene region the apparent higher power for composite likelihood tests must be achieved partly by modeling association at multiple SNPs simultaneously. Alternative approaches that combine data from multiple SNPs include haplotype-based tests [START_REF] Huang | Detecting haplotype effects in genomewide association studies[END_REF] . Such approaches have the advantage of modeling correlations between markers, potentially increasing power, along with the characterization of genetic effects on different haplotypic backgrounds. The disadvantages include the difficulty in deciding how to define haplotype 'windows', the heavy computational requirements, lack of a clearly defined disease interval that is refined with accession of data and the difficulty of combining evidence across samples. Imputation of genotypes, which can usefully increase coverage and potentially provide further increases in power, must also increase the computational and multipletesting burden for haplotype tests, in line with single-SNP based analyses.

Some authors have found that imputing genotypes is rather accurate [START_REF] Hao | Accuracy of genome-wide imputation of untyped markers and impacts on statistical power for association studies[END_REF] but note that power increases only slightly as imputation 'results in modest gain in genetic coverage, but worsens the multiple testing penalties'. This penalty is likely to further erode power when using more comprehensive SNP panels and with imputation at higher densities, as might be achieved (for example) using data from the 1000 genomes project (http://genome.wellcome.ac.uk/doc_WTX047611.html ). Other authors note that the typical imputation error rates of 2 -6% 17 may substantially decrease power and so the utility of this technique may be questioned for single SNP-based analyses.

As individual genetic effect sizes are generally low for common variants involved in complex traits, meta-analysis combining evidence across studies, is an important strategy to increase power and identify novel targets for further follow-up [START_REF] Zeggini | Meta-analysis of genome-wide association data and largescale replication identifies additional susceptibility loci for type 2 diabetes[END_REF] . A composite likelihoodbased approach, in which association evidence from different genome-wide association samples is combined across corresponding regions, provides a test in which individual samples are weighted according to their information ( W , Table 2). This approach also gives an estimate of disease gene location which becomes more precise as further evidence is combined [START_REF] Tapper | Mapping a gene for rheumatoid arthritis on chromosome 18q21[END_REF] . The methods presented here provide a strategy for the analysis of component samples in such a meta-analysis taking advantage of genotype imputation to increase coverage without increasing the multiple-testing penalty.

Polymorphisms in intron 2 of the FGFR2 gene have been implicated as increasing risk of breast cancer in European and Asian populations. Easton et al 1 reported two SNPs, rs2981582 and rs7895676 (at the upstream and downstream boundaries respectively of intron 2), as the most strongly associated and suggested that the latter was most likely to be a causal variant as it showed the strongest association with breast cancer risk. Recently, Boyarskikh et al [START_REF] Boyarskikh | Association of FGFR2 gene polymorphisms with the risk of breast cancer in population of West Siberia[END_REF] , studying a West Siberian population, noted that rs2981582 explained association with disease much more strongly than rs7895676. The authors hypothesized that the actual causal variant lies somewhere within the LD block which includes these two SNPs. Although rs7895676 (location 123323.987 kb) is not represented in the imputationinclusive data set, and rs2981582 does not have the highest single marker chi-square in the sample (Table 3), these markers flank the cluster of associated SNPs in the intron 2 LD block (Figure 2). Given that intron 2 lies within a strong LD block fine mapping to confirm the location of the causal variant will be facilitated by meta-analysis in which the appropriately weighted accessions of data should enable further reduction of the target confidence interval. W: Information weight for location S (LDU).

95% CI: 95% confidence interval for disease gene, size in brackets.

-2lnP: -2 log P-value for association with disease ~2 2 χ .

Titles & Legends to Figures

  evaluates the composite likelihood for two subhypotheses to test the evidence for a disease associated variant in a region. Within a given region the null hypothesis ('model A') assumes only 'background' association and no relationship between the affection status and SNPs with: , 0 p L L M = = . As the null model does not test association with disease there is no location estimate, S . The null model, which estimates no parameters, is contrasted with 'model D', which estimates three parameters: a disease location ( S ), an intercept ( M ), and background association ( L ).

Figure 1 :

 1 Figure 1: The distribution of SNPs in each four-LDU window in the original and imputed

Figure 2 :

 2 Figure 2: Nominal single SNPs chi-squares for association and the LDU map of the FGFR2

Figure 3 :

 3 Figure 3: The difference X in (composite) likelihood, where

  values from single SNPs we determined the relationship between R and T by regression such that a value R could be assigned to each four-LDU region, given T.

				c P	1 (1 = --	n P	) R	. To
	correct P-Regression through the origin gives	s R	=	(0.306239 ) (0.000248 T × -	×	T	2	)	(model R-

squared = 0.96) which enables Bonferroni corrected values c P to be computed. The Bonferroni correction greatly reduces the significance of the nominal P-values. Composite likelihood tests do not require this correction as P-values are based on a permutation test. Following correction of msSNP P-values for the variable number of SNPs in individual four-LDU regions the means (μ) of the corrected 2 2c

Table 3

 3

	is ordered by the 10 most significant regions identified using composite likelihood
	in the imputed data set. The FGFR2 region is highest ranked for both composite likelihood
	and single SNP tests. Power, as indicated by 2 ln P -	(= 2 2 χ ) appears relatively lower in
	these data for single SNP tests compared to the composite likelihood based analysis. There

is quite strong correspondence between ranks in the original and imputed data sets for the five highest ranked regions but less agreement for regions ranked 6-10. There is reduced correspondence between single-SNP and composite likelihood tests, although the NTSR1 gene region has relatively high ranks for both tests.

Table 1 -The distribution of SNPs within four LDU regions

 1

	SNP range	N	T	R
	1 -5	1120	2.22	1.67
	6 -10	1285	7.13	4.58
	11 -14	1243	11.57	6.98
	15 -18	1449	15.52	9.28
	19 -22	1595	19.52	10.46
	23 -25	1223	23.01	12.32
	26 -28	1284	25.97	14.13
	29 -32	1735	29.50	15.83
	33 -36	1647	33.50	17.24
	37 -39	1171	36.96	16.88
	40 -43	1433	40.49	18.74
	44 -48	1723	44.96	20.10
	49 -53	1496	49.98	23.68
	54 -59	1442	55.45	20.91
	60 -66	1436	61.86	24.56
	67 -75	1526	69.81	25.51
	76 -87	1536	80.19	28.04
	88 -104	1391	94.40	27.52
	105 -136	1527	117.71	30.70
	137 -250	1348	171.24	37.33
	251 -550	135	300.16	55.64
	551 -800	5	656.80	97.20
	Total	28750	-	-

N: Total number of four LDU regions containing the number of SNPs in range, T: weighted mean number of SNPs in the range, R: mean effective number of SNPs in the range (Bonferroni correction). Table

contains

regions from original and imputation data sets combined.

Table 2 -Composite Likelihood Analysis in the FGFR2 Region

 2

		Number of SNPs S LDU	W	95% CI in LDU	95% CI in Kb	-2lnP
	Original dataset	33	2379.66	7.65	2378.9-2380.4	123311.2-123357.4	20.63
			(0.36)		(1.5)	(46.2)
	Imputation dataset 54	2379.67	17.55	2379.2-2380.1	123315.6-123347.9	22.68
			(0.24)		(0.9)	(32.3)

Key: S LDU: Maximum likelihood location for disease gene (LDU), standard errors in brackets.

Table 3 -Composite Likelihood: The ten highest ranked regions (imputation-inclusive dataset)

 3

	Chromsomes

"Imputed" here refers to the aggregated dataset consisting of both original and imputed SNPs.

Gene including the msSNP identified.

Only 555,148 SNPs out of the 555,352 originally genotyped SNPs were available in the CGEMS Portal datasets.

Only SNPs with HW χ 2 ≥ 10 in the Nurses' Health Study [NHS] control group were excluded from further analyses.

ACKNOWLEDGEMENTS

The work was funded by the Breast Cancer Campaign.

Conflict of interest.

The authors declare no conflicts of interest. Supplementary Figure 2: Enlarged version of the multidimensional scaling plot (Supplementary Figure 1) in which the four CGEMS outliers are clearly illustrated.

Supplementary Material

[Population abbreviations as above in Supplementary Figure 1.] Supplementary Figure 3: Genotypic coverage achieved across the CGEMS cohort from original genotypes and sufficiently well-imputed genotypes. Each bar in total represents the total number of SNPs per chromosome for which imputation was performed.