Inherited Predisposition of Lung Cancer: A Hierarchical Modeling Approach to DNA Repair and Cell Cycle Control Pathways

The familial relative risk of lung cancer was reported to be ∼2-fold from several registry-based studies (1-3). Although shared environmental factors may explain part of the elevated familial relative risk, a meta-analysis reported a 1.5-fold elevated risk of lung cancer among never-smoking probands with affected first-degree relatives (4), which provides evidence for a genetic component in lung carcinogenesis. Despite a recent linkage analysis of 52 high-risk pedigrees that localized a lung cancer susceptibility locus at chromosome 6q23-25 (5), the exact inheritance mechanism of lung cancer is still largely undefined. The candidate gene approach has guided research on inherited susceptibility to lung cancer in the past decade. However, most of the associations in such studies have not been replicated. The likely reasons for this include lack of study power and false positives due to multiple comparisons, which may be exacerbated by publication bias. We therefore conducted a large-scale genetic association study on lung cancer to ensure adequate statistical power and conducted analyses using hierarchical modeling to include the prior knowledge of sequence variation into the modeling. The present investigation is focused on the sequence variations in DNA repair and cell cycle control pathways, two critical defense mechanisms against human carcinogenesis.

The DNA repair system maintains the integrity of the human genome by removing DNA damage, reducing the mutation frequency of cancer-related genes, minimizing replication errors, and minimizing deleterious rearrangements arising via aberrant recombination (6). DNA repair processes are classified into four main pathways responsible for repairing different classes of DNA damage. First, base excision repair occurs through the excision of a modified base, followed by fill-in repair synthesis using the opposite strand as a template. Second, nucleotide excision repair removes photoproducts formed by UV radiation and other bulky adducts formed by a variety of chemicals. Third, DNA double-strand breaks are repaired either through the nonhomologous end-joining pathway or by homologous recombination. Finally, mismatch repair operates on base mismatches and small loop-outs that arise during replication by misincorporation or slippage on the template strand (7) and direct reversal of methylation is operated by enzymes such as O₆-methylguanine-DNA methyltransferase (MGMT). Cells with damaged DNA must either pause in the cell cycle to allow for repair or succumb to elimination by apoptosis and the activation of cell cycle checkpoints is a critical component of the cellular response to DNA damage (8).

We hypothesized that genetic variants in DNA repair and cell cycle control pathways will influence the predisposition to lung cancer, and studied 27 variants in 17 DNA repair enzymes and 10 variants in eight cell cycle control genes in 1,604 lung cancer patients and 2,053 controls. The majority of variants were selected from published association studies or chosen from the SNP500 project¹

To improve the estimation of risks for specific variants, we applied a Bayesian approach in which we allowed the prior knowledge to be incorporated in a hierarchical model. To develop a coherent scale of prior probability of whether the variant is likely to be deleterious or neutral, we took an evolutionary perspective, which takes account of evolutionary biology and physicochemical properties of the variant. The underlying idea is that the ability of a missense substitution to alter the function of the protein is related to two factors: (a) the degree of sequence conservation observed at the position of the missense substitution in the candidate genes and (b) the chemical difference between the missense substitution and the canonical residue at its corresponding position in the protein (10, 11). Such an approach has been implemented in the classification algorithm SIFT (Sorting Intolerant from Tolerant),²

The study was conducted in 15 centers in six countries of Central and Eastern Europe, including Czech Republic (Prague, Olomouc, Brno), Hungary (Borsod, Heves, Szabolcs, Szolnok, Budapest), Poland (Warsaw, Lodz), Romania (Bucharest), Russia (Moscow), and Slovakia (Banska Bystrica, Bratislava, Nitra). The study details have been previously described (12). Briefly, each center followed an identical protocol and recruited a consecutive group of newly diagnosed cases of lung cancer and a comparable group of controls between 1998 and 2002. Eligible subjects were resident in the study area for at least 1 year. All lung cancer diagnoses had a histologic or cytologic confirmation. Consent for participation was required from both the patient and the patient's physician. This study was approved by the institutions at all study centers, and ethical approval was obtained from the IARC (Lyon, France), the coordinating center.

Controls in all centers except Warsaw were chosen among subjects admitted as in-patients or out-patients in the same hospital as the cases, and were frequency matched with the case group by sex, age (±3 years), center, and referral (or residence) area. The eligible diseases for controls were non–tobacco-related diseases, including minor surgical conditions, benign disorders, common infections, eye conditions (except cataract or diabetic retinopathy), common orthopedic diseases (except osteoporosis), etc. In Warsaw, population controls were selected by random sampling from the Polish Electronic List of Residents.

Both cases and controls underwent an identical face-to-face interview during which they completed a detailed questionnaire including sections on (a) demographic variables, (b) medical history, (c) family history of cancer, (d) tobacco smoking and involuntary smoking, (e) alcohol drinking, (f) consumption of selected food items, and (g) occupational history.

Blood samples were collected at the time of interview. In total, 2,250 cases and 2,289 controls have provided blood samples. Genomic DNA was extracted from blood samples by automated Gentra Systems; DNA concentrations were measured by using PicoGreen dsDNA Quantification kits (Molecular Probes). Genomic DNA was extracted from blood samples with the use of a QIAamp 96 DNA Blood Kit (Qiagen), or with Puregene chemistry (Gentra Systems) on an Autopure instrument (Gentra Systems). Samples that yielded an insufficient amount of DNA at extraction were subjected to whole-genome amplification by use of a phi29-based protocol (GenomiPhi, Amersham Biosciences) or reextracted using Gentra technology. DNA concentrations were measured by using PicoGreen dsDNA quantification kits (Molecular Probes).

A total of 37 single nucleotide polymorphisms (SNP) were genotyped, including APEX Q51H, D148E, LIG3 Ex21-250C>T, OGG1 S326C, PARP V762A, XRCC1, R194W, R280H, R399Q, LIG1 7C>T, ivs 9-21A>G, ERCC1 8092C>A, 354T>C, XPA 23G>A, XPC K939Q, XPD D312N, K751Q, XPF P379S, R415Q, XPG D1104H, MGMT 171C>T, L84F, I143V, K178R, MLH1 I219V, XRCC2 R188H, XRCC3 T241M, XRCC4 S298N, CHEK2 I157T, CCND1 G870A, MDM2 162A>G, MDM2 309T>G, MMP1 insG, p16/CDKN2A A148T, p21/CDKN1A S31R, TP53 R72P, intron 3 16-bp repeats and P73 14C>T (rs numbers are listed in Table 1). Thirty-three of them were analyzed by the 5 exonuclease assay (i.e., TaqMan assay), one using the allele-specific PCR-based Amplifluor assay (XRCC1 R399Q rs25487), and two using the MGB Eclipse Probe System (P73 C14T, rs1801173: MLH1 I219V, rs179977). In addition, the presence of the TP53 intron 3 16-bp duplication was assessed using a gel-based discrimination assay. Designs of genotyping assays for the SNPs were either taken from the Web site of the SNP500 project¹ or were designed by Applied Biosystems, Proligo, or in-house using the Genbank databases and Primer Express (Applied Biosystems). PCR primers and TaqMan probes were synthesized by Applied Biosystems or Proligo. Eclipse probes were obtained from Nanogen and Ampliflour probes were from Serologicals Corporation. Briefly, 10 ng aliquots of DNA were placed into separate wells of a 96-well or 384-well PCR plate along with a PCR cocktail that included fluorescently labeled allele-specific probes. The fluorescence of the PCR products was then plotted and genotype was determined according to the signal of the two probes. Details of primers, probes, and protocols for all genotyping are available upon request. Ten percent of the subjects were randomly selected and regenotyped for each polymorphism to examine the reliability of the genotyping. In total, 1,604 cases (71.3% of all cases) and 2,053 (70.8% of all controls) had complete genetic data for all makers in DNA repair and cell cycle control pathways, and were included in the final hierarchical regression. Primary data of 9 of these 37 variants included in the hierarchical regression have been previously reported, including variants of CHEK2, XRCC1, OGG1, TP53, CCND1, and p16/CDKN2A (12-14).

The frequency distribution of demographic variables and putative risk factors of lung cancer, including country of residence, age, sex, education, and smoking was examined for cases and controls. Tobacco smoking included cigarettes, pipes, and/or cigars. Former smokers were defined as smokers who quit smoking at least 2 years before interview or diagnosis. Cumulative tobacco consumption was calculated as the product of smoking duration and intensity, and expressed as pack-years. Hardy-Weinberg equilibrium of allele distributions was tested in cases and in controls separately.

Genotypes were dichotomized into two categories based on the prior knowledge or frequency distribution of the SNPs as risk versus nonrisk genotypes. Minor alleles were assumed to confer an increased risk unless there was prior evidence indicating otherwise. For markers in linkage disequilibrium, we conducted haplotype analyses based on the expectation-maximization algorithms using the TagSNPs software (15). When there was no evidence of a specific haplotype conferring the primary association with lung cancer risk, we considered subjects who carried any risk genotype for the markers in linkage disequilibrium as the subjects at risk. On the other hand, when there was evidence that a specific haplotype conferred the primary association with lung cancer risk, we used the haplotype dosage as the input variable in the first stage. For example, two variants in ligase 1 (7C>T and inv9-21A>G) are in linkage disequilibrium (D′ = 0.83), and haplotype analysis showed that 7T-inv9G was the primary haplotype associated with lung cancer risk [odds ratio (OR), 1.19; 95% confidence interval (95% CI), 1.02-1.38]. In this case, the haplotype dosage of 7T-inv9G was used as the input value for ligase 1 risk genotype. Similar results were observed for P53 R72P and intron 3 16-bp repeats.

The hierarchical model we used has been described previously (16). Briefly, we applied a two-stage model in which we used the conventional logistic regression as the first stage, and added a second stage with the goal of improving the estimation of the coefficients. The prior matrix in the second-stage model was built on our prior knowledge of the biological function of the gene and the possible functional significance of the variants. The genetic variants were broadly classified into five pathways related to lung carcinogenesis, based on the biological functions of the genes: (a) base excision repair, (b) nucleotide excision repair, (c) mismatch repair and direct reversal, (d) double strand break, and (e) cell cycle control genes. For each pathway, a score was assigned to each genetic marker according to the biological function of the gene and the functional significance of its polymorphism. A score of 0 was assigned to the genetic polymorphism if the gene is not involved in such a pathway (the coefficient of the factor was expected to be zero for such a pathway and the corresponding gene effect was not shrunk toward the other genes in that pathway either). Genes assigned to the same pathway were assumed to be exchangeable; that is, for their effects to arise from a common distribution. The five pathways were not mutually exclusive.

In addition to its role in different pathways, whether a nucleotide change is likely to be deleterious would also influence the prior probability of how likely that the variant would have an effect on cancer risk. One method to predict whether a nucleotide change is likely to be functionally important is to measure how conserved is the region where the nucleotide change occurs: a high conservation level of a genetic region across species indicates that the region is likely to be fundamentally important. Thus, amino acid changes that occur in a more conserved genetic region would have a higher likelihood to be functionally important. We have subsequently expanded the hierarchical modeling by incorporating information on sequence conservation of the variants into the prior model (second stage). Sequence conservation of the variants were obtained by SIFT,² which searches for the sequences of the same gene across species and returns the probability of such an amino acid appearing at the specific position (denoted as P_ca). This probability is used as an indication of how tolerant an amino acid substitution is at the given position (11, 17). For each gene variant, the importance of its functional changes is weighted by the probability that these amino acid substitutions are intolerant or deleterious (which is 1 − P_ca, denoted as D_ca). This estimation can only be obtained for variants that lead to an amino acid substitution; other variants were assigned a D_ca of 0.

The prior matrix used in the second-stage model is reported in Table 1. In summary, the biological implication of this prior matrix is that APEX1, LIG3, OGG1, PARP, XRCC1, and LIG1 play a role in base excision repair pathway, and, LIG1, ERCC1, XPA, XPC, XPD, XPF, and XPG are involved in the nucleotide excision repair pathway, whereas MGMT and MLH1 encode enzymes that are part of the direct reversal or mismatch repair pathway, and XRCC2, XRCC3, XRCC4, and CHEK2 are part of the double strand break pathway or downstream from it. Their designated risk genotypes or haplotypes decrease the DNA repair capacity through various mechanisms. Furthermore, CHEK2, CCND1, MDM2, MMP1, CDKN2A, CDKN1A, TP53, and P73 are crucial in the cell cycle control pathway, and their designated risk genotype or haplotype are hypothesized to lead to aberrant cell cycle control. As described above, D_ca served as an additional covariate in the prior matrix for nucleotide changes that lead to amino acid changes. For markers with multiple genotype or haplotypes, we used the average of D_ca as the covariates.

We applied the hierarchical modeling with the empirical-Bayes approach (allowing τ², the variance of the residual effect to be estimated from the data). When the estimate of τ² was set to 0, which may occur due to model overparameterization, we subsequently took the semi-Bayes approach and specified τ² as 0.05 (which would allow ∼2.4-fold of variation in the effect of the gene variants).

We conducted stratified analyses by histology to investigate the effect of each genetic polymorphism on each major histologic type. We estimated the effect of sequence variants separately for subjects with no family history of any cancer, with family history of lung or upper aerodigestive tract cancers, or with family history of lung cancer. The rationale for analyzing cases with family history separately is that one would expect genetic factors to have a stronger effect on cancer development among such subjects. We also evaluated the modulating effects of tobacco consumption, by stratifying and comparing the stratum-specific risk estimates.

Matching variables and potential confounders, such as country of residence, age (continuous), sex, and smoking pack-year (continuous), were included in the multivariate logistic regression at the first-stage model. Hierarchical modeling analyses were conducted with SAS IML code in conjunction with GLIMMIX macro (18).

The frequency distribution of the demographic factors and the putative risk factors among subjects with complete genetic data are shown in Table 2. Comparing cases and controls on demographic variables, controls were slightly more educated than cases. As expected, smoking prevalence was higher among cases. The percentages of subjects from each country ranged from 6% to 33%, but the distributions were slightly different between cases and controls (P < 0.01) due to differences in DNA availability from different countries. None of the SNP allele frequencies in the control group deviated from the Hardy-Weinberg equilibrium at P < 0.01.

Table 3 shows the main effect of each genetic variant, as well as the estimated effects of each pathway. For the sake of comparison, we presented five estimates for the main effects from (a) a conventional analysis with single marker in each logistic regression, (b) conventional analysis with all markers in one single logistic regression, (c) hierarchical modeling with empirical Bayes approach and pathway indicators only, (d) hierarchical modeling with the empirical Bayes approach and a SIFT covariate D_ca, and (e) hierarchal modeling with semi-Bayes approach and a SIFT covariate D_ca. The single marker analysis suggested that five genes were associated with lung cancer risk (OGG1, LIG1, ERCC1, CHEK2, and TP53), among which two risk estimates were no longer significant when we included all markers at once. The conventional analysis with all genes in one model suggested that four genes (OGG1, MGMT, CHEK2, and TP53) might be associated with lung cancer risk, three of which were supported by the empirical-Bayes hierarchical modeling (OGG1, MGMT, and CHEK2). When we allowed the risk to vary to a larger extent using a prespecified τ² of 0.05 (i.e., 2.4-fold variation), OGG1, MGMT, and CHEK2 remained associated with lung cancer risk. Based on the estimation from the hierarchical modeling, subjects who carried OGG1 326C/326C genotype, MGMT 143V or 178R, and CHEK2 157I allele had ORs of lung cancer of 1.45 (95% CI, 1.05-2.00), 1.18 (95% CI, 1.01-1.40), and 1.58 (95% CI, 1.14-2.17), respectively. The association of CHEK2 157I allele, which has been previously reported and independently replicated (14), seems to be robust in the hierarchical modeling. Although TP53 A2-72P was shown to be associated with increased risk of lung cancer in the conventional analysis, the association was not supported by the hierarchal modeling. Based on the data we have for each pathway, we could estimate pathway effects, although none appeared to be prominent. Stratification by family history of cancer did not seem to be informative, which might be due to small numbers.

Table 4 shows the effect of genetic markers on different histologic subtypes. A semi-Bayes approach was adopted for adenocarcinoma, as the τ² was estimated to be zero in the empirical-Bayes approach. CHEK2 157I allele seemed to have a stronger effect on squamous cell carcinoma with OR of 2.00 (95% CI, 1.16-3.47). On the other hand, the XRCC2 188H allele and XRCC3 214T allele seemed to be more important for small cell carcinoma with an OR of 1.61 (95% CI, 1.15-2.24), and 1.59 (95% CI, 1.10-2.28), respectively. OGG1 seemed to be positively associated with adenocarcinoma (OR, 1.54; 95% CI, 1.03-2.30).

In the analysis stratified by smoking status, we observed a positive association between lung cancer risk and XRCC2 188H and CHEK2 157I among current smokers, with ORs of 1.30 (95% CI, 1.03-1.65) and 1.50 (95% CI, 1.07-2.11), respectively (see Table 5).

In this study, we applied hierarchical modeling in a Bayesian framework to incorporate evolutionary biology into the risk estimates. In general, our study supports the notion that sequence variation may explain at least some of the variation of inherited susceptibility. In particular, sequence variants in CHEK2, MGMT, and OGG1 were suggested to be positively associated with lung cancer risk overall.

MGMT is a DNA repair protein that acts on O⁶-methylguanine, one of the most potent premutagenic lesions (19), by transferring the methyl group from the O⁶-methylguanine moieties to a cysteine residue located at the active site of the MGMT protein. The expression level of MGMT in humans varies up to 40-fold (20), and studies have shown that genetic components might explain part of the interindividual variation (21). Several studies have been conducted to investigate the association between MGMT polymorphisms and lung cancer risk, but the results are inconsistent, mainly due to lack of statistical power, because the majority of the studies had under 400 case-control pairs (22-27). I143V and K178R polymorphisms are in strong linkage disequilibrium. K178R was suggested to be unlikely to lead to altered AGT activity, as human AGT can be truncated from the COOH-terminal to the position 176 without loss of activity (28); nevertheless, the I143V variant was suggested to be functionally important as it is located adjacent to the alkyl acceptor ¹⁴⁵Cys of the active site (29) and was therefore hypothesized to modulate its repair activity and increase cancer risk. The association we observed with lung cancer might reflect the functional consequence of I143V, although further confirmation through experimental studies would be required.

The primary analysis of CHEK2 I157T and OGG1 S326C polymorphisms was reported previously (12, 30), in which we identified CHEK2 157T to be a protective allele and OGG1 326C to be a risk allele. The risk estimate using hierarchical modeling confirmed the positive associations reported in the previous publication, for both the overall main effect and the modifying effect of histology subgroup. Specifically, the prominent associations of CHEK2 I157T variant with squamous cell carcinoma, and of OGG1 S326C with adenocarcinoma were both supported by the risk estimates from hierarchical modeling. However, we did observe a higher point estimate from the conventional model compared with the estimates derived from the hierarchical model. This can be due to either an overestimation of the conventional model or inaccurate assumption of the prior distribution in the hierarchical model. In addition, the precision of the estimates from the hierarchical modeling was substantially higher than in the conventional model: For example, the OGG1 326C variant conferred an OR of 2.24 (95% CI, 1.30-3.86) for adenocarcinoma under conventional regression model, whereas the hierarchical model obtained a more precise estimate of 1.54 (95% CI, 1.03-2.30). The biological implications of the association between CHEK2, OGG1 and lung cancer risk were discussed in the previous reports.

XRCC2 188H was shown to increase risk of lung cancer among current smokers, and the effect was mainly on small cell carcinoma. XRCC2 protein is an essential component of homologous recombination involved in the double-strand breaks. In vitro studies showed that cells with a mutation or deletion of amino acid at codon 188 had a lower survival fraction after mitomycin C–inducted DNA damage (31). The common polymorphism R188H was also shown to have a small effect on XRCC2 function and cell survival after DNA damage (31). Only limited studies have been conducted to assess the association between XRCC2 R188H and lung cancer risk and the results have been inconsistent (26, 32, 33). Our observations of XRCC2 188H allele among current smokers and small-cell carcinoma are compatible to the previous experimental evidence.

TP53 variants, on the other hand, were suggested to increase lung cancer risk in conventional analysis, but the risk estimate was no longer significant in the hierarchical model. This could indicate that the results from the conventional analyses were overly estimated, or that the effect of TP53 variants may be derived from a different distribution compared with the other variants in the cell cycle pathways. To further clarify the role of TP53 variants and lung cancer risk, a more detailed analysis of TP53 based on haplotype tagging SNPs is currently under way.

Recent genetic association studies have been hampered by failures to replicate results, which could be due to multiple comparison or selective reporting. We showed that hierarchical modeling may improve the analyses of genetic data by increasing the precision of the risk estimates (e.g., as evidenced by narrower confidence intervals), and reducing the likelihood of false positive via shrinkage toward to the prior mean. The hierarchical modeling approach also has an appealing advantage of providing the risk estimate of each pathway investigated, although none of the pathways investigated have seemed to be have a more prominent effect than the others. This might be due to the low number of markers included in each pathway.

There are several limitations of our study. The first limitation to be considered is the possible violation of the exchangeability assumption, as the true effect size of the variants in the same pathway might not arise from the same distribution. Some variants in the pathway may have a larger effect, whereas others in the same pathway have very modest effects; or, in a more extreme scenario, the alleles specified in the same pathway may confer opposite effects. In addition, we dichotomized the genotypes into risk and nonrisk genotypes to minimize the number parameters to be modeled. Given that the previous literature on some of these variants remains inconclusive, the assumption of the minor alleles being the risk allele, unless there is prior evidence indicating otherwise, could also potentially lead to the violation of exchangeability. When the assumption of exchangeability is violated, the effect estimation would be shrunk toward a wrong prior mean. Nevertheless, previous sensitivity analyses have shown that hierarchical modeling is relatively robust to the alteration of the prior, and it should provide a more precise risk estimates (compared with the conventional model) under a reasonable specified prior (16, 34, 35). Second, despite our attempts to incorporate the evolutionary biology perspectives into our models in addition to the pathway indicators, we can only obtain prior information on nonsynonymous variants systematically based on the available tool such as SIFT. However, SIFT scores vary depending on depth and the species for which sequence alignment of specific variants is available. Thus, it is important to recognize that these scores are not the absolute measurements of how conserved the variants of interests are. Nevertheless, it provides reasonable representation of the evolution biology and physicochemical properties of the variants for the purpose of our analysis. Finally, we consider it reasonable to assume that the prior probability of a synonymous variant conferring a true effect to be very low. Nevertheless, it is possible that this type of variants can be functional through other mechanisms such as alternation of the splicing site, which currently is difficult to predict. The prior matrix can be further refined when such information become available.

In summary, our data provided robust positive signals of genetic regions that might harbor the inherited predisposition of lung cancer through a Bayesian framework and conservative estimation. To confirm the association found in this study, independent replication is required. Further investigation on OGG1, MGMT, and CHEK2 focusing on the genetic regions where the present markers are located or the haplotype blocks tightly linked with these markers might be warranted.

We thank Valérie Gaborieu for the technical assistance.