An image-to-answer algorithm for fully automated digital PCR image processing

Lab on a Chip PAPER Cite this: Lab Chip,2022,22,1333 Received 25th December 2021, Accepted 23rd February 2022 DOI: 10.1039/d1lc01175h rsc.li/loc An image-to-answer algorithm for fully automated digital PCR image processing† Zhiqiang Yan, ‡ a Haoqing Zhang, ‡ a Xinlu Wang, a Martina Gaňová, b TomášLednický, b Hanliang Zhu, a Xiaocheng Liu, a Marie Korabečná, c Honglong Chang a and Pavel Neužil * a The digital polymerase chain reaction (dPCR) is an irreplaceable variant of PCR techniques due to its capacity for absolute quantification and detection of rare deoxyribonucleic acid (DNA) sequences in clinical samples. Image processing methods, including micro-chamber positioning and fluorescence analysis, determine the reliability of the dPCR results. However, typical methods demand high requirements for the chip structure, chip filling, and light intensity uniformity. This research developed an image-to-answer algorithm with single fluorescence image capture and known image-related error removal. We applied the Hough transform to identify partitions in the images of dPCR chips, the 2D Fourier transform to rotate the image, and the 3D projection transformation to locate and correct the positions of all partitions. We then calculated each partition's average fluorescence amplitudes and generated a 3D fluorescence intensity distribution map of the image. We subsequently corrected the fluorescence non-uniformity between partitions based on the map and achieved statistical results of partition fluorescence intensities. We validated the proposed algorithms using different contents of the target DNA. The proposed algorithm is independent of the dPCR chip structure damage and light intensity non-uniformity. It also provides a reliable alternative to analyze the results of chip-based dPCR systems. Introduction The polymerase chain reaction (PCR) method was initially developed for enzymatic deoxyribonucleic acid (DNA) amplification using thermal cycling. It allowed further qualitative characterization of the amplified DNA fragments, 1 especially for gene mutation analysis and the detection of polymorphic DNA sequences. The PCR technique evolved with the introduction of intercalating dyes or target sequence hybridization with fluorescently labeled probes, leading to real-time monitoring of reaction kinetics and subsequent quantitative analysis of nucleic acids in samples, known as real-time or quantitative PCR (qPCR). 2,3 It was then further developed and, with the advent of micromachining and miniaturization, the hand-held qPCR for point-of-care applications was introduced. 4 The development of the microfluidics technique enabled new digital PCR (dPCR) technology, a variant of end-point detection PCR techniques. A compartmentalized sample forms a binary system of thousands or millions of subsamples, each containing a single copy or no copy of the target DNA. 5 This technique allows absolute quantification of selected target DNA sequences and, most importantly, performs multiplex PCR with minimal interaction between different targets, 6 distinguishing rare DNA targets from abundant ones. 7 The compartmentalization of the original sample is achieved either by filling predefined partitions forming a chip-based dPCR (cdPCR) 8 or by forming individual droplets as microreactors of a droplet-based dPCR (ddPCR). 9,10 The PCR protocol is then conducted with several thermal cycles (typically between 40 and 45), and the fluorescence amplitude (F) from each well or droplet (partition) is determined once the protocol is completed for either single DNA target 11 or several DNA targets thus performing target multiplexing. 12 The partitions with amplified targets are defined as positive partitions (PW), exhibiting an increased value of F when compared with the original value. The partitions Lab Chip,2022,22, 1333–1343 | 1333This journal is © The Royal Society of Chemistry 2022 a School of Mechanical Engineering, Department of Microsystem Engineering, Northwestern Polytechnical University, 127 West Youyi Road, Xi'an, Shaanxi, 710072, P.R. China. E-mail: pav[email protected] b Central European Institute of Technology, Brno University of Technology, Purkyňova 656/123, 612 00, Brno, Czech Republic c Institute of Biology and Medical Genetics, First Faculty of Medicine, Charles University and General University Hospital in Prague, Albertov 4, 128 00 Prague, Czech Republic †Electronic supplementary information (ESI) available. See DOI: 10.1039/ d1lc01175h ‡These two authors contributed equally, and both are considered as first authors. Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online View Journal | View Issue 1334 |Lab Chip,2022,22, 1333–1343 This journal is © The Royal Society of Chemistry 2022 without targets, defined as negative partitions (NW), show significantly lower Fvalues. Image processing methods, including a series of image corrections and Fvalue analysis, determine the reliability of the dPCR results. 13 In this study, we mainly discuss cdPCR-based image processing, typically determining partition positions, error correction, and subsequent Fextraction. Current techniques to determine partition positions are based on pattern recognition, chip structure, or alignment marks. 14 A passive reference dye is added to the dPCR master mix to provide the background fluorescence signal, help the software recognize PW and NW in the chip, and improve statistical precision. 15 The pattern recognition method identifies the partitions with fluorescence using a machinelearning algorithm, 16 a region-growing algorithm, 17 and a clustering algorithm. 18–20 However, chip defects, such as dust and scratches, often result in bright image areas incorrectly identified as PW, while the partitions not filling the sample may be identified as NW, affecting the result's accuracy. The chip structure-based partition positioning method excludes false recognition as only the fluorescence in each partition's area is captured. 15,21,22 This method utilizes the chip's design to determine the position of all partitions, removing fixed pattern noise in the image caused by dust or defects. Recently, we proposed a method using the chip's alignment marks to determine the partition locations by capturing them using a bright-field imaging method, registering their centers, and then capturing the partition area by fluorescence imaging in six blocks, followed by image-stitching and subsequent image processing. 23 The chip structure-based partition positioning method improves recognition precision. However, feature structure damage, such as missing alignment marks, affects the accuracy of partition position determination. In addition, this method does not consider chip illumination inhomogeneity. 21 Once the partitions are recognized, the Fvalues are extracted, and a fixed F(F T ) threshold is defined to distinguish PW from NW. 24–26 However, the chip illumination non-uniformity and the emitted Fvalue affect the reliability of the extracted results. 27 The adaptive F T method reduces the influence of light non-uniformity. 28 First, the image is split into several segments, assuming that each segment's light distribution is uniform. Each segment's threshold is then defined, and the PW number is counted in each segment. Alternatively, a machine-learning algorithm determines whether partitions with fluorescence are PW or NW. 29 Additionally, the details of image processing algorithms in commercial dPCR systems are not available to users; thus, it is impossible to elaborate on their merits or demerits. Only a few papers have described chip-based image processing methods based on deep learning. 30 However, it requires a high computer configuration level and an extensive training dataset, making detection complex and time-consuming. Previous research has proposed a random background transfer-based image processing algorithm to simplify the training dataset significantly, as it only requires three experimental images. 31 However, it is used in augmenting images to lower the effect of uneven illumination rather than to remove it. Here, we develop an image-to-answer processing algorithm that requires only a single fluorescence image capture using a dPCR chip. The process eliminates capturing a set of bright field images, registering the chip location, capturing a set of fluorescence images, stitching them together, and then processing them as previously required. 23 Moreover, the proposed method simplifies the image processing procedure as it does not rely on a training dataset. Additionally, we implement several image-related error corrections, making this proposed dPCR image processing method versatile. Furthermore, we locate the partitions' positions using a new Fig. 1 A dPCR hardware module: (A) schematic design of the dPCR system, comprising a LED, fluorescence filter set, TEC, temperature sensor, Z-stage, X–Ystage, spacer, lens, and SLR camera; (B) photograph of an assembled dPCR system showing details of the TEC with the dPCR chip (inset). Lab on a ChipPaper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online Lab Chip,2022,22, 1333–1343 | 1335This journal is © The Royal Society of Chemistry 2022 method combined with pattern recognition and chip structure, providing an option to help identify the location of the partitions in the chip. Together with the previously published method of dPCR image emulation, 32 both works complement each other for dPCR image emulation and processing. Materials and methods dPCR chips We utilized a 9 ×9mm 2 dPCR chip with 26 448 partitions, each with a target diameter of 50 μm split into six blocks, designed using the script-based Nanolithography Toolbox, 33,34 a family of dPCR chips with the same previously described size and basic structure. 23 The dPCR chip was made of Si using a single lithography step, followed by deep reactive ion etching with a target depth of 30 μm. The chip was first immersed in H 2 SO 4 /H 2 O 2 (piranha) with a 96%/30% concentration and a 4 : 1 ratio solution at a temperature of ≈120 °C for ≈10 min. This piranha solution treatment demineralizes organic materials, removing all traces of previous dPCR experiments and making the dPCR chip reusable. This treatment also forms OH − groups at the Si surface, making it hydrophilic. The prepared PCR master mix subsequently filled the partitions in the dPCR chip, and the chip was sealed with a microscope cover glass of ≈10 × 10 mm 2 and ≈170 μm in size and thickness, respectively. The cover glass was spincoated with ≈50 μm polydimethylsiloxane and then coated with ≈2μm thick parylene C by physical vapor deposition. The dPCR hardware We designed a fluorescence imaging system using a fluorescein isothiocyanate (FITC) cube from an optical Table 1 Sequences of the DNA template and both primers DNA template CGCTGAATCCTGCGGACGACCCTTCTCGGGGTCGCTTGGGACTCTCTCGTCCCCTTCTCCGTCTGCCGTTCCGACCGACCACGGGG CGCACCTCTCTTTACGCGGACTCCCCGTCTGTGCCTTCTCATCTGCCGGACCGTGTGCACTTCGCTTCACCTCTGC Forward primer 5′-GTCGCTTGGGACTCTCTC-3′ Reverse primer 5′-GCAGATGAGAAGGCACAGA-3′ Fig. 2 Schematic diagram of an automated dPCR image processing method: step 1: system calibration: radial lens distortion k, image skew angle θ, and light non-uniformity correction mask Aare determined for dPCR image processing system calibration. Step 2: fluorescence analysis: after the amplification, Fvalues are extracted from the dPCR chip. Then the dPCR image is analyzed to calculate the cn of the target. Lab on a Chip Paper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online 1336 |Lab Chip, 2022, 22, 1333–1343 This journal is © The Royal Society of Chemistry 2022 microscope (Fig. 1A). The cube was mounted on an X–Ystage with a light-emitting diode (LED) with a principal wavelength and a maximum power of 470 nm and 30 W, respectively; the optical power was controlled by an LED controller using the pulse-width modulation (PWM) method. A commercial complementary metal-oxide semiconductor (CMOS) singlelens reflex (SLR) camera, with a 35 ×24 mm 2 full-frame and 20.2 Mpixel imager size, can capture images with 16or 8-bit gray-scale resolution. In this study, we operated the images with 16-bit resolution, which is 256 times more precise than 8-bit resolution. The camera, with a 180 mm focal length macro lens and 3.5 minimum aperture, set to 8.0, one extension tube with a total length of 150 mm, and a 1.4× telephoto image enhancer, was attached to the cube via a tube spacer. The dPCR chip was placed under the cube's center on a brass plate. The plate was soldered onto a thermoelectric cooler (TEC) powered by electrical current pulses from an H-bridge system (Fig. 1B). The TEC was fixed to the top of a copper cooler mounted on a Z-stage and located on two sliders, making chip filling and replacement more convenient. The dPCR chip was placed in a brass holder on top of the TEC to fix the chip's location (Fig. 1B inset), protecting the sealed chip from accidental movement. The holder also allowed fast heat transfer from the TEC to the solution in the partitions during thermal cycling. The brass plate's thermal cycling was monitored with a Pt100 temperature sensor. The temperature was controlled using the PWM method, allowing proportional, integrative, and derivative methods to control the dPCR chip's temperature from a personal computer. We used the same setup as before for verification of the dPCR image generation method. 32 DNA template and PCR protocol The chip was filled with a ≈0.15 μM concentration fluorescein solution, which had ≈0.17 ± 0.01 (mean ± standard deviation from three experiments) of the Fvalue of a double-stranded DNA amplicon in the presence of EvaGreen after 40 cycles of PCR, determined using a commercial qPCR system (ESI†section A). First, the chip Fig. 3 Partition location algorithm: (A) original image and coordinate system; (B) schematic illustration of the image θcalculation algorithm; (C) rotated image; (D) emulated dPCR fluorescence image used for algorithm verification; (E) partitions detected (red) by CHT, the enlarged local map shows the performance in a corner; (F) four corner points (red) are calculated by four boundary lines (green) as control points to locate all partitions based on a predefined chip structure; (G) location of all partitions after employing the correction algorithm. The enlarged local map shows the performance in a corner. Table 2 Value of θcalculated by the proposed algorithm Predefined θ(°)−14 −13 −12 −11 −10 −9−8−7−6−5 Detected θ(°)−14.00 −13.00 −11.99 −11.01 −10.00 −9.00 −8.00 −7.00 −6.01 −5.01 Predefined θ(°)−4−3−2−10 12345 Detected θ(°)−4.01 −3.00 −2.00 −0.99 0 0.99 2.00 3.00 3.99 5.00 Predefined θ(°)67891011121314— Detected θ(°) 6.00 6.99 7.99 8.99 9.99 11.01 12.00 13.00 14.00 — Lab on a ChipPaper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online Lab Chip,2022,22, 1333–1343 | 1337This journal is © The Royal Society of Chemistry 2022 photograph was used as an input calibration dPCR chip image. Then, we chose a synthesized hepatitis B virus gene as a DNA template. 35 The amplicon length was 102 base pairs (bp), and the template and primer sequences are given in Table 1. The PCR master mix consisted of ≈0.96 μLTaq polymerase, ≈5μLbuffer,≈1μL of forward and reverse primers with a final concentration of ≈100 nM, ≈0.4 μL fluorescein in a concentration of ≈0.15 μM helping to increase the contrast between PW and NW, ≈1.2 μL EvaGreen, ≈1μL DNA template, ≈0.94 μLbovine serum albumin at a concentration of 20 mg mL −1 and ≈13.5 μL deionized H 2 O. The polymerase concentration was increased due to the unfavorable surface-to-volume ratio. 36 We performed a three-step PCR protocol, a hot start to activate the polymerase at a set temperature of 95 °C for 60 s followed by 40 PCR cycles, consisting of the following steps: denaturation at a set temperature of 95 °C for 15 s, annealing at 56 °C for 20 s, and elongation at 72 °C for 30 s. Principle of the method We used a single fluorescence image to automatically determine the copy number (cn) of targets in a dPCR chip. The image processing algorithm included system calibration and fluorescence analysis (Fig. 2). The algorithm workflow is shown in ESI†section B. We determined a series of image correction algorithm parameters, including radial lens distortion correction coefficient k, image skew angle θ, and light nonuniformity correction mask A. These parameters were used to correct the dPCR chip's fluorescence image after DNA amplification. Point of note: these correction operations are optional; some can be skipped because optical systems Fig. 4 Global non-uniformity correction: (A) schematic illustration of two light non-uniformity correction methods, generating 3D (left) and 2D (right) fluorescence intensity distribution maps after non-uniformity correction; (B) a comparison of the fluorescence images before and after correction. We increased the contrast and used different colors (from blue to red) to denote different gray-scale values; (C) a comparison of the fluorescence intensity distribution curves before and after correction. A histogram (left) of Fis the distribution value after correction. The enlarged local map (right) shows residual non-uniformity that was not eliminated. Table 3 Non-uniformity index Mand Gauss distribution curve width σbefore and after correction Original image Corrected image using method 1 Corrected image using method 2 M(%) 3.44 0.95 0.40 σ(pixel) 7.96 3.20 1.55 Correction rate of M(%) —72.38 88.37 Correction rate of σ(%) —59.80 80.53 Lab on a Chip Paper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online 1338 |Lab Chip, 2022, 22, 1333–1343 This journal is © The Royal Society of Chemistry 2022 with acceptable global non-uniformity do not require correction. System calibration The image processing algorithm started with several preprocesses from the acquired calibration dPCR chip fluorescence image. First, the original image's Canon raw 2 (CR2) format was converted into tag image file format (TIFF), maintaining 16-bit resolution. Then, the background area of the image was cropped so that the only area of the dPCR chip with partitions remained, and the red, green, and blue (RGB) image was converted into the gray-scale format. Image conversion Radial distortion can be caused by the lens, 37,38 including barrel and pincushion types, and is a common form of image error that negatively affects the determination of partition positions. If such a lens is used, correction of this distortion is essential. Therefore, we mainly used a commercial SLR camera with automated lens recognition and a barrel distortion correction program. However, if the image is captured using a different camera, such as a smartphone camera, the kimplementation algorithm should be employed (Fig. S2†). Details of this algorithm are shown in ESI† section C. Detection of partitions and image skew correction The camera is always rotated (skewed) by an unpredictable value of θwhen taking a fluorescence image of the dPCR chip (Fig. 3A), as it is impossible to place the dPCR chip under the camera with a value of θ=0°. If the rotation angle is not compensated, the subsequent location correction is affected, leading to a false recognition. The image skew correction allows the subsequent algorithm to detect the image boundary easily. Also, the chip layout we processed was based on an array with a hexagonal pattern rather than a square. This pattern makes chip boundary recognition more complicated if the skew is still present. We proposed a fast Fourier transform (FFT)-based skew correction algorithm (Fig. 3B) to determine the image's θvalue. This algorithm started with partition detection using a circle Hough transform (CHT) (Fig. 3B, I 1 ), 39,40 determining the circle Fig. 5 Example of fluorescence analysis of an amplified dPCR chip, showing the workflow of the algorithm: (A) image of original dPCR chip filled with fluorescein; (B) image of a light non-uniformity correction mask. In this image, the pixel color was linked to its correction coefficient; image of a dPCR chip (C) before and (D) after skew and non-uniformity correction; (E) four boundary lines (green), four corner points (black), and all partitions (white) in the image of the chip; (F) the locations of all partitions were corrected based on the chip design shown in red; (G) the occurrence histogram as a function of the Fvalue is built based on the extracted Fvalues. Table 4 Processed results of chips with a set partition diameter of 50 μm. Each experiment was conducted once Calculated λ C (copies per partition) 0.1 0.2 0.3 0.4 0.7 0.8 Calculated percentage of PW number (%) 9.52 18.13 25.92 32.97 50.34 55.07 Extracted λ E (copies per partition) 0.105 0.173 0.237 0.347 0.564 0.676 Extracted percentage of PW number (%) 9.87 15.92 20.16 26.24 39.47 46.01 Error in λ(%) 0.49 2.67 6.31 5.33 13.58 12.41 Lab on a ChipPaper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online Lab Chip, 2022, 22, 1333–1343 | 1339This journal is © The Royal Society of Chemistry 2022 indexed iwith a radius r i and center location (x i ,y i ) defined by the pixel coordinate (Fig. 3A). The partitions' radius and layout are two essential parameters for the dPCR chip image processing we describe in our work. Because we used the dPCR chip we designed and fabricated, we had access to these vital design parameters. Should we process somebody else's design, we first need to extract these parameters, which is not difficult using modern optical microscopes with size measurement. Afterward, I 1 was converted into a binary image I 2 (Fig. 3B). The gray level with a radius of three pixels in the circle areas centered at (x i ,y i ) was set to 1, and the gray level in other areas was set to 0. This method removed the nonuniformity of subsequently used spectrum image I 3 . Thus, we did not need to use a multi-level segmentation algorithm. We then applied 2D FFT and generated the spectrum image I 3 ; the gray value of each pixel showed periodicity at a unique frequency and direction (Fig. 3B). We decreased the contrast of I 3 and binarized the generated spectrum image with an empirical threshold value of 0.5. Subsequently, a line Hough transform (HTM) was employed to detect the first eight distinct lines on I 3 . The median rotation angle of these lines was regarded as θin the fluorescence image. Finally, we rotated the original fluorescence image to get the skew-corrected image I rot (Fig. 3C). The partitions' previously detected center position was also updated using a transform formula: xnew ¼xbase þxorig −xbase  ·cos −θðÞþyorig −ybase  ·sin −θðÞ ynew ¼ybase −xorig −xbase  ·sin −θðÞþyorig −ybase  ·cos −θðÞ 8 > < > : ; (1) where (x new ,y new )inI rot corresponds to (x orig ,y orig ) in original image I 1 , and (x base ,y base ) denotes the coordinates of the center of I 1 . After that, we verified the algorithm to determine the θ value of the dPCR chip image, using a series of emulated images (Fig. 3D) with predefined θvalues from −14°to 14°. 32 These emulated images contained different normalized F values, such as partitions with 80% defined as PW with DNA, and a lower level of 30%, representing NW and unfilled (empty) partitions. A comparison between a predefined angle in an emulated image and the calculated angle extracted from it is shown in Table 2. The mean value of the absolute residual angle was only 0.0057°; this was considered an acceptable value for the subsequent image processing algorithm. Location correction Some partitions, especially the ones near the chip corners, were not detected using the CHT method (Fig. 3E) due to fabrication defects on the chip and incomplete filling of the chip with the solution. Therefore, we proposed a 3D projection transformation method based on an automatic partition location correction algorithm to locate all partitions on the chip. We analyzed the image processed by the CHT method, locating the upper, lower, left, and right boundary partitions. First, we searched for characteristics of partitions neighboring the analyzed partition. Once we found that the analyzed partition had no neighbors on the left-hand side, the partition was regarded as a potential left boundary partition. Thus, the four boundary lines were fitted to the Table 5 Processed data from commercial dPCR chips Calculated cn 2 20 Calculated percentage of PW number (%) 0.04 0.41 Extracted cn 214 Extracted percentage of PW number (%) 0.04 0.29 Error in cn (%) 0 30% Fig. 6 Fvalue analysis of a single block of an amplified dPCR chip with a partition diameter of 20 μm as an example: (A) a fluorescence image of a single block after thermal cycling; (B) the histogram of occurrence as a function of Fvalue. Lab on a Chip Paper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online 1340 |Lab Chip, 2022, 22, 1333–1343 This journal is © The Royal Society of Chemistry 2022 corresponding potential boundary partitions using the random sample consensus (RANSAC) method. 41 This method fitted lines on the boundaries, eliminating wrongly identified boundary partitions (ESI†section D). Next, the coordinates of the four corner points such as top left, top right, bottom left, and bottom right were calculated as cross points of the four boundary lines (Fig. 3F). These four corner points were used as control points to calculate the 3D projection transformation matrix from the chip image to a predefined chip structure. Finally, the positions of all partitions were calculated by performing inverse 3D projection transformations (Fig. 3G) (ESI†section E). Global non-uniformity corrections Non-uniform illumination on dPCR chips negatively affects the precise detection of Fvalues. In addition, illumination non-uniformity could result in overlapping fluorescence distribution of PW and NW, making it difficult to distinguish the two types. Moreover, the fluorescence system readout might not report uniform detection across the imaging device; these global non-uniformities should be corrected. A dPCR illumination system was proposed to reduce the influence of light non-uniformity based on the improved illumination setup. 27 However, the demonstrated device did not completely eliminate the related problem. Therefore, we proposed a global non-uniformity correction algorithm, making a mask from a set of images of a dPCR chip filled with fluorescein solution, a method similar to the fixed pattern noise reduction utilized in infrared focal plane arrays. 42 This algorithm automatically generates a correction mask and suppresses global non-uniform illumination. We used two methods: the first was based on the approximation of non-uniformity by a 3D parabolic function extracted from a set of images, and the second was by the smoothing method. The first steps were almost identical in both methods; generating the correction mask was the only difference. The correction algorithm's first step was to extract Fvalues in each partition. Here, we averaged the gray value in the circle area defined by the center locations (x i ,y i ) and mean radius r=r i _to obtain value F i . This gray value was used to quantify the Fvalues of partition i. Afterward, the relative F values F iN of partition iwere calculated for each partition. As for the PW, F iN was calculated by dividing F i by the mean value of positive partitions F 2 , while for the negative partitions, the F iN value was set to 1. The F iN was affected by structure defects, incomplete filling, and image noise caused by optical system. We employed a median fitting method to smooth F iN in the partitions i. The filtered gray value was set to the median gray value in the partitions, including its neighbor partitions and itself. Finally, a 3D map (Fig. 4A, S 1 ) showed the illumination intensity distribution on the chip, constructed based on F iN and the partitions' location. Method 1: 3D parabolic function fitting method This method was conducted by performing a surface function fitting on the 3D map of S 1 using the function: f(x,y)=P 00 +P 10 x+P 01 y+P 20 x 2 +P 11 xy +P 02 y 2 , (2) where (x,y) is the coordinate of a pixel and f(x,y) is the light non-uniformity coefficient for location (x,y). Then, a new 3D map showing the 3D parabolic function was constructed (Fig. 4A, S 2 ). After obtaining the polynomial surface function, a correction mask matrix with the same size fluorescence image was generated by calculating f(x,y) in each pixel (x,y) (Fig. 4A, M 1 ). This matrix is used to correct the subsequent actual dPCR fluorescence image by dividing its gray value in each pixel. It should also be noted that the non-uniformity of fluorescence intensity is caused by many factors, including non-uniform illumination, non-uniform film coating, and unpredictable image random noise. Therefore, this 2D parabola fitting method is designed only for systematic errors, such as non-uniform illumination and fluorescence capture. This method cannot correct other defects, such as random glass coating non-uniformity. Method 2: the image smoothing method We performed surface adjacent pixel smoothing and created a 3D map, exhibiting non-uniformity caused by uneven illumination, fluorescence imaging, and the parylene coating of the glass coverslip. First, we performed a Gauss filtering algorithm to smooth the light non-uniformity mask S 1 initially extracted from the fluorescein-filled chip using a set standard deviation of 20 pixels (Fig. 4A, S 3 ). Then, a correction mask M 2 , similar to M 1 , was generated by extracting a relative gray value on S 3 in each pixel (x,y). The smoothing method was used for mask smoothing while preserving its intact edges. We calculated the averaged non-uniformity index Mto quantify the effect of our non-uniformity correction algorithm, as described below: M¼ FiN−1jj ·100%:(3) A comparison of the fluorescence intensity images before and after correction is shown in Fig. 4B (above). Then, we extracted the gray value for each partition and built a distribution curve (Fig. 4C). Furthermore, we compared Mto the Gauss distribution width σvalue extracted by performing Gauss distribution function fitting before and after correction (Table 3). The adjacent smoothing method is significantly better than the 3D parabolic function fitting, but it is specific to the glass used for the experiment. The first method based on 3D surface fitting is less precise but more universal. We combined feature-based and chip-structure-based methods to locate all partitions on the chip. The chip was Lab on a ChipPaper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online Lab Chip,2022,22, 1333–1343 | 1341This journal is © The Royal Society of Chemistry 2022 filled with fluorescein solution with a concentration of ≈0.15 μM having ≈17% of Fvalues of a PCR master mix after completion of PCR. We first captured the chip's fluorescence image using a 20.2 Mpixel commercial camera and close lookup lens. Next, CHT method was performed to detect partitions with Fvalues in the detectable range in the captured image. Then, we determined the image barrel distortion and corrected it. The image's skew angle was found using the FFT method. Using the extracted partitions' locations, we determined the coordinates of the four corners (control points). Finally, we performed a 3D projection transformation using these control points and the known design to locate all partitions, including those with Fvalues out of range. According to our experiment, all partitions on the chip were located with accuracy better than ≈5μm. Subsequently, we performed either 3D surface fitting to create a mask or formed the mask by direct 3D smoothing and generated a non-uniformity correction mask for both methods. Applying this mask lowered the non-uniformity by ≈72.4% and 88.4% using 3D parabolic surface fitting and adjacent smoothing, respectively. This technique was applied before in conjunction with DNA melting curve analysis for dPCR chip temperature non-uniformity determination. 43 Here we conducted the following measurements to recognize and identify both NW and PW: adding the fluorescein with a concentration of ≈0.15 μM in the dPCR master mix, capturing the dPCR image using 16-bit resolution, and performing global non-uniformity corrections to reduce the non-uniform Fvalue distribution of both PW and NW on the chip, thus inhibiting the false identification caused by non-uniform illumination. Results and discussion Copy number calculation methods We created a histogram of an average gray-scale value extracted from each partition, using the image with 16-bit gray-scale resolution with a BIN number of 256. A histogram from extracted Fvalues may show several clusters by adjusting the image contrast, even with a single target. 44 An ideal histogram of dPCR results based on a single target test consists of two major clusters, one having a lower Fvalue representing the NW and the other with a higher Fvalue representing the PW. Thus, the Gauss distribution function can be used to fit these clusters in a histogram: fFðÞ¼A1e−2F−F1 ðÞ 2 w12þA2e−2F−F2 ðÞ 2 w22;(4) where f(F) is the number of occurrences of partitions with F values, A 1 and A 2 denote the areas of each peak, F 1 and F 2 are the mean Fvalues for NW and PW, and w 1 and w 2 are the half-widths of each peak of Gauss distribution. Then, the F distribution functions of NW and PW can also be determined: fNW FðÞ¼A1e−2F−F1 ðÞ 2 w12(5) and fPW F ðÞ ¼A2e−2F−F2 ðÞ 2 w22;(6) where f NW (F) and f PW (F) are the Gauss distribution functions for NW and PW, respectively. Then, the numbers for NW and PW can be calculated using the function: NW FðÞ¼ ðfNW·dF Nb (7) and PW F ðÞ ¼ðfPW·dF Nb ;(8) where N b is the BIN number of the histogram, here N b = 256. Unfortunately, there is a chance that multiple targets can land in the same partitions. 45 In this case, the data extraction described above using two Gauss distribution functions would not yield correct results. The probability P(k) of having ktargets in a partition is given by Poisson distribution: 23 PkðÞ¼ λk k!e−λ;(9) where λis the average cn in each partition. Here, we set k values in a range from one to five, as the P(k) value for k>5 can be neglected, and P(k= 0) determines the probability of partitions having no target. The total number of partitions (N) is defined as the sum of NW and PW: N¼NW þPW ¼Ne−λþλ 1!e−λþλ2 2!e−λþλ3 3!e−λþλ4 4!e−λþλ5 5!e−λ  ¼NX k¼5 k¼0 PkðÞ (10) We can extract NW from the histogram by integrating f(F)in the selected range and dividing it by the BIN number, calculating the value of λ: NW ¼N·e−λ;thus λ¼−ln NW N(11) Once we know λ, we can calculate the total value of cn: cn ¼NX k¼n k¼1 k·PkðÞ (12) Experimental results First, we filled a dPCR chip with fluorescein as a calibration chip (Fig. 5A) and created a fluorescence non-uniformity correction mask for the subsequent calibration (Fig. 5B). Lab on a Chip Paper Open Access Article. Published on 24 February 2022. Downloaded on 3/9/2023 6:27:49 AM. This article is licensed under a Creative Commons Attribution 3.0 Unported Licence. View Article Online