PDOs can be established from pre-treated advanced colorectal cancer patient samples

A total of 34 patients with advanced CRC were included in the present study between November 2018 and November 2020, and 50 samples were collected, mainly from liver metastasis. The primary tumor was biopsied only when metastases were not accessible. In the presence of more than one liver metastasis, a sample was taken from each one. We included both naïve and heavily pre-treated patients. Organoids’ growth was observed in 29 lines from 22 patients, with an overall success rate of 59%. All sample characteristics are shown in Supplementary Table S1. Most biopsies were performed at surgery. PDOs establishment was successful in 90% of cases without chemotherapy during the previous 6 months, dropping to 51% if preoperative chemotherapy had been administered, although the difference was not statistically significant. The PDOs’ establishment rate was not affected by primary tumor location (PTL, left vs. right vs. rectum), tissue biopsy site (primary vs. metastasis), or RAS/RAF mutational status. In 19 samples from 12 patients no organoids could be grown, this culture failure owing mainly to lack of initial growth, arrested growth or initial contamination. We observed heterogeneous growth behavior among our PDOs. Some could be established as long-term cultures (more than 3 months in culture and over 10 passages) while some others could only grow as short-term cultures (1–3 months in culture and 1–9 passages) [38] (Supplementary Fig. S1). Long-term cultures displayed an exceptional capacity for multiple freeze thaw cycles. No specific molecular features were identified as predicting long-term versus short-term culture establishment.

PDOs recapitulate morphology and immunohistochemistry and could be derived from tissues with low cellularity

PDOs were stained with hematoxylin and eosin (H&E) to show cellular architecture, which faithfully reproduced the morphology seen in culture (Fig. 1A, top and middle panel). Cellular architecture presents similarities with corresponding tissues (Fig. 1A, bottom panel).

Fig. 1

PDOs recapitulate morphologic features of original tissue and can be obtained from low cellularity biopsies. (A) Culture images of some PDOs lines (top panel; scale bar 50 μm, 20X), compared with corresponding H&E staining (hematoxylin and eosin) (middle panel; scale bar 50 μm, 20X) and with matched tissue stained with H&E (bottom panel; scale bar 100 μm, 10X). (B) Corresponding percentage of tumor cells assessed as percentage of neoplastic cells with respect to the total amount of viable cells in the tissue sample from which the culture was derived, assessed by the pathologist. (C) IHC for CDX2 and CK20 antibodies in a selection of PDOs. CNO75 is a normal mucosa organoid from a CRC patient, used as control for both markers. (Scale bar 50 μm, 20X). (D) Lack of expression of MLH1 and PMS2 in a core biopsy from a liver metastasis of patient 50 and its corresponding derived PDOs. RTO2 is a pMMR model used as positive control. (Scale bar 50 μm, 20X). T: tissue

To elucidate whether tumor cell percentage could impact on the capacity to obtain organoids, each tissue section destined for organoid culture was revised by a dedicated pathologist. A cut-off of 30% tumor cells was applied to classify our samples. Our data indicate that PDOs can be grown independently of tumor cell proportion found in the original tissue (two-tailed Fisher’s exact test p = 0.4568, Fig. 1B).

IHC staining shows that PDOs are positive for CDX2 and CK20, intestinal markers that are both employed in standard diagnosis of CRC (Fig. 1C). The only exception was seen in mCTO66S3, which tested CK20 negative. The original tissue displayed weak and parceled expression, with most cells clearly negative (Supplementary Fig. S2). Indeed, loss of CK20 expression may be considered a negative prognostic marker in some settings [39]. CDX2 and CK20 expression of corresponding tissues is showed in Supplementary Fig. S2, where a weak expression of CK20 and no CDX2 staining in mCTT47 was observed, probably due to low cellularity, high level of fibrosis and intratumoral heterogeneity. On the other hand, our PDO mCTO50 and mCTO50B, derived from a patient with MSI high showed lack of expression of MLH1 and PMS2 proteins, as was seen in the original tissue (Fig. 1D, Supplementary Fig. S2).

PDOs recapitulate the genomic and transcriptomic profile of the original tumor

To genomically compare organoids and their original tissues, we employed NGS analysis using a customized panel [19]. A high correlation was found regarding germline variants, single nucleotide variants (SNVs), short insertions and deletions (r > 0.8, p < 0.0001, Supplementary Fig. S3A-B). The percentage of tumor cells, necrosis, mucinous proportion, and tumor stroma ratio were calculated by a dedicated pathologist. These indices did not interfere with the ability of PDOs to reproduce the original tissue genomic profile (Supplementary Fig. S3C), except when low cellularity was present (Supplementary Fig. S3B).

As expected, the most prevalent mutations affected TP53, KRAS and APC (Fig. 2A), and besides these, PIK3CA, PTEN and SMAD4 mutations were also detected, encompassing all relevant mutations in CRC. The main driver oncogenic mutations and hotspot variants are represented in Table 1. In most cases, allelic fraction (AF) was higher in PDOs (paired t test, p < 0.0001), reflecting their enrichment in epithelial cells. Moreover, CytoScan HD SNP-array was conducted on PDOs and matched fresh tissue to detect copy number alterations and LOH. Likewise, PDOs recapitulate the overall copy number variation profile, and are significantly enriched (Fig. 2B). CytoScan results were validated by ddPCR (Supplementary Fig. S3D).

Fig. 2

PDOs recapitulate genomic and transcriptomic profile of original tissues. (A) Concordance of SNVs and InDels between PDOs and corresponding tissues. On the left percentage of genomic alterations detected across the study subjects. (B) Whole Genome View representation of long-term PDOs and corresponding tissues according to Cytoscan HD. Data is expressed as the weighted log2 ratio of the copy number on the left Y-axis, and the chromosome number on the X-axis. 0 corresponds to diploid, upper and lower spikes indicate gain and loss regions, respectively. mCTO50 tissue was not available. (C) Heatmap showing the Pearson correlation coefficient (color key) between tumors (rows) and organoids (columns) based on the normalized counts as described in materials and methods section. Dendrograms show the hierarchical clustering based on the complete method and Euclidean pair-wise distance. Different passages from the same organoid culture are included

Table 1 Allelic fraction (AF) distribution in PDOs and matched tissues of main driver, hotspot and not hotspot mutations with potential relevance. Each mutation has been manually reviewed in Integrative Genome Viewer (IGV). Only mutations with an AF of at least 5% in PDOs have been included

RNA-seq was also performed to compare gene expression between PDOs and corresponding tissue and to analyze whether it was stable over time in culture. This is a key issue, considering PDOs as in vitro “avatars” of patient tumors. Three biological replicates deriving from distinct culture passages were employed for each organoid line, except for RTO2 and RTO7 for which two replicates were employed, due to technical issues. Overall, PDOs and tissues cluster separately because of the co-existence of more cell types in tissues and stromal genes (Supplementary Fig. S3E). However, narrowing down the analysis to non-stromal genes contribution (supplementary methods), significant expression correlation can be observed between each PDOs line and its corresponding tissue (Fig. 2C). The two exceptions were mCTO50B and mCTO47, which were derived from a tissue with low cellularity (less than 10%). Interestingly, PDOs from different passages exhibit a similar gene expression profile, indicating that the influence of culture condition is minimal and that this profile is stable over time.

PDOs show differential anti-tumor response to standard and experimental treatments

To study whether PDOs could be an appropriate model to assess differential response to standard therapies, we exposed our long-term models to several approved drugs. The reproducibility of drug treatment results in PDOs was confirmed across several experiments (Supplementary Fig. S4A). Table 2 summarizes the clinical and molecular characteristics of patients whose models were used.

Table 2 Clinical and molecular characteristics of patients from which long-term cultures have been generated

All the established PDOs exhibited a heterogeneous response to chemotherapy (Fig. 3A). As an example, mCTO50B (a PDO obtained after chemotherapy) had a worse response than mCTO50, which was generated from the initial tumor before any therapy was given, so the former could represent a potential model of chemo-resistant disease.

Fig. 3

PDOs response to drug screening. (A) Log transformed dose–response curves in selected standard drugs and non-standard (C) drugs. (B) Z-score Ln-AUCs heatmap (red: no response; blue: good response) for standard treatment. (D) Z-score Ln-AUCs heatmap (red: no response; blue: good response) for non-standard treatment (lower panel), matched with genomic data. In the upper panel are depicted main gene mutations, in the following loss of heterozygosity (LOH) and in the last copy number variations

SN-38 (the active metabolite of irinotecan) was highly active in vitro, and only RTO2 displayed lower sensitivity compared to the other lines in our panel. Adding oxaliplatin to 5-FU resulted in a general decrease in PDOs viability (Supplementary Fig. S4B). Along these lines, combination with SN-38 seemed to be more effective than single agent exposure (Supplementary Fig. S4B, not statistically significant). Anti-EGFR treatment with erlotinib showed modest activity in all lines, including the RAS wild type RTO7 PDO. Regorafenib also showed modest in vitro activity, whereas trifluridine/tipiracil (TAS-102) was more active in all PDOs except mCTO66S3 (Supplementary Fig. S4C). A matrix representation of normalized Z-score Ln-AUCs (Fig. 3B) identifies relative resistance/sensitivity to a single drug among different PDOs.

To further explore sensitivity to other anti-cancer agents, our PDOs were exposed to different drugs not currently included in the standard of care for advanced CRC. Agents targeting pathways altered in our models were selected. Response was matched with genomic alterations in all models (Fig. 3C-D). Differential responses were seen for palbociclib, trametinib, alpelisib, the BET inhibitor birabresib (Fig. 3C). For example, RTO7 showed a dramatic response to palbociclib, a selective inhibitor of CDK4/6 kinases, in comparison with the other PDOs, particularly RTO2 and mCTO66S3. Considering genomic data, the lower responder mCTO66S3 has the higher copy number level of both CDK4 and CDK6, reported as a mechanism of resistance [40]. Additionally, RTO7 has the highest copy number of c-MYC, potentially indicating cell cycle activation [41,42,43]. Moreover, this line has a SMAD4 p.P356S missense variant, in a hotspot region with LOH, thus explaining the nearly 100% mutant allelic fraction (Table 1). Despite the lack of functional validation, it is predicted to be oncogenic by several databases, leading to a loss-of-function of the protein [44] and is associated with increased c-MYC activity [45]. Indeed, the high response to palbociclib observed in RTO7 certainly warrants further research, as all these data point to a “MYC-addicted” phenotype of RTO7.

The BET inhibitor birabresib (Fig. 3C, middle) showed good activity in all lines, particularly in SMAD4 loss RTO7 and mCTO66S3 (Fig. 3D, Supplementary Fig. S4D, S5A), confirming synthetic lethality via restoration of MYC inhibition [45]. MEK-inhibitor trametinib showed modest activity in all lines except for KRAS mutant RTO2 (Fig. 3C), while PI3Kα-inhibitor alpelisib showed modest activity in all PDOs except for a slightly more in PIK3CA mutant line RTO2 (Fig. 3C, upper right). Intriguingly, the PI3K-mTORC1/2 inhibitor GSK458 was significantly active, not only in PIK3CA mutant RTO2, but also in mCTO66S3 (Fig. 3C, lower right). In contrast, mCTO50B had a significantly worse response. Indeed, this PDO has a frameshift mutation of the PTEN gene (p.K267Rfs*9) reported as likely oncogenic in OncoKB™, and we showed that it leads to reduced gene expression and a loss of PTEN expression per IHC (Supplementary Fig. S5B). We also performed additional drug testing without observing relevant activity (Supplementary Fig. S4D).

PDOs differentially replicate patient response depending on the site of disease

To determine whether our PDOs would also reproduce patient response to treatment, we retrieved clinical information and matched drug response observed for mCTO50 and RTO7 (Fig. 4), in which follow up data events were available. mCTO50 was generated from liver metastasis of a patient with MSI-high CRC before initiating any anti-tumor therapy. This patient received cytoreductive chemotherapy with FOLFOXIRI (a combination of 5-FU, leucovorin, oxaliplatin and irinotecan), presenting appropriate tumor shrinkage and allowing liver metastases resection. When tumor response was tested in vitro, mCTO50 was significantly more sensitive to 5-FU and oxaliplatin than mCTO50B, which was generated from the remaining residual disease tissue obtained at surgery. The drug treatment effect in this PDO was comparable to that observed in the patient as mCTO50. The patient had disease progression shortly after surgery, indicating that chemotherapy failed to exert a long-term effect. A similar degree of sensitivity to irinotecan was observed in both PDOs, however, prompting us to hypothesize that resistance mechanisms to this drug were not yet established.

Fig. 4

Patient’s response matched with PDOs drug assays, RTO7 in the upper panel and mCTO50 in the lower panel. Each CT scan response evaluation is compared with corresponding treatment administered in the corresponding PDOs. ≠ means discordant response, = means concordant response. Mo(s): months in terms of progression-free survival; St: stage; FOLFIRI: 5-fluorouracil and irinotecan; FOLFOX: 5-fluorouracil and oxaliplatin; FOLFOXIRI: 5-fluorouracil, oxaliplatin and irinotecan; BSC: best supportive care; SD: stable disease; PR: partial response; PD: progressive disease; R1: microscopic residual tumor

RTO7 exhibited a good response to irinotecan-based chemotherapy, received as part of the patient’s first line treatment regimen (Table 2). Conversely, patient response was poor (showing progression at first CT scan). The only concordant response observed was for the anti-EGFR agent and for TAS-102 (i.e. lack of response). However, it should be taken into consideration that the organoid line derived from the primary location rather than liver metastases, which could not be biopsied. This highlights the intrinsic biological heterogeneity of CRC and the need to biopsy metastatic sites in the context of translational studies.

Proteotranscriptomic analysis in PDOs in relation to oxaliplatin or palbociclib sensitivity

Our aim was to uncover functional correlates between baseline proteotranscriptomic expression and drug sensitivity in PDOs. To assess the steady-state protein abundance, a SWATH-MS based label-free quantitative proteomics analysis was performed as a more reproducible strategy in a translational context. The relative level of the 1157 proteins quantified in all samples (FDR < 1%) helped us to distinguish each PDO by principal component analysis (PCA) (Supplementary Fig. S5C) and these were consistently reproduced at different culture passages. We focused on understanding the molecular mechanisms involved in sensitivity to palbociclib and the lack of response to oxaliplatin, where the most striking differences were observed among PDOs. Three biological replicates proceeding from different culture passages were used for each PDOs line.

RTO2, RTO7 and mCTO50 showed better response to oxaliplatin than mCTO50B and mCTO66S3 (unpaired t-test p = 0.0038) (Fig. 3A). A proteomics analysis (Supplementary file S1) identified 95 differentially expressed proteins (51 up-regulated, 44 down-regulated) in non-responder compared to responder organoids. Among the proteins identified with the greatest positive and negative fold change (Fig. 5A), we found several previously associated with oxaliplatin resistance. As an example, resistant PDOs showed higher levels of PRDX6, which is a negative regulator of ferroptotic cell death [46], a process that enhances CRC sensitivity to oxaliplatin [47]. Another up-regulated protein was ALDH9A1, a member of the ALDH family of proteins which are involved in aldehyde detoxification which in turn is associated with acquired chemoresistance in colorectal cancer cells [48]. These same models exhibited lower levels of NDRG1, a replication stress response protein which can inhibit epithelial-mesenchymal transition (EMT), a process that has been associated with phenotypes of chemoresistance [49], and of CDH17, which has been reported as a marker of good response to 5-fluorouracil and oxaliplatin-based chemotherapy [50].

Fig. 5

Proteotranscriptomic characterization of oxaliplatin responder lines in comparison with no-responder ones. (A) Network mapping of 95 proteins showing the 4 sub-clusters composed of more than 3 proteins. Relevant proteins with no associations to others are represented as isolated nodes. Colors are depicted according to the protein abundance (log2ratio) compared to responder PDOs (left panel). Bar-chart of GO terms represented as percentage of annotated proteins using the same color coding (right panel). (B) Hierarchical clustering heatmap of differentially expressed genes per transcriptomic. (C) GSEA hallmarks analysis of RNA-seq data. NES: normalized enrichment score. FDR: false discovery rate. R: responder; NR: non responder

To gain insight into alternative mechanisms of oxaliplatin resistance, a functional protein–protein interaction network of differentially expressed proteins was evaluated by STRING database (http://string-db.org, version 11). Using a high confidence score, some interesting interactions were identified, with some subclusters remaining at the highest confidence score (Fig. 5A). Gene ontology (GO) analysis revealed an increase of proteins related with translation in non-responder models, with an interesting subcluster enriched in tRNA-aminoacylation process that included different aminoacyl-tRNA synthetase (ARSs) and an auxiliary protein AIMP2, joined to form the cytoplasmic multi-tRNA synthetase complex [MSC] [51]. Beyond their role in protein synthesis, ARSs were related to the induction of unfolded protein response, which in turn was associated with escape from chemotherapy induced senescence [52].

Another enriched process in non-responder models was mitochondrial import, including among other proteins two subunits of the ATP synthase (ATP5A and ATP5B), a mitochondrial membrane protein complex mediating ATP synthesis, and citrate synthase (CS), one of the key enzymes in the tricarboxylic acid cycle (TCA), indicating a switch from a glycolytic based to a mitochondrial metabolism with oxidative phosphorylation (OXPHOS) as main source of energy. This could represent a way for some cancer cells to repair platinum-induced DNA damage more efficiently [53].

In the RNA-seq analysis, 597 differentially expressed genes were detected (Fig. 5B). GSEA study showed that non-responding PDOs were enriched in G2M checkpoint, TGFbeta and DNA repair hallmarks (Fig. 5C), somewhat expected as regards oxaliplatin response. Indeed, it acts by inducing DNA adducts, therefore an increase in the capacity to repair damaged DNA could help the cancer cell to survive. A key moment of DNA repair occurs during G2 to M phase transition. In addition, TGFbeta contributes to oxaliplatin resistance in CRC [54]. We noted that non-responder PDOs showed enriched unfolded protein response and PI3K-Akt-mTOR and mTORC1 signaling hallmarks (Fig. 5C). As previously indicated, the former could be involved in the escape from chemotherapy induced senescence [58].

After observing high sensitivity to palbociclib in RTO7 (unpaired t-test p = 0.0002), and to understand the molecular mechanisms involved in this, a comparison of protein expression profiles was made, revealing 245 significantly changing proteins in RTO7 versus non-responder models (RTO2 and mCTO66S3) (Supplementary file S1). Functional protein–protein interaction networks were evaluated using the STRING database with a highest confident score (Fig. 6A). Gene ontology (GO) functional enrichment analysis matched the largest identified cluster (55 of 111 proteins) to biological processes related to protein synthesis, folding and degradation, and with mRNA splicing process via spliceosome, with most of the proteins upregulated in our RTO7 PDO as a palbociclib- responder model (Supplementary Fig. S5). We identified a subcluster comprising all the eight proteins that form the T-complex protein Ring Complex (TRiC) (Supplementary Fig. S5D), an essential eukaryotic molecular chaperonin that aids in the folding of ~ 10% of the proteome including oncoproteins [55].

Fig. 6

Proteotranscriptomic characterization of Palbociclib responder line in comparison with no-responder ones. (A) Network mapping of 111 proteins showing the 4 clusters composed of more than 3 proteins. Proteins with no associations to others were removed and nodes are colored according to the protein abundance (log2ratio) compared to responder PDO (left panel). Bar-chart of GO terms represented as percentage of annotated proteins using the same color coding (right panel). (B) Hierarchical clustering heatmap of differentially expressed genes per transcriptomic. (C) GSEA hallmarks analysis of RNA-seq data. NES: normalized enrichment score. FDR: false discovery rate. R: responder; NR: non responder

Other top processes included in the remaining clusters were related to regulation of cellular component organization, cell adhesion, and metabolic processes, among others (Fig. 6A). All these processes point towards a high proliferative phenotype supported by the up regulation of different mechanisms that preserve proteome integrity.

Transcriptional-wide changes were also observed in palbociclib sensitive model (Fig. 6B). GSEA analysis showed that RTO7 is enriched in MYC targets hallmarks, fatty acids and lipid biosynthesis (which together are a consequence of MYC activation [56]) and unfolded protein response (Fig. 6C), a process that indicates an alteration of protein homeostasis. In the latter hallmark, NPM1, a c-MYC activator [57], is one of the most significantly enriched genes. We could not quantify Myc protein, by SWATH-MS analysis, maybe due to relative scarcity being it a transcription factor. Nevertheless, we were able to detect higher levels of NPM1 among non-clustering proteins.

We observed a weak correlation between protein abundance and transcript quantification among differentially expressed proteins in both comparisons (palbociclib: r = 0.1477, p < 0.01; oxaliplatin: r = 0.1446, p = 0.02; Supplementary Fig. S6A), indicating the occurrence of complex post-transcriptional regulation. For instance, the relevant proteomic data regarding TRiC and ARSs roles in palbociclib responder and oxaliplatin non-responders’ models respectively, were not confirmed by RNA-seq (Supplementary Figure S6B). Nevertheless, some relevant processes related to protein folding, biosynthesis and proliferative features were captured by both RNA and protein analysis.

Integrative functional network (IFNA) analysis of proteotranscriptomic data

To gain deeper insight into the processes involved in the differential drug response in PDOs we performed an integrative analytical approach [58]. Thus, we fused the transcriptomics and proteomics differentially expressed datasets and extracted the common functional context through an integrative network analysis by STRING application via Cytoscape.

Data fusion and integrative functional network analysis (IFNA) were conducted for both oxaliplatin and palbociclib drug response comparisons. The functional clusters shown in the integrative network, representing a high level of interaction and integration of the proteomic and transcriptomic data, are referred as modules in the rest of the manuscript. High confidence interacting proteins are shown in Fig. 7 and Fig. 8, where relevant modules are marked with a shading. In both cases, large networks containing approximately 30–45% of the differential proteins were observed. Interestingly, most of the proteomic clusters and RNA-seq hallmarks (Fig. 5B and Fig. 6B) are now contained in the networks, including some previously isolated clusters. The functional enrichment analysis of the different modules observed in IFNA are listed in Additional File 2 for oxaliplatin and palbociclib comparisons.

Fig. 7

Integrative functional network analysis of oxaliplatin RNA/protein fused dataset. High confidence interactions (score 0.7) showing 7 modules. Nodes in red correspond to proteins identified as differentially expressed by RNA-seq only, in blue those identified by proteomics only and in green those identified by both omics. Key for modules: 1: Glutathione metabolism & antioxidant activity; 2: Nitrogen compound metabolism; 3: Translation & tRNA aminoacylation; 4: TCA & Oxidative phosphorylation; 5: Cell cycle regulation; 6: Antigen presentation; 7: Cell adhesion

Fig. 8

Integrative functional network analysis of palbociclib RNA/protein fused dataset. Highest confidence interactions (score 0.9) showing modules. Nodes in red correspond to proteins identified as differentially expressed by RNA-seq only, in blue those identified by proteomics only and in green those identified by both omics. Grey line marks on the left 7 metabolism related modules with the highest integration. Key for metabolic modules: 1: Oxidative phosphorylation; 2: Pentose phosphate pathway, and Glycolysis; 3: Metabolism of nucleotides, and purinergic nucleotide receptor signaling pathway; 4: PPARA signaling pathway; 7: Lipid biosynthetic process; 9: Oxidoreductase activity; 10: Organonitrogen compound metabolic process. Key for modules on the right: 1: Regulation of signal transduction; 2: Regulation of cell communication and motility; 3: Translation; 4: RNA splicing; 5: Transcription and DNA repair; 6: Cell cycle; 7: Membrane trafficking; 8: Proteasome; 9: Protein ubiquitination; 10: Interferon signaling; 11: ER-Golgi vesicle-mediated transport; 12: TRiC/CCT complex. Details of the modules are indicated in Additional file 2

In the context of oxaliplatin response, the whole dataset contained 12 proteins commonly represented by both omics, 83 only by proteomics and the remaining by RNA-seq. Only high confidence interactions are shown in Fig. 7 containing the main network and two isolated clusters. Globally it contained 191 proteins, 50 of them were detected only by protein, 136 only by RNA while 5 were commonly detected.

Within the main network, 5 modules were selected to perform a functional enrichment analysis that confirmed the processes highlighted in individual omics profiles and now appear connected into this new integrative network. For instance, the relevant translation and tRNA aminoacylation processes appear in module 3 which is now complemented with RNA data. Another interesting finding is module 4, where the proteomic cluster mitochondria protein import is now complemented with TCA cycle and Oxidative phosphorylation processes by RNA data. Proteasome proteomic cluster and DNA repair RNA-seq hallmark are now enriched processes in module 5.

Some previously unconnected differential proteins now appear integrated into the network. Such is the case of CDH17 in module 7, and, interestingly, PRDX6 that is now part of module 1 enriched in redox metabolism, which may be involved in oxaliplatin resistance.

Regarding palbociclib sensitivity, the fused dataset contained 172 proteins detected by proteomics, 73 detected by both proteomic and RNA-seq, and the remaining detected only by RNA-seq. Figure 8 shows the large network of the highest confidence interacting proteins in palbociclib comparison. Certain modules, related to splicing,