RNA-Sequencing, Differentially Expressed Gene (DEG), Ingenuity Pathway Analysis (IPA), and Enrichment Analyses

cDNA libraries were prepared using the TrueSeq RNA Sample Preparation Kit (Illumina) from RNA samples isolated from eight non-SS controls and nine pSS patients obtained from the OMRF (see Figure 1A and Table 1 for sample ID numbers and additional patient information). The cDNA libraries were then sequenced on an Illumina NovaSeq sequencer (50 cycle paired-end). After initial quality control metrics were determined using FASTQC v0.4.3, the raw reads were then mapped to the reference genome (Hg38 build) with Hisat2 (25) v2.1.0. using Bowtie (26) v2.2.6 as the underlying aligner. Reads mapping uniquely to each gene in the reference genome were then quantified using featureCounts (27) from the Subread (28) package. The resulting counts matrix was imported into R for read count normalization and differential gene expression analysis using the DESeq2 (29) package. An adjusted p-value < 0.05 based on Benjamini-Hochberg method was used as a cut-off for determining differentially expressed genes. Gene ontology analysis for identification of enriched pathways was performed using the Database for Analysis Visualization and Integrated Discovery (DAVID) v6.8 (30, 31). The databases were queried by providing the list of official names of genes of interest. The resulting KEGG Pathway (32–35) table was then imported into R for generating bar plots. ClueGO v2.5.7 (36) add-on of the Cytoscape Platform v3.72 (37) was used to generate Figures 1C, D . For Figure 5A , a table containing official names and corresponding foldchange values of the 80 genes that are common across all three datasets was used as input for the Ingenuity Pathway Analysis (IPA) software (Qiagen). The core analysis function of the software was used to interpret the data and generate the graph ( Figure 5A ) of canonical pathways that correspond to the changes in gene expression. Figures 5B–D were generated from the genes enriched in the indicated canonical pathways using STRING v11 (38–41), and the clusters were defined based on likelihood of protein-protein interaction using the STRING’s implementation of the Markov cluster algorithm (MCL) clustering (42). Sequencing data generated for this study has been deposited in the Gene Expression Omnibus (GEO) database under the accession number {"type":"entrez-geo","attrs":{"text":"GSE157159","term_id":"157159"}}GSE157159.

Comprehensive Transcriptomic Analysis of Minor Salivary Glands of Patients with Sjögren’s Syndrome. (A) Plot demonstrates principal component analysis (PCA) coordinates for each minor salivary gland for eight non-pSS controls and nine pSS patient samples. Blue and red circles represent pSS and non-pSS, respectively. (B) Heatmap visualization of the top 100 differentially expressed genes (DEGs) in minor salivary glands of non-SS controls and pSS patient samples. Network visualization of enriched biological processes of the top 100 upregulated (C) and top 100 downregulated (D) DEGs in pSS.

Gene Networks and Canonical Pathways Generated by Ingenuity Pathway Analysis (IPA). (A) Top 21 significantly affected canonical pathways for DEGs identified in our dataset based on IPA. Bars denote the different pathways based on Z-scores (2.0 < Z score < −2.0). Orange color indicates pathway activation, while blue denotes suppression. (B–D) Gene network analysis using IPA for pathways selected from panel A above. Solid lines denote high correlation between members of a defined cluster and dotted lines denote weak inter-cluster correlation.