Optimizing sequence design strategies for perturbation MPRAs: a computational evaluation framework

Abstract

The advent of perturbation-based massively parallel reporter assays (MPRAs) technique has facilitated the delineation of the roles of non-coding regulatory elements in orchestrating gene expression. However, computational efforts remain scant to evaluate and establish guidelines for sequence design strategies for perturbation MPRAs. In this study, we propose a framework for evaluating and comparing various perturbation strategies for MPRA experiments. Within this framework, we benchmark three different perturbation approaches from the perspectives of alteration in motif-based profiles, consistency of MPRA outputs, and robustness of models that predict the activities of putative regulatory motifs. While our analyses show very similar results across multiple benchmarking metrics, the predictive modeling for the approach involving random nucleotide shuffling shows significant robustness compared with the other two approaches. Thus, we recommend designing sequences by randomly shuffling the nucleotides of the perturbed site in perturbation-MPRA, followed by a coherence check to prevent the introduction of other variations of the target motifs. In summary, our evaluation framework and the benchmarking findings create a resource of computational pipelines and highlight the potential of perturbation-MPRA in predicting non-coding regulatory activities.

Graphical Abstract

Figure 1.

An outline of the framework for evaluation of perturbation-based massively parallel assays technique. In the ‘Three sequence designing approaches’ box, we used the ‘GATA_known9’ motif as an example. In detail, the GATA motifs are a group of sequences conforming to the consensus WGATAR (W = A or T and R = A or G) (marked by the wavy underline), that can be recognized and bound by GATAbinding transcription factors (45).

Figure 2.

Evaluations of perturbation-wise metrics. (A) Examples of ‘hit’ and ‘fail’ sequences. Please refer to the Supplementary Notes for the full perturbation sequences. (B) A comparison of hit rates among three perturbation approaches. The ‘N/A’ category represents the sequences that are excluded from this study because their barcodes failed the sequencing quality check (Supplementary Notes). (C) Examples of ‘perturbed’ and ‘non-perturbed’ sequences. (D) A comparison of perturbation rates among three perturbation approaches.

Figure 3.

Evaluations of motif-based metrics. (A) An example of calculating perturbation specificity. Refer to the Supplementary Notes for the full perturbation sequences. (B) A comparison of perturbation specificity among three perturbation approaches. Significant P values (P < 0.05) are shown in red. (C) An example of calculating ‘newly introduced target motifs per sequence’. Please refer to the Supplementary Notes for the full perturbation sequences. (D) A comparison of ‘newly introduced target motifs per sequence’ among three perturbation approaches.

Figure 4.

Evaluation of general alteration in the number of motifs. (A) Toy examples of calculating general alteration in the number of motifs. (B–D) The results for motif perturbations: (B) the number of gained motifs, (C) the number of lost motifs and (D) the net change in the number of motifs. Significant P values (P < 0.05) are shown in red. (E–G) The results for random perturbation sequences: (E) number of gained motifs, (F) number of lost motifs and (G) net change in the number of motifs. Significant P values (P < 0.05) are shown in red.

Figure 5.

Assessment of the important features representing perturbation sequences. (A) The number of important features shared by three perturbation approaches. (B) Top 30 important features of each perturbation approach. The names of features that are shared by at least two perturbation approaches are marked in bold. (C) Gene ontology enrichment analysis of the top 2500 genes represented by the TF binding factors.

Figure 6.

Consistency of MPRA outputs among three perturbations. (A) Number of sequences that share the same FRS identities. The bars are colored by activators (red) and repressors (blue). In the ‘intersection type’ matrix. The percentages are row-normalized, indicating the proportion of sequences belonging to different intersection types within each perturbation approach. (B) The correlation of Log₂FC between motif_PERT1 and motif_PERT2. Each dot is a perturbation sequence and is colored by the time point. (C) The correlation of Log2FC between motif_PERT1 and motif_PERT3. (D) The correlation of Log2FC between motif_PERT2 and motif_PERT3.

Figure 7.

Comparison of the Log₂FC among three perturbation approaches. The Log₂FC values are separated by time point before being compared among three perturbation approaches.

Figure 8.

Performance of classification models. (A) The area under the receiver-operating characteristic curve (AUROC) of different classification models. Asterisks/ns indicate levels of statistical significance, calculated by pairwise Wilcoxon rank sum tests (P-value < 0.05*, < 0.01**, < 0.001***, < 0.0001****; ns, non-significant). (B) A summary of the mean ± standard deviation values for AUROCs of classification models.

Figure 9.

Performance of regression models. (A) The Pearson correlation coefficients of different regression models. Asterisks/ns indicate levels of statistical significance, calculated by pairwise Wilcoxon rank sum tests (P-value < 0.05*, < 0.01**, < 0.001***, < 0.0001****; ns, non-significant). (B) A summary of the mean ± standard deviation values for Pearson correlation coefficients of regression models.