Introduction
Allogenic hematopoietic stem cell transplantation (HSCT) is the predominant treatment used to cure malignant and non-malignant hematological disorders. The number of HSCTs conducted has increased due to an overall improvement in the safety of the procedure resulting from reduced-intensity conditioning regimens, and the availability of new donor sources, available in the national registries (1,2).
Reduced-intensity conditioning is widely used to avoid the complications of myeloablative conditioning to prepare for HSCT in the adult population. However, this procedure is associated with a high risk of complications that may result in graft loss. To prevent this, it is important to monitor chimerism for early intervention (3). Evaluation of chimerism status at regular intervals is useful to prevent risk of early graft rejection and relapse in patients suffering from malignant diseases. Quantification of the chimerism percentage is also a potential marker of minimal residual disease (MRD) for patients without suitable MRD markers and it provides useful information on graft vs. host disease and graft vs. tumor effects (4–6).
Chimerism analysis is a tool that allows to determination of the genotypic origin of post-transplantation hematopoiesis. Subsequent to HSCT, a patient presenting with 100.0% donor-origin cells during follow up is considered to have the status of complete chimerism (CC), patients in which the donor- and recipient-origin cells coexist have the status of mixed chimerism (MC) (7). Informative genetic markers are used to discriminate between recipient and donor genomes in order to detect the chimerism status (8).
At present, different approaches based on polymerase chain reaction (PCR) amplification of polymorphic DNA sequences (short tandem repeat, STR; single nucleotide polymorphism, SNP; and insertion/deletion, INDEL) are used for chimerism analysis. In the vast majority of laboratories, semi-quantitative fluorescent PCR of STRs is the procedure of choice for diagnostic purposes. The key advantage offered by this method is the highly polymorphic nature of the STR markers, which allow for a high probability of two-genome discrimination. Laboratories currently use different commercial multiplex kits for forensic identifications or in house assays, however have formed consortiums to standardize procedures and set guidelines for the correct interpretation of results, in order to improve this intrinsically semi-quantitative platform with a sensitivity of 1.0–3.0% (9).
Next generation sequencing (NGS) technologies are an innovation in human and animal genomics research, as they are capable of producing 100-fold more data than the most powerful Sanger based capillary sequencers; thus enabling researchers to investigate the large number of queries that remain to addressed (10).
NGS generates hundreds of giga-bases of nucleotide sequences per instrument run and produces this data at a lower cost, thus motivating researchers to use NGS for various purposes: To identify rare variations on the whole genome or on a target sequence, to analyze transcriptome profiling of cells, tissues and organisms and to identify epigenetic markers for disease diagnosis. Progress in the optimization of procedures, in addition to further reduction of costs, are the key factors that will lead to a more extensive uptake of this technique in diagnosis and for practical clinical applications.
NGS provides qualitative and quantitative data. Quantitative data depends on the depth of sequence data collected on each sample and on the quality of the target to expose. For samples with a lower abundance target, many more sequence reads are required to achieve accurate quantification (11). A previous study demonstrated that NGS exhibits sensitivity comparable to that of quantitative PCR (qPCR) in the evaluation of MRD in B cell disorders (12).
In the present study, an Ion AmpliSeq custom chimerism (ACCh) panel and a custom bioinformatic pipeline was created for chimerism quantification by NGS. The first aim was to detect the existence of cells of two origins in chimera samples and then to evaluate the capability of NGS to determine the percentage of the recipient cells.
Materials and methods
DNA sample preparation
The Ethics Committee of the Institute for Maternal and Child Health, IRCCS 'Burlo Garofolo' approved the present study (approval number: Prot. 18/2015, Cl. M/11). Written informed consent was obtained from all the participants.
Total peripheral blood was collected from 10 volunteer donors (V01-V10; 4 males and 6 females) with ages ranging from 20–50 years, and from 2 pediatric patients that underwent allogeneic HSCT (pR1, male, 5 years; and pR2, male, 12 years) and their donors (pD1, male, 9 years; and pD2, female, 25 years; Table I). Written informed consent was obtained from all the participants.
All DNA samples were isolated using the QIAamp DNA Blood kit according to manufacturer's protocol (Qiagen GmbH, Hilden, Germany). The DNA status was evaluated using the NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Inc., Waltham, MA, USA) and the Qubit dsDNA HS Assay kit and the Qubit fluorometer (Thermo Fisher Scientific, Inc.). A DNA stock solution was prepared for all DNA samples at 20.0 ng/μl.
DNA samples from the 10 volunteer donors were randomly paired. A total of 5 artificial chimeric DNA mixtures, as the donor/recipient chimera, were created by diluting DNA with its paired DNA at several percentages of the first DNA for each artificial DNA mixture (aCh).
Ion AmpliSeq custom chimerism panel design
A multi phase strategy was employed to evaluate the main characteristics of the ACCh panel: i) The panel average heterozygosity was assessed around 0.5 for the European population (HapMap Phase 3 CEU population); ii) two SNPs per somatic chromosome, termed 'main SNPs' (mSNPs), were selected and located in two different regions of the same chromosome; iii) the amplicon composition was evaluated according to the following requirements: a) GC percentage ranging between 40.0 and 60.0%; b) presence of one mSNP inside each amplicon; c) mSNP location preferably in the centre of the amplicon; d) absence of INDEL SNPs; e) absence of homopolymers and potential homopolymer generation from SNP variants and their flanking regions; f) absence of flanking SNPs to the mSNPs.
In total, 44 single-nucleotide, biallelic, polymorphisms were selected from the NCBI dbSNPs database (http://www.ncbi.nlm.nih.gov/SNP/, build 138, last database update 28.03.2014; Table II). A total of 4 base sequences including mSNPs were used as target regions for primer design. The primer pool, intended for DNA library construction through multiplex PCR, was defined by Ion AmpliSeq Designer software, version 3.0.1 (Thermo Fisher Scientific, Inc.). A single-tube, 44 primer pair pool was purchased from Life Technologies (Thermo Fisher Scientific, Inc.).
Ion torrent library preparation and sequencing
DNA sample library preparation was performed according to the AmpliSeq Library Preparation protocol (Life Technologies; Thermo Fisher Scientific, Inc.). For each DNA sample, a library was constructed using 10.0 ng genomic DNA through the Ion AmpliSeq Library kit, version 2.0. Subsequently, according to the library preparation protocol, each DNA library was indexed using the Ion Xpress Barcode Adapters kit (Thermo Fisher Scientific, Inc.) and was purified using AMPure XP magnetic beads (Beckman Coulter, Inc., Brea, CA, USA). Each DNA library was then quantified by qPCR using the thermo-cycler 7900HT Fast Real-Time PCR system with the Ion Library TQMN Quantification kit (Thermo Fisher Scientific, Inc.). Template Ion Sphere Particles were arranged using the Ion Personal Genome Machine (PGM) Template OT2 200 kit (Thermo Fisher Scientific, Inc.) and a single end 200 base-read sequencing run was conducted using the Ion Torrent PGM system. Libraries were pooled at 8 pM using the following rates: Donor/recipient, 1:1; and chimera/chimera, 1:1. The recipient/chimera rate was fixed at 1:40 in order to obtain an average coverage of the above libraries around 250X:10,000X. Library pools were sequenced on ion 314 and 316 chip (Table III).
Hotspot panel bed file
A hotspot panel bed file was created using the UCSC Genome Browser (https://genome.ucsc.edu). All SNPs located in the central region of each amplicon were included (NCBI dbSNPs build 138; 'Common SNPs' =286). All the INDELs present across the amplicons, and the SNPs near the 5′ and 3′ ends of the amplicons were excluded from the file. All the above 44 selected SNPs were marked as 'mSNP' and the SNPs belonging to the same amplicon were indexed with the same chromosome/amplicon ID number. Finally, a hotspot panel bed file was created: 'HP286SNPs'.
Genotyping and quantification
Genotyping of all DNA was performed automatically, together with the quantification of all the chimeras, using a custom bioinformatics tool. The code of our tool was written using the Shiny package in R, a web framework to build interactive web applications (https://cran.r-project.org/web/packages/shiny/index.html). A functional diagram of the code is presented in Fig. 1 and the full code is available on request. The code is based on dependencies of the Bioconductor package that must be pre-installed for proper tool functionality.
The custom tool requires as input the sequencing bam files of 'Donor', ' Recipient' and 'Chimeric' patient. Briefly, it uses readGAlignments and pileLettersAt functions, from the GenomicAlignments package (13), to read bam files and extracts the letters/nucleotides into a set of individual genomic positions defined from the bed file. Thresholds for 'Donor' and 'Recipient', homozygous and heterozygous genotyping calls, are settled in base counts frequency ranges of 90.0–100.0% and 30.0–60.0%, respectively. Genotyping calls not included in the thresholds ranges were excluded as unreliable; users can modify the thresholds according to their needs from the user interface. Genotypes from each library were crosschecked to select only SNPs comparable in all conditions. Selected SNPs from donors and recipients were labelled as informative recipient alleles (IRA) according to the following schema: Donor homozygous and recipient heterozygous [Donor (AA) and Recipient (Aa); Donor (aa) and Recipient (Aa)]; donor and recipient homozygous for different alleles [Donor (AA) and Recipient (aa); Donor (aa) and Recipient (AA)].
Only the IRA SNPs tagged as informative were used to calculate the chimera's donor:recipient ratio as median of the allele frequency ratio, while standard error was used to calculate confidence intervals of prediction at 95.0%.
To cross validate the tool, genotyping of all donor and recipient samples was also performed manually obtaining the variant data from the Ion Torrent plugin Variant Caller, version 4.4 using the 'Generic-PGM-Germ Line-Low Stringency' configuration coupled by the HP286SNPs hotspot bed file.
Microsatellite analysis and patient data validation
Multiplex PCR amplification of V01-V10 and aCh1–13 samples, in addition to the patient samples pD2, pR2 and pCh1–6, was performed according to the manufacturer's instructions of the AmpFlSTR Identifiler Plus PCR Amplification kit (Thermo Fisher Scientific, Inc.). Amplicons were resolved on a Genetic Analyzer 3130 and analyzed with GeneMapper software, version 4.1 (Life Technologies; Thermo Fisher Scientific, Inc.).
Patient samples (pCh1–6) were also analyzed by qPCR (data not shown), as previously investigated by Bai et al (14). This analysis was performed as an additional validation method of NGS data, where a discrepancy between NGS and STR data was present.
Results
Ion chips and ACCh panel performance
A total of 7 library pools were loaded and sequenced on ion chip 314. Each pool was comprised of one donor, one recipient and one chimera. In addition, three additional pools were run on ion 316 chip (Table III). The mean values of performance of the 10 runs and of all samples are summarized in Table IV.
NGS genotyping performances using the ACCh panel
A total of 14 DNA samples were genotyped on Ion Torrent PGM using the ACCh panel with the HP286SNPs bed file. Regarding the mSNPs, the Variant Caller output identified that 2 mSNPs (rs121893 and rs12612347) were assigned as 'No Call' in over 50.0% of the genotyping runs due to low quality. The remaining 42 mSNPs were successfully genotyped. Concerning the remaining 242 SNPs, 27 SNPs were assigned as 'No Call', with an average of 10 SNPs per patient. The call of these SNPs failed in two Variant Caller filtering steps: 'Maximum common signal shift' and 'minimum coverage on either strand'.
Genotyping of all samples was additionally performed using our tool with the HP286SNPs bed file (Table V). To identify the IRAs, the data of the donor was compared with the recipient using the tool and manually cross-validated with the Variant Caller genotypes. Inside the genotyping calls of 242 SNPs, a small bias was present between these 2 tools; this is due to the high conserved filters of variant caller, dedicated predominantly for standard sequencing applications, and due to the absence of these filters in our custom tool.
NGS linearity, detection limit and accuracy with the ACCh panel
In order to test the linearity of Ion Torrent PGM with the ACCh panel in a fixed detection range (0.5–100.0%), a series of DNA mixtures was developed, diluting a DNA with its paired DNA at several percentages of the original. In order to increase the genetic marker variability in addition to the biological variability, a total of 12 artificial chimeras (aCh1–12) were prepared from 5 different DNA pairs. Finally a pure DNA (V09) was run as 100.0% DNA (aCh13). Subsequent to Ion Torrent sequencing, using the custom tool, quantitative data for all IRAs of each artificial chimera were obtained.
In addition, to increase the putative points in the dynamic range, the informative alleles of both DNA in the chimeras aCh12 (40.0% of V09 and 60.0% of V10) and aCh11 (20.0% of V07 and 80.0% of V08) were calculated and quantified.
Least-squares analysis of the above putative points identified a clear linearity (R2=0.999; Y=1.008X-0.005) between NGS and the reference values (Fig. 2).
Analyzing the artificial chimeras aCh6–13 by capillary electrophoresis using the STRs markers, chimerism ranging from 4.0–100.0% was detected. Least-squares analysis identified a clear linearity (R2=0.999; Y=1.012X-0.009) between NGS and STRs values (Fig. 3).
In addition, the background of Ion Torrent generated by the ACCh panel was estimated. In this case, the custom tool was used, considering the 'donor' samples as chimera. The average background value at each SNP was estimated at 0.3% (range, 0.0–0.8%) and with a 95.0% confidence interval between 0.1 and 0.5%. Considering the background values and the reported literature on the error rates at each base of NGS technologies (range from 0.04–1.0%) (15), the detection limit of the NGS protocol with the ACCh panel was set at 1.0%, although an artificial chimera was detected at 0.5%.
Finally, considering that the method determined each chimera, calculating the average value of all IRA, the average standard error was used as an indirect marker of accuracy, using the data of all artificial chimeras ranging from 1.0–99.0% (excluding 0.5 and 100.0%). For the dynamic range of 1.0–20.0% of chimeras, the average standard error was calculated at 0.3% with a deviation at 0.2%. For higher values of chimeras, up to 99.0%, the standard error increased up to 1.8% with a maximum deviation at 2.0%.
Patient chimerism evaluation on the NGS platform
Considering the linearity between NGS and reference values and between NGS and STR values using the standard NGS workflow and the custom tool, our workflow was tested in 6 samples of the same patient (pCh1–6) in which the chimerism quantification report, previously performed by microsatellite analysis at different times in an external laboratory (Department of Molecular Medicine, University of Padova, Padova, Italy), presented at least one CC between two MCs (Table I).
NGS analysis detected a mixed chimerism in all samples (pCh1, pCh2, pCh3, pCh4, pCh5 and pCh6) while the microsatellite analysis only in 4 of them (pCh1, pCh3, pCh4 and pCh5). For these 4 MC samples, the percentage of predicted chimerism was equal between the 2 methods of analysis (Fig. 4). In regard to the pCh2 and pCh6 samples, NGS analysis evaluated a mixed chimerism at 3.0 and 2.0%, respectively (Fig. 4).
To confirm the obtained NGS data for the pCh2 and pCh6 samples, the NGS analysis was repeated, and qPCR was additionally perfomed in all patient samples. The results of qPCR analysis were in agreement with that of the NGS and microsatellite data for the pCh1, pCh3, pCh4 and pCh5, and confirmed the NGS data for the pCh2 and pCh6 samples.
Discussion
NGS technologies have revolutionized the field of genomics, and its application has been extended to different fields such as clinical diagnostics and forensic science (15–17). As a result of the continuous development of NGS, several applications previously performed on Sanger sequencing with capillary electrophoresis have been transferred onto the NGS platform, enabling fast and cost-effective generation sequence data with high resolution and accuracy. For this reason, different panels are being developed for the sequencing of genetic mutations involved in human diseases (e.g. MiSeqDx Cystic Fibrosis Clinical Sequencing Assay; Illumina, Inc., San Diego, CA, USA) or in cancer (Ion AmpliSeq BRCA1 and BRCA2 Panel; Life Technologies; Thermo Fisher Scientific, Inc.).
In the field of chimerism quantification by NGS platforms, Debeljak et al (18) reported an innovative and well performed study, using haplotype counting. In addition, Kim et al (19) briefly reported a relative quantification analysis of SNP markers by NGS in one human bone marrow chimerism sample. However, the study by Kim et al (19) was conducted in a 4.9% chimerism sample without any detail concerning the limit of detection, the technical error or additional important technical information and validation data of NGS application in chimerism quantification.
In the present study, a full workflow was designed, and the proposed protocols and a bioinformatics tool were tested for chimerism quantification by NGS. A 44-amplicon custom chimerism panel based on Ion AmpliSeq technology was designed, and in addition a bioinformatics tool dedicated to the genotyping and quantification of NGS data was coded. These resources were created in order to provide a novel tool for the evaluation of the chimerism following allogenic HSCT, thus potentially increasing the number of clinical analyses supported on NGS platforms.
The ACCh proposed panel is composed of 44 amplicons, containing 44 selected mSNPs, of which 2 mSNPs are located in different regions of each somatic chromosome. It is suggested that the different mSNP locations in all somatic chromosome may be useful to avoid predominantly false negatives results caused by chromosomal deletions characteristic of certain malignancies (20). In addition, the bed file uploaded in the custom tool, containing all targeted SNPs, can be modified in order to exclude the SNPs present in chromosome target regions subjected to deletions in a specific patient.
The panel average heterozygosity was assessed around 0.5 for the European population in order to obtain different informative markers for each transplanting pair, for a more precise and robust quantification. The theoretical panel informativity for unrelated donor:recipient values, calculated according to the data present on the NCBI dbSNPs database, was estimated to be approximately 16/42 mSNPs, while for siblings the informativity was estimated at 50.0% (approximately 8/42 mSNPs). In order to increase the informativity of the ACCh panel, an additional 242 selected SNPs present in the targeted regions were included in the bed file. This addition of SNPs experimentally increased the average informativity (Table V).
The ACCh panel reached the limit of detection on the Ion Torrent PGM platform of 0.5%, however, this was updated to the conservative value of 1.0% for two reasons: i) The Ion Torrent error is defined to be between 0.04 and 1.0% (21); and ii) the background of the ACCh panel, based on the IRA data of our experiments, ranges between 0.1 and 0.5%.
Regarding the timing of chimerism analysis, the UK NEQAS Consortium has recommended that results should be assessed in 5 working days from the reception of the sample and in 3 working days for urgent requests (6). The protocol suggested in the current study is feasible in 2 days; the first day for library preparation and quantification and the second for template preparation, run sequence and data analysis. In addition, due to the fact the ACCh protocol suggested in the current study does not present any differences from the standard AmpliSeq Library Preparation and Ion Torrent PGM Run Sequence protocols, it is possible to introduce it to a standard manual library preparation workflow or in a library preparation workstation.
Considering the AmpliSeq library construction protocol (based on multiplex PCR) from Life Technologies (Thermo Fisher Scientific, Inc.), the same concept could be tested on additional NGS platforms, such as Illumina or Roche, according to their library preparation protocols. Concerning the custom bioinformatics tool, any bam and bai file coupled by a bed file, generated from any platform could be used.
At present, the cost of NGS analysis, compared with microsatellite methods, remains high, however considering the continuously reducing cost per NGS run, an NGS-based method for chimerism quantification could be evaluated in the future for its adoption in laboratories with a high volume of activity, and with NGS platforms already in use for other purposes. Although the present study reported a clear correlation between NGS and STR methods and identifed important technical details, further experimental replications are required in order for the NGS protocol to be validated for future laboratory use.