A Bayesian approach to accurate and robust signature detection on LINCS L1000 data

Abstract

Motivation: LINCS L1000 dataset contains numerous cellular expression data induced by large sets of perturbagens. Although it provides invaluable resources for drug discovery as well as understanding of disease mechanisms, the existing peak deconvolution algorithms cannot recover the accurate expression level of genes in many cases, inducing severe noise in the dataset and limiting its applications in biomedical studies.

Results: Here, we present a novel Bayesian-based peak deconvolution algorithm that gives unbiased likelihood estimations for peak locations and characterize the peaks with probability based z-scores. Based on the above algorithm, we build a pipeline to process raw data from L1000 assay into signatures that represent the features of perturbagen. The performance of the proposed pipeline is evaluated using similarity between the signatures of bio-replicates and the drugs with shared targets, and the results show that signatures derived from our pipeline gives a substantially more reliable and informative representation for perturbagens than existing methods. Thus, the new pipeline may significantly boost the performance of L1000 data in the downstream applications such as drug repurposing, disease modeling and gene function prediction.

Availability and implementation: The code and the precomputed data for LINCS L1000 Phase II (GSE 70138) are available at https://github.com/njpipeorgan/L1000-bayesian.

Supplementary information: Supplementary data are available at Bioinformatics online.

Fig. 1.

Illustration of pipeline for robust L1000 perturbagen signature detection

Fig. 2.

(a) The calibrated reference expression values of the invariant set compared with the original ones. The error bars show the ranges corresponding to ±1 SD. The data points are horizontally offset for better visibility. (b) The distribution of the slope and ${\chi }^2$ of the fitted relationship between unscaled expression values and calibrated reference values

Fig. 3.

The median and 68% confidence interval of the best fit DOF of Student’s t-distribution. Each data point corresponds to the DOFs of all peaks within its neighborhood (±0.01 in terms of log2 peak location). Our choice of three DOF is shown by a gray line

Fig. 4.

The relationship between the scale parameter of the Student’s t-distribution and the log2 expression level, i.e. the center of an isolated peak

Fig. 5.

Hexplots between the true peak and recovered peak positions from Bayesian MLE, k-means and AGMM. The darkness of the hexagon indicates the number of points in it. The Pearson correlation coefficient and MSE for each method are shown in each figure

Fig. 6.

Typical examples where discrepancies happens between Bayesian MLE, L1000 level 2 data and AGMM. For panel (a–d), the results from L1000 peak deconvolution and AGMM are shown in red and green arrows, where the thick and thin arrows indicate the peaks with high (2/3) and low (1/3) abundances, respectively. The results from our Bayesian method are shown as probability distributions in thick and thin blue curves, respectively. The peak locations used in Bayesian MLE are the fluorescent intensities where the probability distributions reach their maxima. The examples are from well REP.A028_MCF7_24H_X2_B25_D11

Fig. 7.

(a) Receiver operating characteristic (ROC) curves for replicates identification. Expression profiles from our Bayesian pipeline, Bayesian MLE and L1000 standard pipeline are tested with GSEA under different query sizes. The three methods are labeled as Bayesian, Bayesian MLE and L1000 level 4 in the figure. (b) The comparison of the median quantile (FPR at $\text{TPR}=0.5$)

Fig. 8.

(a) ROC curves for similar perturbagen recognition based on combined z-scores from Bayesian, Bayesian MLE and L1000 level 5 data. The area under the curve is also shown for each method. (b) The comparison of the median quantile (FPR at TPR = 0.5)