Targeted next generation sequencing (NGS)

Targeted NGS sequencing was performed on the 199 patient samples, using an Illumina MiSeq [18] instrument. Initially 12 runs were performed to sequence the samples; later 15 samples were re-sequenced due to low quality of the initial sequencing. The target regions (all introns, exons, 5’ and 3’ UTRs) of 22 MMR genes (MLH1, MSH2, MSH6, PMS2, MSH3, PMS1, MLH3, EXO1, RFC1, RFC2, RFC3, RFC4, RFC5, PCNA, LIG1, RPA1, RPA2, RPA3, POLD1, POLD2, POLD3 and POLD4) with a total size of 1.213 Mb were captured using 6961 probes and the Illumina Nextera Rapid Capture Enrichment Kit (custom, 96 samples). An overview of these 22 genes, their function and associated phenotypes are shown in Table 1. Sequencing was performed at the Medical Genetics Laboratory at Hunter Medical Research Institute (HMRI), University of Newcastle, Australia.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. List of genes in target panel.

https://doi.org/10.1371/journal.pone.0235613.t001

Bioinformatic analysis

Raw reads (.fastq files) generated by the sequencer were processed by the following three major steps:

1. Data pre-processing: Raw reads were aligned to the reference genome (version hg19), and sequence alignment maps were generated. These alignment maps were used for read visualization and to call variants.
2. Variants discovery: The alignment maps generated from previous steps were compared against the reference genome to generate a list of nucleotide variants.
3. Variants annotation: Variants were annotated using different databases and tools.

A pipeline was constructed to perform the above-mentioned steps of analysis. Detailed overview of pipeline and tools used can be found as S1 File. Schematic overview of the pipeline is shown in Fig 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Schematic overview of the bioinformatics pipeline.

https://doi.org/10.1371/journal.pone.0235613.g001

Filtration of variants

All called variants were annotated by using Alamut-batch [19] before filtering. Filtus [20] was used for filtering variants. All variants were classified into 4 region-wise categories; exons, UTRs, introns, and splice sites (variant distance ≤ 10 nucleotides from nearest splice site). In the first stage of filtering, variants from all these four regions were filtered based on frequencies of variants in the gnomAD database [21]. Exonic variants, intronic variants, and variants near splice sites were filtered-in for frequencies less than 0.1% (or no frequency). UTR variants were filtered-in for frequencies less than 0.01% (or no frequency). In further stages of filtering, different strategies were adopted for every region. See Fig 2 for the workflow. Detailed filtering steps can be found in S2 File.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Filtering workflow and number of genetic variants detected.

https://doi.org/10.1371/journal.pone.0235613.g002

Validation of variants

Sanger sequencing was performed for validation of selected variants. The fragments were amplified using AmpliTaq Gold® 360 MasterMix and 360 GC Enhancer (Life Technologies). Cycle sequencing reaction was performed with BigDye® Terminator v3.1 (Life Technologies) and subsequent capillary electrophoresis was performed on the ABI 3130xl or ABI 3730 (Life Technologies). List of primer sequences can be provided upon request. Sanger sequencing data was analysed using SeqScape Software v3.0 (Life Technologies). Some variants have not been verified by Sanger sequencing, partly due to unavailability of primers for some of these gene, but also due to logistic issues. But variants were thoroughly inspected in BAM files to assure they were likely to be true positive variants (enough coverage and an allele fraction of about 50%, between 30 and 75%).

Interpretation and classification of DNA variants

The remaining variants after filtering were classified into 5 classes according to the American College of Medical genetics (ACMG) guidelines [22]. To determine whether these variants had been detected before, literature and databases including LOVD/InSIGHT (https://www.insight-group.org/variants/databases/) and ClinVar (https://www.ncbi.nlm.nih.gov/clinvar/) were searched. Potential pathogenicity of missense variants was interpreted using Alamut batch (annotation) [19] and Alamut Visual (interpretation) [23].

Exonic variants

A total of 207 variants were called in exonic regions of the target panel, over all samples. The variants were filtered by removing cases according to their frequency in gnomAD (> 0.1%) and annotation in ClinVar (benign/likely-benign) [49]. Of the 22 exonic variants (S3 File) that remained after filtering, there were 2 putative pathogenic variants, 7 variants of unknown significance (VUS) and 13 variants without any information (NO_Info) according to ClinVar (only non-synonymous variants).

Among these 22 variants there were 2 variants in MLH1 (NM_000249.3). Both MLH1 variants were classified as class 3 in pathogenicity, according to ACMG guidelines [22]. The variant c.453G>A p.(Thr151 =) is found in the last nucleotide of exon 5. It may alter the ligation of adjacent exons 5 and 6 and is predicted to be splice site deactivating by prediction tools (SSF [50], MES [51]) (nearest-SS-change score: -0.29). The first and the last three positions of the exon are an integral part of the 3’ and 5’splice site consensus sequences [52], the variant position is highly conserved (PhastCons score: 0.99), and predicted as pathogenic by UMD-predictor [53]. According to ClinVar it is classified as a likely pathogenic / VUS variant, with multiple submissions in ClinVar where many of them has a HNPCC/Lynch syndrome phenotype. With strong evidences for being a pathogenic variant, it is a candidate for further RNA/functional studies. The variant c.2009A>G p.(Lys670Arg) has no frequency in the gnomAD database, but has recently been reported in ClinVar (as VUS) and in other databases. This variant has been associated with a HNPCC phenotype and hereditary cancer-predisposing syndrome, according to ClinVar. The variant position is highly conserved (PhastCons: 1, phyloP: 4.6) and lies in a helix secondary structure of the protein. It has been predicted as pathogenic (UMD-prediction, MutationTaster).

There were also four exonic variants in MSH6 (NM_000179.2). Two of the variants, c.335A>G p.(Asn112Ser) and c.2203C>A (p.Leu735Ile), have been classified as VUS by ClinVar. These two variants have previously been associated with Lynch syndrome, HNPCC and hereditary cancer-predisposing syndrome-like phenotype according to ClinVar and other databases. Both variants are at highly conserved positions and both have been predicted as damaging by prediction tools. We classified these two variants as class 3. Variants c.1409C>G p.(Ser470*) and c.3802-4_3825dup p.(Glu1276*) both code for “STOP gain” and are disease causing. None of these variants have entries in ClinVar or frequency in gnomAD. We classify these as class 5 variants.

Three exonic variants in MSH2 was identified (NM_000251.2), c.97A>C p.(Thr33Pro), c.1228G>T p.(Gly410Cys) and c.2732T>G p.(Leu911Arg), with all three classified as VUS by ClinVar. All three have phenotypic association to Lynch syndrome/HNPCC and hereditary cancer-predisposing syndrome and have been predicted as pathogenic/disease-causing by many prediction tools (UMD-prediction, PolyPhen, SIFT and MutationTaster). All three variant positions are highly conserved (with high scores in PhastCons and phyloP). Variant c.97A>C p.(Thr33Pro) was identified from a low quality sample (coverage depth at variant position 4X and sample coverage 7X), but was verified as a true variant by Sanger sequencing. It has been scored with a high value for decreasing protein stability (SNPs3D [54] score: -1.08) and has been suggested as a cause of reduced mismatch binding/release efficiency compared to wild-type protein in previous studies by Ollila et al. [26, 27]. Variant c.1228G>T p.(Gly410Cys) has no frequency in gnomAD, but has been reported to ClinVar and is located in a helix secondary structure of the protein. It has a high score for structural change (Grantham-Distance: 159), but it is predicted not to alter protein stability (SNPs3D: +3.43). Variant c.2732T>G p.(Leu911Arg) was also identified in a low quality sample (coverage depth at variant position 5X, sample 10X). It lies in a helix secondary structure, has high score for structural change (Grantham-Distance: 102) and for decreased protein stability (SNPs3D: -1.08). These three MSH2 variants have been classified as class 3.

A missense exonic variant in PMS2 was also detected (NM_001322014.), c.137G>T p.(Ser46Ile), which was found in two samples. It was classified as likely pathogenic according to ClinVar, and is reported to be a founder mutation [55]. The protein region has helix-like secondary structure (UniProt [56]), and the position is highly conserved (PhastCons [57] score:1; phyloP [58] score: 6.178). It has been classified as pathogenic by several prediction tools (UMD-predictor, PolyPhen [59], SIFT [60], and MutationTaster [61], the variant has been referred to in many previous studies, and it has been considered for strongly decreased DNA mismatch repair activity. This variant was classified as class 4.

One exonic variant in POLD1 was identified (NM_001308632.1), c.952G>A p.(Glu318Lys), was classified as VUS according to ClinVar. It was called in a low-quality sample (coverage depth at variant position 11X, sample’s mean coverage depth 22X). The position is highly conserved (PhastCons: 1, phyloP:3.9). The variant is in the DNA binding cleft of the exonuclease active domain of POLD1, it has a high score for decreased protein stability (SNPs3D: -2.68), and is predicted as damaging by prediction tools. A previous study has predicted it to be disease causing [62]. However, functional studies are needed to confirm pathogenicity, and therefore it was classified as class 3.

Exonic variants were also found in five other genes; MSH3, LIG1, RFC1, EXO1 and RPA3. All these variants were classified as class 3. In MSH3 (NM_002439.4) variants were c.1394A>G p.(Tyr465Cys), c.845C>T p.(Thr282Ile), c.1232G>A p.(Arg411His), c.1396A>G p.(Ser466Gly), c.2041C>T p.(Pro681Ser). In LIG1 (NM_000234.2) variants were c.1159C>T p.(Arg387Cys) and c.692T>G p.(Phe231Cys). In RFC1 (NM_001204747.1) this was c.3445T>A p.(*1149Argext*15), which introduces a “STOP loss” and extension of 15 amino acids in the product protein and c.2017G>T p.(Val673Leu); in EXO1 (NM_006027.4) c.409G>T p.(Ala137Ser); In RPA3 (NM_002947.4) it was c.295T>C p.(Tyr99His).

Intronic variants

Among all detected variants, 9,260 were identified as intronic, which was ~97% of all variants. Intronic regions of human DNA, being extraordinarily larger in comparison to other regions, it is expected to find most of the variants in these non-coding regions. After frequency-based filtering (< 0.1%), this list was reduced to 4,197 variants, which was further reduced by splice site related filtering, using strict filtering criteria to reduce the large number of variants. These variants were filtered for two categories, first for “New Donor/Acceptor site” and then for “Cryptic Donor/Acceptor Site STRONG activation” (see S2 File for filtering details). We found in total five variants, with four variants in the first category and one in the second (see S3 File). Two of these variants were in RPA1 and RPA3, and have been predicted as new acceptor site, two were in RFC1 and RFC4 and have been predicted as new donor sites, and one was in RPA3 and has been predicted to give strong activation of a cryptic donor site.

According to the Human Gene Mutation Database (HGMD) more than 10% of all disease-causing hereditary mutations are splice site altering [63–65]. Variants in vicinity of exon-intron junctions were therefore studied. After filtering (see Supporting Material), we found four variants of interest in the vicinity of splice sites (see S3 File). Two of these were in MSH6 (NM_000179.2). The variant c.628-7C>A has been classified as VUS by us and ClinVar, The variant c.3439-1G>T, at the last nucleotide of 5^th intron, has been classified as pathogenic by ClinVar, it has been linked to LS/HNPCC phenotype, and has a maximum score (-1) for splice-site deactivation. We classified it as class 5 and hence disease causing.

An intronic variant in MLH1 (NM_000249.3: c.306+4A>G) is found close to a splicing junction and was predicted for splice site deactivation. It is in a highly conserved position (PhastCons:1, phyloP:4.2), and has been classified as VUS in ClinVar. Experimental studies have shown that this variant results in the activation of a cryptic donor site and skipping of exon 3 in an ex-vivo splicing minigene assay [24], but as no studies have verified this in patient samples, we classified it as class 3 variant. An intronic variant in MSH2 (NM_000251.2:c.1277-7C>A), previously classified as likely benign, we classified as a class 3 variant.

UTR variants

There were 140 variants identified in UTR regions. Due to limitations of annotation tools and databases, any effects of most mutations in these regions are hard to predict. Hence, a relatively strict filtering compared to standard (for diagnostics) [66] was used for variants in these regions, to reduce the number of variants to a manageable size. After frequency-based filtering (< 0.01%) this list reduced to 28 variants, of which 9 variants were in 5’ and 19 variants were in 3’ UTR.

Variants in 5’ UTR were annotated for transcription factor binding sites (TFBSs), using the UniBind database [67]. Among the 9 variants in 5’ UTR, three had significant hits in the database, where each of these three variants was found inside a potential binding site for at least one transcription factor (TF) according to UniBind data (see S2 Table). One variant in MSH6 (NM_000179.2: c.-76G>T) had overlap with potential binding sites for the TFs CTCF, STAT3, E2F7 and E2F1. For CTCF there is a high frequency of the reference allele (G) compared to the alternate allele (T) at the variant position, which can indicate a strong preference for the reference variant, and possibly a significant effect of the alternate variant on TFBS specificity (frequency matrices from the JASPAR database [68, 69] were used for this analysis). According to ChIP-seq data visualized with the UCSC genome browser [70] there are relatively strong signals for CTCF at this position (see S1 Fig) compared to other potential TFs. Mutations in CTCF binding sites have for example been associated with chromosomal instability and aberration and have been found in gastric and colorectal cancer [71], which strengthens the possibility that this variant may have an effect through altered binding of CTCF. A variant in RFC3 (NM_002915.3: c.-106A>G) had hits for the TFs GABPA, JUN, CREM, JUND, ATF1, MITF, NR3C1, ATF7 and CREB1. Among these hits, 6 TFs (JUN, CREM, JUND, R3C1, ATF7 and CREB1) had a very high frequency of the reference allele (A) compared to the alternate allele (G) at the variant position. ChIP-seq data shows strong signals for CREB1 (see S2 Fig.), which may indicate a potential for significant effects due to alteration in the binding site. A variant in RFC4 (NM_002916.3: c.-90C>T) had a hit for the TF AR.

Nineteen variants in 3’UTR were annotated using TargetScan v6.2 [72] and a two-step SVM prediction of micro-RNA (miRNA) target sites [73]. A SVM score normalization method [74] was used to normalize the score and miRNA data were taken from MirBase v22 [75]. Only a variant in gene EXO1 (NM_130398.3:c.*84T>G) was predicted as a likely true candidate for affecting miRNA binding, for miR-370-3p and miR-93-3p (see S3 Table). Several studies have shown the importance of EXO1 in replication, DNA repair pathways, cell cycle checkpoints and its association to cancer [76], and GWAS studies have identified specific mutations in EXO1 gene as risk alleles for different types of cancer [77, 78]. SNPs in miRNA binding sites have been associated with CRC [79]. For the two miRNAs predicted to be affected by variation in their binding site, miR-370-3p has been identified as a tumour suppressor in EC via endoglin regulation [80]. The miR-93-3p can be considered as an important factor for CRC suppression and inhibition of tumorigenesis [81], as a previous study has associated the down-regulation of miR-93 with unfavourable clinicopathologic features and short overall survival of CRC patients [82].

Implications of the study

In this study we found 35 significant variants (22 exonic, 4 UTR, 9 intronic), with 15 variants in the 4 MMR genes known to cause LS (MLH1, MSH2, MSH6, PMS2) and 20 in the additional MMR genes included in this study (MSH3, POLD1, RFC1, RFC3, RFC4, LIG1, EXO1, RPA1, RPA3). This helped in identification of variants in less studied genes, as well as polygenic variations (although none of the 199 samples in this particular study showed polygenic variants of interest for further investigations). This study also used the complete genomic regions of the genes, which very few previous studies have done [83, 84].

Though all known genes in the MMR pathway were studied, there will always be a possibility of additional genes and associated variants with similar disease effects, e.g., POLE mutations in EC cases contributing towards Polymerase Proofreading Associated Polyposis (PPAP) [85, 86], or germline deletions in another gene (EPCAM1) leading to silencing of the MSH2 gene, causing Lynch syndrome [87]. These limitations can only be removed by expanding the panel by including more genes, up to the extent of the whole genome. However, this will also increase the potential of noise and complexity of the analysis, by including more genes and variants that are less likely to be relevant in a given study. Another limitation is associated with the PMS2 gene in this panel, which has a pseudo-gene (PMS2CL), and where 6 exons (exon 9, 11, 12, 13, 14 & 15) are highly similar to PMS2CL. This creates challenges in alignment of correct reads at these exons and creates artefacts during variant calling. This limitation has also been mentioned in a pilot study [83]. This makes it important to manually check reads and coverage in a genomic viewer, and to do Sanger verification of variants, as we did for the PMS2 gene.

The current study emphasises the importance of including non-coding intronic regions. These regions will often have splice site variants, which may contribute to 10% of all disease-causing hereditary mutations according to HGMD [63–65], and deep intronic variants (e.g., in branch-point sequences, U2 type introns) which also contribute towards disease, most frequently by creating new pseudo-exons by activating non-trivial splice sites or by changing splicing regulatory elements. Intronic variants can also disrupt transcription regulatory motifs and non-coding RNA genes [88]. However, it is challenging to annotate these intronic variants due to limitations of annotation databases and tools. In a clinical setting, these variants can easily be missed unless RNA studies are performed to check for exon skipping, generation of new donor sites or cryptic site activation. Considering the potential importance of such variants, the current study included all intronic regions in order to search for this type of variant. Among 10 significant intronic variants we found four in the splice site vicinity and six in deep intronic regions.

NGS was performed, aiming at a data quality greater than 100X (average read coverage depth) for all samples. However, only 23 samples achieved this coverage (highest among them 169X), whereas 50 samples had coverage of less than 30X. These low-quality samples were included in the study, with the aim of exploring the value of low-quality data when searching for true positive (TP) variants. Using low quality data (i.e., with low coverage) led to a higher fraction of false positive (FP) variants, as 16 variants identified from the data analysis were subsequently identified as false positive variants by Sanger sequencing. Most of them had low coverage at the variant position (between 14X to 6X coverage), whereas others were in repeat regions. FP variants in MLH1, MSH6 & PMS2 genes were in repeat regions, and had low coverage (except a PMS2 variant with 84X coverage), which possibly led to their false SNV call. On the other hand, we cannot rule out that some were not verified due to SNPs in the primer binding site (allelic dropout). However, we also found many true positive variants in low coverage regions, as we found and confirmed 6 true positive variants in regions with low coverage (between 16X and 4X). Among these were two class 3 variants (MSH2), one class 3+ variant (POLD1) and a class 5 variant (MSH6). This shows the potential for finding true variants of significance even in low-quality samples, given that the variants can be verified.

Our initial aim also included identification of CNVs. CNVs can occur in both exonic and intronic regions of protein coding genes, with intronic CNVs being more frequent [89], and both types can contribute towards disease. However, due to the limitations of data quality (non-uniform and low coverage depth), it was not possible to do reliable CNV calling. Also, there is no availability of MLPA kits (MRC-Holland) for detecting CNVs for many genes in this panel.

To associate variants with possible effects we utilized in silico resources and tools, in addition to published literature. Effect prediction and annotation of all variants was done using multiple tools as mentioned in the methods section. Also, multiple potential factors and effects, like conservation in variant position or structural changes at protein level, were checked for each variant. This consensus-like approach (multiple tools, multiple potential effects) increases the robustness of predictions and annotations of the variants, although we also had cases of contradictory predictions, which illustrates the challenge of using in silico prediction tools.

Among the 199 EC patient samples, we identified variants of interest (for further investigations) in 34 patients. Among these, we found 3 patients with class 5 variants (in MSH6 gene) and two patients with the same class 4 variant (in PMS2 gene); ~2,5% of patients had pathogenic variants representing a very likely cause of cancer in these five patients. This is in accordance with other studies. One meta-analysis of 53 studies concluded the prevalence of LS in EC patients to be approximately 3% [90]. These studies have only looked into the coding part of the four MMR genes MLH1, MSH2, MSH6 and PMS2. We found class 3 variants in another 29 patients, some of which are highly suspicious of being pathogenic variants. This indicates potential causes of their disease, although further studies are required to confirm their actual significance. It is an important limitation for further interpretation of these class 3 variants that we lack information about the patients debut age of cancer, and results from tumour analyses (MSI status and immunohistochemistry of MMR genes). For the remaining 164 patients, we did not find any significant variant to explain their disease. Expansion of panel size with more genes, improved annotation (particularly of variants in non-protein-coding regions), and improved data quality may help in explanation of some of these cases. However, since the study cohort consist of consecutive EC patients, most of the cancers will be sporadic with no underlying high penetrant genetic cause.