25class FitLadderModel:
 26    def __init__(self, ladder_assigner: PeakLadderAssigner) -> None:
 27        """
 28        Initialize FitLadderModel object with ladder_assigner object and its attributes.
 29
 30        Args:
 31            ladder_assigner (PeakLadderAssigner): Object of PeakLadderAssigner class.
 32
 33        Returns:
 34            None.
 35        """
 36        self.n_knots = 2
 37        self.ladder_assigner = ladder_assigner
 38        self.fsa_file = self.ladder_assigner.fsa_file
 39        self.best_combination = ladder_assigner.assign_ladder_peak_sizes().reshape(
 40            -1, 1
 41        )
 42
 43        self.fit_model()
 44        self.mse, self.r2 = self.model_score()
 45        self.adjusted_baisepair_df = self.generate_adjusted_trace_df()
 46
 47    def fit_model(self) -> None:
 48        """
 49        Fit model based on ladder type.
 50
 51        Returns:
 52            None.
 53        """
 54        ladder = self.fsa_file.ladder
 55        if ladder == "ROX":
 56            self._fit_ROX_ladder()
 57        elif ladder == "LIZ" or ladder == "ORANGE":
 58            self._fit_LIZ_ladder()
 59        else:
 60            raise LadderNotFoundError(f"'{ladder}' is not a valid ladder")
 61
 62    def model_score(self) -> tuple[float, float]:
 63        """
 64        Calculate mean squared error and R-squared score of the fitted model.
 65
 66        Returns:
 67            Tuple (mse, r2).
 68        """
 69        true_Y = self.fsa_file.ref_sizes
 70        predicted_Y = self.model.predict(self.best_combination)
 71        mse = mean_squared_error(true_Y, predicted_Y)
 72        r2 = r2_score(true_Y, predicted_Y)
 73
 74        return mse, r2
 75
 76    def generate_adjusted_trace_df(self) -> pd.DataFrame:
 77        """
 78        Generate a dataframe with adjusted basepairs and peaks.
 79
 80        Returns:
 81            Dataframe with columns time, peaks and basepairs.
 82        """
 83        # Continue loop until all basepairs are unique by changing n_knots in SplineTransformer
 84        for _ in range(10):
 85            df = (
 86                pd.DataFrame({"peaks": self.fsa_file.trace})
 87                .reset_index()
 88                .rename(columns={"index": "time"})
 89                .assign(
 90                    basepairs=lambda x: self.model.predict(
 91                        x.time.to_numpy().reshape(-1, 1)
 92                    )
 93                )
 94                .loc[lambda x: x.basepairs >= 0]
 95            )
 96
 97            if df.shape[0] == df.basepairs.nunique():
 98                logging.info(f"Ladder fitting model: {self.model}")
 99                return df
100            # If not all bp are unique
101            self.n_knots += 1
102            self.fit_model()
103
104        # reaches only if there is a problem
105        raise ModelFittingError(
106            "There is a problem with the fitting of the model to the ladder"
107        )
108
109    def _fit_ROX_ladder(self) -> None:
110        """
111        Fit model with ROX ladder.
112
113        Returns:
114            None.
115        """
116        self.model = make_pipeline(
117            SplineTransformer(
118                degree=4, n_knots=self.n_knots + 4, extrapolation="continue"
119            ),
120            LinearRegression(fit_intercept=True),
121        )
122
123        X = self.best_combination
124        y = self.fsa_file.ref_sizes
125
126        self.model.fit(X, y)
127
128    def _fit_LIZ_ladder(self) -> None:
129        """
130        Fit model with LIZ ladder.
131
132        Returns:
133            None.
134        """
135        # NB!
136        # Changed the SplineTransformer(n_knots=2) instead of 3!!
137        self.model = make_pipeline(
138            SplineTransformer(degree=3, n_knots=self.n_knots, extrapolation="continue"),
139            LinearRegression(fit_intercept=True),
140        )
141
142        X = self.best_combination
143        y = self.fsa_file.ref_sizes
144
145        self.model.fit(X, y)

 13class PeakLadderAssigner:
 14    """
 15    A class that assigns peak sizes to a ladder for use in DNA fragment size analysis.
 16    The peak sizes are generated based on a given size standard and the input parameters.
 17
 18    Attributes:
 19    fsa_file (FsaFile): an object of the FsaFile class that contains parameters for analysis
 20
 21    Methods:
 22    assign_ladder_peak_sizes(): Assigns peak sizes to a ladder based on the given parameters and returns the best combination.
 23    get_peaks(size_standard): Finds peaks in the size standard based on given parameters and returns the peaks.
 24    generate_graph(peaks): Generates a graph based on the peaks and returns the graph.
 25    generate_combinations(graph): Generates combinations of peaks based on the graph and returns the combinations.
 26    get_best_fit(combinations, method): Finds the best fit for the given combinations and returns the best combination.
 27    _polynomial_model_inv_r2_score(ladder, comb): Returns a score based on the polynomial model inverse r2 score.
 28    _max_fractional_deviation_score(ladder, comb): Returns a score based on the maximum fractional deviation.
 29    _max_first_derivative_score(combination): Returns a score based on the maximum first derivative.
 30    _max_second_derivative_score(combination): Returns a score based on the maximum second derivative.
 31    _max_spline_second_derivative_score(combination): Returns a score based on the maximum spline second derivative.
 32    """
 33
 34    def __init__(self, fsa_file: FsaFile) -> None:
 35        """
 36        Initializes the class with an object of the FsaFile class.
 37
 38        Args:
 39        fsa_file (FsaFile): an object of the FsaFile class that contains parameters for analysis
 40        """
 41
 42        self.fsa_file = fsa_file
 43
 44    def assign_ladder_peak_sizes(self):
 45        """
 46        Assigns peak sizes to a ladder based on the given parameters and returns the best combination.
 47
 48        Returns:
 49        np.array: the best combination of peak sizes
 50        """
 51        peaks = self.get_peaks(self.fsa_file.size_standard)
 52        graph = self.generate_graph(peaks)
 53        combinations = self.generate_combinations(graph)
 54        best_combination = self.get_best_fit(combinations)
 55
 56        return best_combination
 57
 58    def get_peaks(self, size_standard) -> np.array:
 59        """
 60        Finds peaks in the size standard based on given parameters and returns the peaks.
 61
 62        Args:
 63        size_standard (np.array): the size standard array
 64
 65        Returns:
 66        np.array: the peaks array
 67        """
 68
 69        peaks_obj = signal.find_peaks(
 70            size_standard,
 71            distance=self.fsa_file.min_interpeak_distance,
 72            height=self.fsa_file.min_height,
 73        )
 74
 75        peaks = peaks_obj[0]
 76        heights = peaks_obj[1]["peak_heights"]
 77
 78        df = pd.DataFrame({"peaks": peaks, "heights": heights})
 79
 80        # TODO
 81        # Dropped the below, then it worked with the new test files
 82        # idxmax = df["heights"].idxmax()
 83        # df = df.drop(idxmax)
 84
 85        peaks_adj = df.nlargest(self.fsa_file.max_peak_count, ["heights"])
 86
 87        return peaks_adj["peaks"].sort_values().to_numpy()
 88
 89    def generate_graph(self, peaks) -> nx.DiGraph:
 90        """
 91        Generates a graph based on the peaks and returns the graph.
 92
 93        Args:
 94        peaks (np.array): the peaks array
 95
 96        Returns:
 97        nx.DiGraph: the graph object
 98        """
 99        G = nx.DiGraph()
100
101        for p in peaks:
102            G.add_node(p)
103
104        i = 0
105        while i < peaks.size:
106            j = i + 1
107            while j < peaks.size:
108                diff = peaks[j] - peaks[i]
109                if diff <= self.fsa_file.max_ladder_trace_distance:
110                    G.add_edge(peaks[i], peaks[j], length=diff)
111                j += 1
112            i += 1
113
114        return G
115
116    def generate_combinations(self, graph):
117        """
118        Generates combinations of nodes from a given directed acyclic graph (DAG)
119        that satisfy certain conditions.
120
121        Parameters:
122        graph (networkx.DiGraph): a DAG representing a state machine
123
124        Returns:
125        A numpy array representing a sequence of nodes that satisfies the
126        conditions of the DAG.
127        """
128
129        # Get start nodes that have zero in-degree
130        start_nodes = [node for node in graph.nodes if graph.in_degree(node) == 0]
131
132        # Get end nodes that have zero out-degree
133        end_nodes = [node for node in graph.nodes if graph.out_degree(node) == 0]
134        if len(start_nodes) == 0 or len(end_nodes) == 0:
135            raise ValueError("Graph does not have start or end nodes")
136
137        # Get all simple paths from start node to end node
138        all_paths = []
139        for start_node in start_nodes:
140            for end_node in end_nodes:
141                paths = nx.all_simple_paths(graph, start_node, end_node)
142                all_paths.extend(paths)
143
144        if len(all_paths) == 0:
145            raise ValueError("No paths found from start to end nodes")
146
147        # Generate combinations of nodes that satisfy certain conditions
148        for p_arr in all_paths:
149            for i in range(0, len(p_arr) - self.fsa_file.ref_count + 1):
150                yield np.array(p_arr[i : i + self.fsa_file.ref_count])
151
152    def get_best_fit(self, combinations, method="2nd_derivative"):
153        """
154        Calculates the best fit for a given set of combinations using a specified method.
155
156        Parameters:
157        combinations (iterable): an iterable containing combinations of nodes
158        method (str): a string indicating the method to use. The default method is
159            "2nd_derivative".
160
161        Returns:
162        A numpy array representing the combination of nodes with the best fit.
163        """
164
165        # Create an empty dataframe to store combinations and scores
166        df = pd.DataFrame()
167
168        # Store the combinations in the dataframe
169        df["combination"] = list(combinations)
170
171        # Calculate the score for each combination using the specified method
172        if method == "2nd_derivative":
173            df["score"] = np.vectorize(self._max_spline_second_derivative_score)(
174                df["combination"]
175            )
176
177        if method == "first_derivative":
178            df["score"] = np.vectorize(self._max_first_derivative_score)(
179                df["combination"]
180            )
181
182        # Sort the dataframe by score in ascending order
183        df = df.sort_values(by="score", ascending=True)
184
185        # Get the best combination (i.e., the one with the lowest score)
186        best = df.head(1)
187
188        # Return the best combination as a numpy array
189        return best.combination.squeeze()
190
191    @staticmethod
192    def _polynomial_model_inv_r2_score(ladder: np.array, comb: np.array) -> float:
193        """
194        Calculates the inverse of the R^2 (coefficient of determination) score for a polynomial regression model.
195
196        Parameters:
197        - ladder (np.array): array containing the ladder positions of a dataset
198        - comb (np.array): array containing the corresponding comb heights of a dataset
199
200        Returns:
201        - float: the inverse of the R^2 score for the polynomial model
202
203        Notes:
204        - This method fits a polynomial regression model of degree 3 to the dataset using NumPy's polyfit function.
205        - It then calculates the predicted values of the comb heights using NumPy's poly1d function.
206        - Finally, it calculates the R^2 score between the predicted comb heights and the actual comb heights using the r2_score function
207          from the scikit-learn library, and returns the inverse of that score.
208
209        """
210        # fit polynomial regression model
211        fit = np.polyfit(ladder, comb, 3)
212
213        # create predicted values using fitted model
214        predict = np.poly1d(fit)
215
216        # calculate R^2 score between predicted and actual values, and return its inverse
217        return 1 - r2_score(comb, predict(ladder))
218
219    @staticmethod
220    def _max_fractional_deviation_score(ladder: np.ndarray, comb: np.ndarray):
221        """
222        Calculates the maximum fractional deviation score for a dataset.
223
224        Parameters:
225        - ladder (np.ndarray): array containing the ladder positions of a dataset
226        - comb (np.ndarray): array containing the corresponding comb heights of a dataset
227
228        Returns:
229        - float: the maximum fractional deviation score
230
231        Notes:
232        - This method calculates the interval lengths between adjacent ladder positions and corresponding comb heights using NumPy's
233          diff function and min-max scaling with a feature range equal to the range of ladder positions.
234        - It then calculates the fractional deviation between the interval lengths and the scaled comb heights, and returns the
235          maximum deviation.
236
237        """
238        # calculate interval lengths and scaled comb height intervals
239        l_intervals = np.diff(ladder)
240        c_intervals = np.diff(minmax_scale(comb, feature_range=(ladder[0], ladder[-1])))
241
242        # calculate fractional deviation between interval lengths and scaled comb height intervals
243        frac_deviation = np.abs(l_intervals - c_intervals) / l_intervals
244
245        # find index of maximum deviation and return the corresponding value
246        max_frac_deviation_idx = np.argmax(frac_deviation)
247        return frac_deviation[max_frac_deviation_idx]
248
249    def _max_first_derivative_score(self, combination: np.ndarray) -> float:
250        """
251        Calculates the maximum first derivative score for the given combination of features.
252
253        Args:
254            combination (np.ndarray): A 1-dimensional numpy array containing a combination of features.
255
256        Returns:
257            float: The maximum first derivative score for the given combination of features.
258        """
259        # Scale the combination to the range of reference sizes
260        comb_scaled = minmax_scale(
261            combination,
262            feature_range=(self.fsa_file.ref_sizes[0], self.fsa_file.ref_sizes[-1]),
263        )
264
265        # Calculate the differences between consecutive values of the scaled combination and the reference sizes
266        diff_intervals = np.diff(comb_scaled) - np.diff(self.fsa_file.ref_sizes)
267
268        # Calculate the absolute gradient of the differences between consecutive values
269        abs_diff_intervals_gradient = np.abs(np.gradient(diff_intervals))
270
271        # Find the index of the maximum absolute gradient
272        max_abs_diff_intervals_gradient_idx = np.argmax(abs_diff_intervals_gradient)
273
274        # Return the maximum absolute gradient
275        return abs_diff_intervals_gradient[max_abs_diff_intervals_gradient_idx]
276
277    def _max_second_derivative_score(self, combination: np.ndarray) -> float:
278        """
279        Calculates the maximum absolute second derivative score of the given combination.
280
281        Args:
282            combination (np.ndarray): An array of feature values for a given combination.
283
284        Returns:
285            float: The maximum absolute second derivative score of the given combination.
286
287        Raises:
288            None.
289
290        """
291        # Scale the combination to match the range of reference sizes.
292        comb_scaled = minmax_scale(
293            combination,
294            feature_range=(self.fsa_file.ref_sizes[0], self.fsa_file.ref_sizes[-1]),
295        )
296
297        # Calculate the difference between the scaled combination and reference sizes.
298        diff_intervals = np.diff(comb_scaled) - np.diff(self.fsa_file.ref_sizes)
299
300        # Calculate the absolute second derivative of the difference intervals.
301        abs_second_derivative = np.abs(np.gradient(np.gradient(diff_intervals)))
302
303        # Find the index of the maximum absolute second derivative score.
304        max_second_derivative_idx = np.argmax(abs_second_derivative)
305
306        # Return the maximum absolute second derivative score.
307        return abs_second_derivative[max_second_derivative_idx]
308
309    def _max_spline_second_derivative_score(self, combination: np.ndarray) -> float:
310        """
311        Calculates the maximum absolute value of the second derivative of a given combination of values using a UnivariateSpline.
312
313        Args:
314        - self (object): an instance of the class containing the function
315        - combination (np.ndarray): an array of values representing a combination
316
317        Returns:
318        - max_value (float): the maximum absolute value of the second derivative of the given combination
319
320        """
321        # create a UnivariateSpline object using the reference sizes and the given combination
322        spl = UnivariateSpline(self.fsa_file.ref_sizes, combination, s=0)
323
324        # calculate the second derivative of the spline
325        der2 = spl.derivative(n=2)
326
327        # evaluate the second derivative at each reference size and return the maximum absolute value
328        max_value = max(abs(der2(self.fsa_file.ref_sizes)))
329
330        return max_value

16def baseline_arPLS(y, ratio=0.99, lam=100, niter=1000, full_output=False):
17    L = len(y)
18
19    diag = np.ones(L - 2)
20    D = sparse.spdiags([diag, -2 * diag, diag], [0, -1, -2], L, L - 2)
21
22    H = lam * D.dot(D.T)  # The transposes are flipped w.r.t the Algorithm on pg. 252
23
24    w = np.ones(L)
25    W = sparse.spdiags(w, 0, L, L)
26
27    crit = 1
28    count = 0
29
30    while crit > ratio:
31        z = linalg.spsolve(W + H, W * y)
32        d = y - z
33        dn = d[d < 0]
34
35        m = np.mean(dn)
36        s = np.std(dn)
37
38        w_new = 1 / (1 + np.exp(2 * (d - (2 * s - m)) / s))
39
40        crit = norm(w_new - w) / norm(w)
41
42        w = w_new
43        W.setdiag(w)  # Do not create a new matrix, just update diagonal values
44
45        count += 1
46
47        if count > niter:
48            print("Maximum number of iterations exceeded")
49            break
50
51    if full_output:
52        info = {"num_iter": count, "stop_criterion": crit}
53        return z, d, info
54    else:
55        return z

 14class FsaFile:
 15    def __init__(
 16        self,
 17        file: str,
 18        ladder: str,
 19        normalize: bool = False,
 20        trace_channel: str = "DATA1",
 21        size_standard_channel: str = None,
 22        min_interpeak_distance: int = None,
 23        min_height: int = None,
 24        max_ladder_trace_distance: int = None,
 25    ) -> None:
 26        """
 27        Constructs an FsaFile object with the given parameters.
 28
 29        Args:
 30            file (str): The path to the sequencing file.
 31            ladder (str): The name of the ladder used for sequencing.
 32            normalize (bool): Whether to normalize the data using the arPLS algorithm.
 33            trace_channel (str): The channel to extract the trace data from.
 34            size_standard_channel (str): The channel to extract the size standard data from.
 35            min_interpeak_distance (int): The minimum distance between peaks in the ladder.
 36            min_height (int): The minimum height of peaks in the ladder.
 37            max_ladder_trace_distance (int): The maximum distance between the ladder and the trace.
 38
 39        Returns:
 40            None
 41        """
 42        peak_count_padding = 3
 43        self.file = Path(file)
 44        self.file_name = self.file.parts[-1]
 45
 46        # Extract data from the sequencing file
 47        if ladder not in LADDERS.keys():
 48            raise LadderNotFoundError(f"'{ladder}' is not a valid ladder")
 49        self.ladder = ladder
 50        self.fsa = SeqIO.read(file, "abi").annotations["abif_raw"]
 51        self.trace_channel = trace_channel
 52        self.normalize = normalize
 53
 54        # Extract data from the ladder reference
 55        self.ref_sizes = LADDERS[ladder]["sizes"]
 56        self.ref_count = self.ref_sizes.size
 57
 58        # Use default values if nothing is given
 59        self.size_standard_channel = size_standard_channel or LADDERS[ladder]["channel"]
 60        self.min_interpeak_distance = (
 61            min_interpeak_distance or LADDERS[ladder]["distance"]
 62        )
 63        self.min_height = min_height or LADDERS[ladder]["height"]
 64        self.max_ladder_trace_distance = (
 65            max_ladder_trace_distance or LADDERS[ladder]["max_ladder_trace_distance"]
 66        )
 67        self.max_peak_count = self.ref_count + peak_count_padding
 68
 69        # Normalize data if requested
 70        if normalize:
 71            self.size_standard = np.array(
 72                baseline_arPLS(self.fsa[self.size_standard_channel])
 73            )
 74            self.trace = np.array(baseline_arPLS(self.fsa[trace_channel]))
 75        else:
 76            self.size_standard = np.array(self.fsa[self.size_standard_channel])
 77            self.trace = np.array(self.fsa[trace_channel])
 78
 79    def __str__(self):
 80        """
 81        Returns a string representation of the FsaFile object.
 82
 83        Returns:
 84            str: A string representation of the object.
 85        """
 86        return f"""
 87            FsaFile-object with following parameters:
 88            
 89            File: {self.file}
 90            Filename: {self.file_name}
 91            Size standard channel: {self.size_standard_channel}
 92            Ladder name: {self.ladder}
 93            Number of ladder steps: {self.ref_count}
 94            Minimum interpeak distance: {self.min_interpeak_distance}
 95            Minimum height: {self.min_height}
 96            Minimum Ladder trace distance: {self.max_ladder_trace_distance}
 97            Maximum peak count: {self.max_peak_count}
 98            Normalized data: {self.normalize}
 99            Trace Channel: {self.trace_channel}
100            Ladder Sizes: {self.ref_sizes}
101            """

 25class PeakAreaDeMultiplex:
 26    def __init__(
 27        self,
 28        peaks: PeakFinder,
 29        cutoff: float = None,
 30    ) -> None:
 31        self.peaks = peaks
 32        self.cutoff = cutoff or None
 33        self.peaks_dataframe = self.peaks.peaks_dataframe
 34        self.peak_information = self.peaks.peak_information
 35        self.file_name = self.peaks.file_name
 36        self.peaks_index = self.peaks.peaks_index
 37        self.found_peaks = self.peaks.found_peaks
 38        self.area_plots = []
 39
 40        self.find_peak_widths()
 41        # divade peaks based on their assay they belonging
 42        self.divided_assays = self.divide_peak_assays()
 43        # how many assays does this sample contain?
 44        self.number_of_assays = len(self.divided_assays)
 45        # divide all peaks in each assay into separate dataframes
 46        self.divided_peaks = [self.divide_peaks(x) for x in self.divided_assays]
 47        # logging
 48        logging.info(
 49            f"""
 50        Number of assays found: {self.number_of_assays}
 51        Number of peaks found: {self.peak_information.shape[0]}
 52        """
 53        )
 54
 55    def __iter__(self):
 56        return PeakAreaDeMultiplexIterator(self.number_of_assays)
 57
 58    def find_peak_widths(self, rel_height: float = 0.95):
 59        widths = peak_widths(
 60            self.peaks_dataframe.peaks,
 61            self.peaks_index,
 62            rel_height=rel_height,
 63        )
 64
 65        df = pd.DataFrame(widths).T
 66        df.columns = ["x", "peak_height", "peak_start", "peak_end"]
 67
 68        self.peak_widths = (
 69            df.assign(peak_start=lambda x: np.floor(x.peak_start).astype(int))
 70            .assign(peak_end=lambda x: np.ceil(x.peak_end).astype(int))
 71            .assign(peak_name=lambda x: range(1, x.shape[0] + 1))
 72            .merge(self.peak_information, on="peak_name")
 73        )
 74
 75    def divide_peak_assays(self) -> list[pd.DataFrame]:
 76        """
 77        Divide the peaks belonging to different assays based on their assay number
 78        """
 79        df = self.peak_widths
 80        return [df.loc[df.assay == x] for x in df.assay.unique()]
 81
 82    def divide_peaks(self, assay: pd.DataFrame, padding: int = 4) -> list[pd.DataFrame]:
 83        # add some padding to the left and right to be sure to include everything in the peak
 84        return [
 85            self.peaks_dataframe.iloc[x.peak_start - padding : x.peak_end + padding]
 86            for x in assay.itertuples()
 87        ]
 88
 89    def fit_lmfit_model(self, peak_finding_model: str, assay_number: int):
 90        if assay_number >= self.number_of_assays:
 91            raise IndexError(
 92                f"""
 93                The sample only contains {self.number_of_assays} assays. 
 94                Use a number inside of the range.
 95                Indexing starts at 0.
 96                """
 97            )
 98
 99        if peak_finding_model == "gauss":
100            model = GaussianModel()
101        elif peak_finding_model == "voigt":
102            model = VoigtModel()
103        elif peak_finding_model == "lorentzian":
104            model = LorentzianModel()
105        else:
106            raise NotImplementedError(
107                f"""
108                {peak_finding_model} is not implemented! 
109                Options: [gauss, voigt, lorentzian]
110                """
111            )
112
113        fit_df = []
114        fit_params = []
115        fit_report = []
116        for df in self.divided_peaks[assay_number]:
117            df = df.copy()
118            y = df.peaks.to_numpy()
119            # CHANGED to time instead of basepair
120            x = df.time.to_numpy()
121
122            params = model.guess(y, x)
123            out = model.fit(y, params, x=x)
124
125            fit_df.append(df.assign(fitted=out.best_fit, model=peak_finding_model))
126            fit_params.append(out.values)
127            fit_report.append(out.fit_report())
128
129        # Update the instances of the model fit
130        self.fit_df = fit_df
131        self.fit_params = fit_params
132        self.fit_report = fit_report
133
134    # TODO
135    # Fix so that the cutoff is a range or something else.
136    def calculate_quotient(
137        self,
138    ) -> None:
139
140        """
141        :Params:
142        """
143        areas = np.array([x["amplitude"] for x in self.fit_params])
144
145        right_by_left = (
146            self.cutoff is None
147            or pd.concat(self.fit_df).basepairs.mean() >= self.cutoff
148        )
149        # if there only is 1 peak, return 0
150        if len(areas) == 1:
151            self.quotient = 0
152            return
153
154        # if there only are 2 peaks, return the quotient
155        if len(areas) == 2:
156            # left peak divided by right peak
157            if not right_by_left:
158                self.quotient = areas[0] / areas[1]
159                return
160            # right peak divided by left peak
161            self.quotient = areas[1] / areas[0]
162            return
163
164        # return the last peak divided by the mean of the peaks to the left of it
165        self.quotient = areas[-1] / areas[:-1].mean()
166
167    def peak_position_area_dataframe(
168        self, assay_number: int, name: str
169    ) -> pd.DataFrame:
170        """
171        Returns a DataFrame of each peak and its properties
172        """
173        dataframes = []
174        for i, _ in enumerate(self.fit_df):
175            report = self.fit_report[i]
176            r_value = float(re.findall(r"R-squared *= (0\.\d{3})", report)[0])
177            df = (
178                self.fit_df[i]
179                .loc[lambda x: x.peaks == x.peaks.max()]
180                .assign(area=self.fit_params[i]["amplitude"])
181                .assign(r_value=r_value)
182                .assign(peak_name=f"Peak {i + 1}")
183                .drop(columns="time")
184                .reset_index(drop=True)
185                .rename(
186                    columns={
187                        "peaks": "peak_height",
188                        "fitted": "fitted_peak_height",
189                    }
190                )
191                .drop_duplicates("peak_name")
192                .assign(file_name=self.file_name)
193                .assign(assay_name=name)
194            )
195            dataframes.append(df)
196
197        self.assay_peak_area_df = pd.concat(dataframes).assign(
198            quotient=self.quotient,
199            peak_number=lambda x: x.shape[0],
200            assay_number=assay_number + 1,
201        )
202
203    def plot_areas(self, peak_finding_model: str):
204        fig_areas, axs = plt.subplots(
205            1, len(self.fit_df), sharey=True, figsize=(20, 10)
206        )
207        # if there is only one peak
208        if len(self.fit_df) == 1:
209            axs.plot(self.fit_df[0].basepairs, self.fit_df[0].peaks, "o")
210            axs.plot(self.fit_df[0].basepairs, self.fit_df[0].fitted)
211            axs.set_title(f"Peak 1 area: {self.fit_params[0]['amplitude']: .1f}")
212            axs.grid()
213        # if more than one peak
214        else:
215            for i, ax in enumerate(axs):
216                ax.plot(
217                    self.fit_df[i].basepairs,
218                    self.fit_df[i].peaks,
219                    "o",
220                )
221                ax.plot(self.fit_df[i].basepairs, self.fit_df[i].fitted)
222                ax.set_title(
223                    f"Peak {i + 1} area: {self.fit_params[i]['amplitude']: .1f}"
224                )
225                ax.grid()
226
227        fig_areas.suptitle(f"Quotient: {self.quotient: .2f}")
228        fig_areas.legend(["Raw data", "Model"])
229        fig_areas.supxlabel("basepairs")
230        fig_areas.supylabel("intensity")
231        plt.close()
232
233        return fig_areas
234
235    def fit_assay_peaks(
236        self,
237        peak_finding_model: str,
238        assay_number: int,
239        name: str,
240    ) -> None:
241        """
242        Runs fit_lmfit_model, calculate_quotient and peak_position_area_dataframe
243        """
244        self.fit_lmfit_model(peak_finding_model, assay_number)
245        self.calculate_quotient()
246        self.peak_position_area_dataframe(assay_number, name)
247        # area plots
248        area_plot = self.plot_areas(peak_finding_model)
249        self.area_plots.append(area_plot)
250        return self.assay_peak_area_df
251
252    def assays_dataframe(self, peak_finding_model: str = "gauss"):
253        dfs = []
254        for i in self:
255            if self.peaks.custom_peaks is not None:
256                name = self.peaks.custom_peaks.name.unique()[i]
257            else:
258                name = ""
259            dfs.append(self.fit_assay_peaks(peak_finding_model, i, name))
260        df = pd.concat(dfs, ignore_index=True)
261        # initialize attribute self.final_df
262        self.final_df = df
263        return df