Consistent analysis of differentially expressed genes across 7 cell types in papillary thyroid carcinoma

Abstract

Single-cell transcriptome sequencing (scRNA-seq) provides a higher resolution of cellular differences than bulk RNA-seq, enabling the dissection of cell-type-specific responses to perturbations in papillary thyroid carcinoma (PTC). However, cellular genomic features are highly heterogeneous and have a large number of genes without any expression signals, which hinders the statistical power to identify differentially expressed genes and may generate many false-positive results. To overcome this challenge, we conducted an integrative analysis on two PTC scRNA-seq datasets and cross-validated consistent differential expression. By combining results from 32 common cell types in the two studies, we identified 31 consistently differentially expressed genes (DEGs) across seven cell types, including B cells, endothelial cells, epithelial cells, monocytes, NK cells, smooth muscle cells, and T cells. Functional enrichment analysis revealed that these genes are important for the adaptive immune response and autoimmune thyroid diseases. The additional disease-free survival analysis also confirmed that these 31 genes significantly affected patient survival time in large scale thyroid cancer cohort. Furthermore, we experimentally validated one of the top consistent DEGs as a potential biomarker gene of PTC epithelial cells, KRT7, which may be a upstream gene for the NF-κB signaling pathway. The result shows that KRT7 may promote thyroid cancer metastasis through the epithelial-mesenchymal transition and NF-κB signaling pathway. In summary, our single-cell transcriptome integration-based approach may provide insights into the important role of NF-κB in the underlying biology of the PTC.

Keywords: EMT; KRT7; NF-κB; Papillary thyroid carcinoma; ScRNA-seq.

Graphical abstract

Fig. 1

The workflow to identify and validate the consistently differentially expressed genes in thyroid cancer. The comprehensive literature review was conducted and filtered out two scRNA-seq datasets related to PTC. The raw matrix files were downloaded to run the quality control and data normalization. By further data scale and regression, we detected the highly variable genes and define the cell clusters. The cell annotation was performed by against human primary cell atlas database. The gene differential expression analysis was conducted by using Seurat for the two datasets. The non-parametric Wilcoxon rank sum tests on the two datasets were conducted to identify those consistent differential expression. A total of 31 consistent DEG markers were identified across 7 cell types. By mapping the 31 key marker genes to the functional and clinical features, we explored the potential prognostic application in thyroid cancer based on survival data. An epithelial marker gene was validated by experimental evidence in thyroid cancer cell lines.

Fig. 2

Cell clustering and annotation for two PTC scRNA-seq datasets. The t-SNE visualization of cells colored by cluster membership for (A) GSE184362 and (B) GSE191288. The cluster ID were marked as numbers in the corresponding t-SNE map. The pie-chart plots depict the proportion of cell types for each cell cluster in (C) GSE184362 and (D) GSE191288.

Fig. 3

The expression heatmap of 31 genes consistently differentially expressed in two PTC scRNA-seq datasets: GSE184362 (A) and GSE191288 (B). The cell types exhibited distinct coloring.

Fig. 4

The violin plots ofKRT7andCLUpresent the expression level across multiple cell types for two PTC scRNA-seq datasets. The expression of the KRT7 and CLU are summarized for (A) GSE184362 and (B) GSE191288.

Fig. 5

The clinical and network feature for the 31 consistently differentially expressed genes. (A) Based on the TCGA thyroid cancer dataset, the tumor mutational burden (TMB) patterns for those samples with and without genetic mutations in the 31 genes. (B) Analysis of disease-free survival for patients with and without mutations in 31 genes. (C) Based on the 31 genes, a functional subnetwork was reconstructed. The distribution of the 31 genes across cell types is depicted in the bar chart on the right. The network node color is identical to the color code used in the bar chart.

Fig. 6

The functional enrichment analysis for the 31 consistently expressed genes in two PTC datasets. (A) the top over-represented gene ontology terms. (B) the IPA network summary for the top over-represented canonical pathways, regulators, and the 31 DEGs. (C) the tissue specificity based on the large-scale gene expression data. (D) the transcription factor regulation enrichment.

Fig. 7

Knockdown ofKTR7undermines cell motility and reverses EMT process. (A) Transwell assay was preformed to show that KD KRT7 attenuated significantly migration abilities of BCPAP. Scale bars, 20 µm. The statistic was analyzed on the right individually. (B) siKRT7 was transfected to BCPAP for 48 h and protein levels of N-CADHERIN, E-CADHERIN, VIMENTIN, P65 were analyzed by Western blot. (C) siKRT7 was transfected to BCPAP for 48 h and mRNA levels of KRT7, P65, VIMENTIN, SNAIL, N-CADHERIN, were analyzed by RT-qPCR. (D and E) The levels of VIMENTIN (D) and F-ACTIN (E) were assessed by immunofluorescence in siKRT7 treated BCPAP cells. Scale bars, 20 µm. Relative fluorescence intensities of VIMENTIN and F-ACTIN were shown on the right. Data in this figure, mean ± SD, *P < 0.05, * *P < 0.01, * **P < 0.001.