Biomaterials serve as the basis of implants, tissue engineering scaffolds, and multiple other biomedical therapeutics. New technologies, such as single cell RNA sequencing (scRNAseq), are enabling characterization of the biomaterial response to an unprecedented level of detail, facilitating new discoveries in the complex cellular environment surrounding materials. We performed scRNAseq and integrated data sets from multiple experiments to create a single cell atlas of the biomaterials response that contains 42,156 cells from biological extracellular matrix (ECM)-derived and synthetic polyester (polycaprolactone, PCL) scaffold biomaterials implanted in murine muscle wounds. We identified 18 clusters of cells, including natural killer (NK) cells, multiple subsets of fibroblasts, and myeloid cells, many of which were previously unknown in the biomaterial response. To determine intra and intercellular signaling occurring between the numerous cell subsets, including immune-stromal interactions in the biomaterial response, we developed Domino (github.com/chris-cherry/domino), a computational tool which allows for identification of condition specific intercellular signaling patterns connected to transcription factor activation from single cell data. The Domino networks self-assembled into signaling modules and cellular subsets involved in signaling independent of clustering, defining interactions between immune, fibroblast, and tissue-specific modules with biomaterials-specific communication patterns. Further compilation and integration of biomaterials single cell data sets will delineate the impact of materials chemical and physical properties and biological factors, such as anatomical placement, age, or systemic disease, that will direct biomaterials design.

Abstract Background Tumor response to therapy is affected by both the cell types and the cell states present in the tumor microenvironment. This is true for many cancer treatments, including notably immune checkpoint inhibitors (ICIs). While it is well-established that ICIs promote T cell activation, their broader impact on other intratumoral immune cells is unclear; this information is needed to identify new mechanisms of action and improve ICI efficacy. Many preclinical studies have begun to use single cell analysis to delineate therapeutic responses in individual immune cell types within tumors. One major limitation to this approach is that therapeutic mechanisms identified in preclinical models have failed to fully translate to human disease, restraining efforts to improve ICI efficacy in bench to bedside research. Method We previously developed a computational transfer learning approach to identify shared biology between independent high-throughput single-cell RNA sequencing (scRNA-seq) datasets. In the present study, we test this framework’s ability to identify conserved and clinically relevant transcriptional changes in complex tumor scRNA-seq data and further expand its application beyond comparison of scRNA-seq datasets into comparison of scRNA-seq datasets with additional data types such as bulk RNA-seq and mass cytometry. Results We found a conserved signature of NK cell activation in anti-CTLA-4 responsive mice and human tumors. In human melanoma, we found that the NK cell activation signature correlates with longer overall survival and is predictive of anti-CTLA-4 (ipilimumab) response. Additional molecular approaches to confirm the computational findings demonstrated that human NK cells express CTLA-4 and bind anti-CTLA-4 independent of the antibody binding receptor (FcR), and that similar to T cells, CTLA-4 expression by NK cells is modified by cytokine-mediated and target cell-mediated NK cell activation. Conclusions These data demonstrate the ability of our transfer learning approach to identify cell state transitions conserved in preclinical models and human tumors. This approach can be adapted to explore many immuno-oncology questions, enhancing bench to bedside research and enabling better understanding and treatment of disease. Graphical Abstract

Abstract Non-negative matrix factorization (NMF) is an unsupervised learning method well suited to high-throughput biology. Still, inferring biological processes requires additional post hoc statistics and annotation for interpretation of features learned from software packages developed for NMF implementation. Here, we aim to introduce a suite of computational tools that implement NMF and provide methods for accurate, clear biological interpretation and analysis. A generalized discussion of NMF covering its benefits, limitations, and open questions in the field is followed by three vignettes for the Bayesian NMF algorithm CoGAPS (Coordinated Gene Activity across Pattern Subsets). Each vignette will demonstrate NMF analysis to quantify cell state transitions in public domain single-cell RNA-sequencing (scRNA-seq) data of malignant epithelial cells in 25 pancreatic ductal adenocarcinoma (PDAC) tumors and 11 control samples. The first uses PyCoGAPS, our new Python interface for CoGAPS that we developed to enhance runtime of Bayesian NMF for large datasets. The second vignette steps through the same analysis using our R CoGAPS interface, and the third introduces two new cloud-based, plug-and-play options for running CoGAPS using GenePattern Notebook and Docker. By providing Python support, cloud-based computing options, and relevant example workflows, we facilitate user-friendly interpretation and implementation of NMF for single-cell analyses.

Abstract Advanced age is strongly correlated with both increased cancer incidence and general immune decline. The immune tumor microenvironment (ITME) has been established as an important prognostic of both therapeutic efficacy and overall patient survival. Thus, age-related immune decline is an important consideration for the treatment of a large subset of cancer patients. Current studies of aging-related immune alterations are predominantly performed on non-cancerous tissue, requiring additional study into the effects of age on tumor immune infiltration. We leverage large scale transcriptional data sets from The Cancer Genome Atlas and the Genotype-Tissue Expression project to distinguish normal age-related immune alterations from age-related changes in tumor immune infiltration. We demonstrate that while there is overlap between the normal immune aging phenotype and that of the ITME, there are several changes in immune cell abundance that are specific to the ITME, particularly in T cell, NK cell, and Macrophage populations. These results suggest that aged immune cells are more susceptible to tumor suppression of cytotoxic immune cell infiltration and activity than normal tissues, which creates an unfavorable ITME in older patients in excess of normal immune decline with age and may inform the application of existing and emerging immunotherapies for this large population of patients. We additionally identify that age-related increases in tumor mutational burden are associated with decreased DNA methylation and increased expression of the immune checkpoint genes PDL1, CD80, and LAG3 which may have implications for therapeutic application of immune checkpoint blockade in older patients.

Abstract Natural killer (NK) cells play a critical role in physiologic and pathologic conditions such as pregnancy, infection, autoimmune disease and cancer. In cancer, numerous strategies have been designed to exploit the cytolytic properties of NK cells, with variable success. A major hurdle to NK-cell focused therapies is NK cell recruitment and infiltration into tumors. While the chemotaxis pathways regulating NK recruitment to different tissues are well delineated, the mechanisms human NK cells employ to physically migrate are ill-defined. We show for the first time that human NK cells express fibroblast activation protein (FAP), a cell surface protease previously thought to be primarily expressed by activated fibroblasts. FAP degrades the extracellular matrix to facilitate cell migration and tissue remodeling. We used novel in vivo zebrafish and in vitro 3D culture models to demonstrate that FAP knock out and pharmacologic inhibition restrict NK cell migration, extravasation, and invasion through tissue matrix. Notably, forced overexpression of FAP promotes NK cell invasion through matrix in both transwell and tumor spheroid assays, ultimately increasing tumor cell lysis. Additionally, FAP overexpression enhances NK cells invasion into a human tumor in immunodeficient mice. These findings demonstrate the necessity of FAP in NK cell migration and present a new approach to modulate NK cell trafficking and enhance cell-based therapy in solid tumors. Graphical Abstract

Abstract The etiology of diseases driven by dysregulated mRNA metabolism can be elucidated by characterizing the responsible RNA-binding proteins (RBPs). Although characterizations of RBPs have been mainly focused on their binding sequences, not much has been investigated about their preferences for RNA structures. We present nearBynding, an R/Bioconductor pipeline that incorporates RBP binding sites and RNA structure information to discern structural binding preferences for an RBP. nearBynding visualizes RNA structure at and proximal to sites of RBP binding transcriptome-wide, analyzes CLIP-seq data without peak-calling, and provides a flexible scaffold to study RBP binding preferences relative to diverse RNA structure data types.

Abstract Senescent cells (SnCs) contribute to normal tissue development and repair but accumulate with aging where they are implicated in a number of pathologies and diseases. Despite their pathological role and therapeutic interest, SnC phenotype and function in vivo remains unclear due to the challenges in identifying and isolating these rare cells. Here, we developed an in vivo -derived senescence gene expression signature using a model of the foreign body response (FBR) fibrosis in a p16 Ink4a - reporter mouse, a cell cycle inhibitor commonly used to identify SnCs. We identified stromal cells (CD45 - CD31 - CD29 + ) as the primary p16 Ink4a expressing cell type in the FBR and collected the cells to produce a SnC transcriptomic signature with bulk RNA sequencing. To computationally identify SnCs in bulk and single-cell data sets across species and tissues, we used this signature with transfer learning to generate a SnC signature score (SenSig). We found senescent pericyte and cartilage-like fibroblasts in newly collected single cell RNAseq (scRNASeq) data sets of murine and human FBR suggesting populations associated with angiogenesis and secretion of fibrotic extracellular matrix, respectively. Application of the senescence signature to human scRNAseq data sets from idiopathic pulmonary fibrosis (IPF) and the basal cell carcinoma microenvironment identified both conserved and tissue-specific SnC phenotypes, including epithelial-derived basaloid and endothelial cells. In a wound healing model, ligand-receptor signaling prediction identified putative interactions between SnC SASP and myeloid cells that were validated by immunofluorescent staining and in vitro coculture of SnCs and macrophages. Collectively, we have found that our SenSig transfer learning strategy from an in vivo signature outperforms in vitro -derived signatures and identifies conserved and tissue-specific SnCs and their SASP, independent of p16 Ink4a expression, and may be broadly applied to elucidate SnC identity and function in vivo .

Abstract Motivation Genomics features with similar genomewide distributions are generally hypothesized to be functionally related, for example, co-localization of histones and transcription start sites indicate chromatin regulation of transcription factor activity. Therefore, statistical algorithms to perform spatial, genomewide correlation among genomic features are required. Results Here, we propose a method, StereoGene, that rapidly estimates genomewide correlation among pairs of genomic features. These features may represent high throughput data mapped to reference genome or sets of genomic annotations in that reference genome. StereoGene enables correlation of continuous data directly, avoiding the data binarization and subsequent data loss. Correlations are computed among neighboring genomic positions using kernel correlation. Representing the correlation as a function of the genome position, StereoGene outputs the local correlation track as part of the analysis. StereoGene also accounts for confounders such as input DNA by partial correlation. We apply our method to numerous comparisons of ChIP-Seq datasets from the Human Epigenome Atlas and FANTOM CAGE to demonstrate its wide applicability. We observe the changes in the correlation between epigenomic features across developmental trajectories of several tissue types consistent with known biology, and find a novel spatial correlation of CAGE clusters with donor splice sites and with poly(A) sites. These analyses provide examples for the broad applicability of StereoGene for regulatory genomics. Availability The StereoGene C++ source code, program documentation, Galaxy integration scripts and examples are available from the project homepage http://stereogene.bioinf.fbb.msu.ru/ Contact favorov@sensi.org Supplementary information Supplementary data are available online.

Identifying potential mechanisms of resistance while tumor cells still respond to therapy is critical to delay acquired resistance. We generated the first comprehensive multi-omics, bulk and single cell data in sensitive head and neck squamous cell carcinoma (HNSCC) cells to identify immediate responses to cetuximab. Two pathways potentially associated with resistance were focus of the study: regulation of receptor tyrosine kinases through the transcription factor TFAP2A, and epithelial-to-mesenchymal transition (EMT) process. Single cell RNA-seq demonstrates heterogeneity, with cell specific TFAP2A and VIM expression profiles in response to treatment. RNA-seq and ATAC-seq reveal global changes within five days of cetuximab therapy, suggesting early onset of mechanisms of resistance; and corroborates cell line heterogeneity, with different TFAP2A targets or EMT markers affected by therapy. Lack of TFAP2A reduces HNSCC growth and is enhanced by cetuximab and JQ1. Regarding the EMT process, short term cetuximab therapy has the strongest effect on inhibiting migration. TFAP2A silencing does not affect cell migration, supporting an independent role for both mechanisms in resistance. Overall, we show that immediate adaptive transcriptional and epigenetic changes induced by cetuximab are heterogeneous and cell type dependent; and independent mechanisms of resistance arise while tumor cells are still sensitive to therapy.

Motivation: Current bioinformatics methods to detect changes in gene isoform usage in distinct phenotypes compare the relative expected isoform usage in phenotypes. These statistics model differences in isoform usage in normal tissues, which have stable regulation of gene splicing. Pathological conditions, such as cancer, can have broken regulation of splicing that increases the heterogeneity of the expression of splice variants. Inferring events with such differential heterogeneity in gene isoform usage requires new statistical approaches. Results: We introduce Splice Expression Variability Analysis (SEVA) to model increased heterogeneity of splice variant usage between conditions (e.g., tumor and normal samples). SEVA uses a rank-based multivariate statistic that compares the variability of junction expression profiles within one condition to the variability within another. Simulated data show that SEVA is unique in modeling heterogeneity of gene isoform usage, and benchmark SEVA's performance against EBSeq, DiffSplice, and rMATS that model differential isoform usage instead of heterogeneity. We confirm the accuracy of SEVA in identifying known splice variants in head and neck cancer and perform cross-study validation of novel splice variants. A novel comparison of splice variant heterogeneity between subtypes of head and neck cancer demonstrated unanticipated similarity between the heterogeneity of gene isoform usage in HPV-positive and HPV-negative subtypes and anticipated increased heterogeneity among HPV-negative samples with mutations in genes that regulate the splice variant machinery. Conclusion: These results show that SEVA accurately models differential heterogeneity of gene isoform usage from RNA-seq data. Availability: SEVA is implemented in the R/Bioconductor package GSReg.