{'fraction_non_alphanumeric': 0.03533026113671275, 'fraction_numerical': 0.016732499451393462, 'mean_word_length': 5.14391641388608, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Identification of alternative splicing and lncRNA genes in pathogenesis of small cell lung cancer based on their RNA sequencing

Funding sources The study was supported by the Scientific Research Projects of Institutions of Medical and Health Institutions in Yunnan Province The Role of NOD-like Receptors and Inflammatory Bodies in the Development of Xuanwei Lung Cancer (project No. 2016NS017); Number of Training Objects for Leaders of Medical Disciplines in Yunnan Province (project No. D2017013); Doctoral Research Fund of the First Affiliated Hospital of Kunming Medical University The research about PDGFRB functions on lung squamous cell carcinoma progression and its potential usage as a clinical lung squamous cell carcinoma maker.


Introduction
Small-cell lung cancer (SCLC) is an aggressive neuroendocrine carcinoma. Lung cancer is the leading cause of cancer-related deaths worldwide. About 15% of lung cancer cases are SCLC cases. 1 Its prognosis is poor, and has a trend toward metastasizing to the brain at an early stage, despite advances in treatment modalities that have improved survival. 2 Histological diagnosis can provide information valuable in planning the treatment program, but it is insufficient for predicting clinical outcomes. 3 At present, suitable molecular biomarkers for SCLC diagnosis are needed to assess the patient's prognosis.
Recently, noncoding RNAs (ncRNAs) have received considerable attention due to their key role in affecting shared gene expression programs. 4 These RNAs are usually microRNAs or long-non-coding RNAs (lncRNAs). With the development of high-throughput sequencing technology, ncRNAs have been discovered in a wide range of biological processes. However, the functions of lncRNAs in SCLC are not well characterized and the identification of lncRNA biomarkers is challenging. Typically, lncRNA functions are characterized using a "guilt by association" strategy. 5 Recent studies have demonstrated that lncRNAs compete with endogenous RNAs (ceRNAs) in regulating expression level of other transcripts. 6,7 For example, the lncRNA HULC acts as an endogenous ceRNA and plays an important regulatory role in lung cancer. 8 In addition, a recent study using network analysis identified a lncRNAassociated ceRNA network (across 12 different cancers) and prognostic markers of lncRNAs. 7 Another study proposed that the ceRNA interaction network of SCLC could reveal canonical oncogenic pathways. 9 Collectively, this data underscores the importance of lncRNAs interaction with ceRNAs and indicates that the integration of expression profiles with network analysis could identify the risk of lncRNAs and the underlying tumor pathology.
Current studies have suggested that the alternative splicing of the RNAs, as an essential layer of gene expression regulation that expands the cell proteome, plays complex and critical role in breast cancer. 10 Additionally, the alternative splicing of EGFR, CD44, PIK3C3, RRAS2, MAP-KAP1, and FGFR2 was shown to raise the risk of stratification in non-small cell lung cancer (NSCLC) patients. 11 Alternative splicing of S6K1 gene can regulate the mammalian target of rapamycin (mTOR) complex 1 signaling pathway. This has been strongly linked with the onset and progression of NSCLC. Alternative placing of this gene may, therefore, constitute a relevant and innovative biomarker. The short isoforms of S6K1 are activated mTORC1, leading to increased 4E-BP1 phosphorylation. 12,13 In this research, a multi-step computational approach was used for constructing a functional lncRNAs-mediated ceRNA network (LMCN) in SCLC. We utilized the highthroughput profiles of SCLC samples, which were collected from the ArrayExpress database. Intergradations of alternative splicing, genome-wide lncRNA/mRNA expression profiles, comprehensive miRNA-target interactions, and functional analyses were entailed in the systematic analysis. Based on LMCN, we found that several lncRNAs exhibited specific topological features, which are consistent with coding mRNAs regulatory associated with SCLC pathology. The alternative RNA splicing involved in RNA processing significantly related with cell cycle functions was also identified. Additionally, by exploring a competitive sub-network with in the LMCN we determined that MCM3AP-AS and MIR17HG were involved in a ceRNA regulatory module and synergistically competed with MATR3, XPO1 and ZCCHC14. These analyses demonstrated that the integration of alternative splicing profile with cancer-associated lncRNA-mediated ceRNA networks could be used to promote the discovery of molecular biomarkers and SCLC therapeutics.
For data preprocess, Illumina reads were trimmed for adaptor sequence, and low-complexity or low-quality sequence was removed, then mapped to Human genome reference hg19 via http://bowtie-bio.sourceforge.net/bowtie2/ index.shtml (data was normalized using DESeq in R package to analyse count data for differential expression(https:// bioconductor.org/packages/release/bioc/html/DESeq. html). For data analysis, RNA-seq data was analyzed using the set of open source software programs of the Tuxedo suite: TopHat v. 1.3.3 and Cufflinks (http://ccb.jhu.edu/ software/tophat/tutorial.shtml and http://cole-trapnelllab.github.io/cufflinks/). 14,15 The R platform available as free software (http://cran.r-project.org/) was applied. In brief, PolyAplus RNA (RNA transcripts) was purified from 100 ng of total RNA with oligo-dT beads, and then the samples were fragmented with divalent cations under high temperature. First strand synthesis was performed with random hexamer and reverse transcriptase, and 2 nd strand synthesis with RNase H and DNA polymerase I. Following cDNA synthesis, the double-stranded products were end-repaired, a single "A" was added and then the Illumina Paired End Adapters (https://support.illumina. com/content/dam/illumina-support/documents/documentation/chemistry_documentation/experiment-design/ illumina-adapter-sequences-1000000002694-09.pdf) were ligated on to the cDNA products. The ligation products were purified using gel electrophoresis. The target size range for these libraries was ~300 bp, such that the final library for sequencing would have cDNA inserts with sizes of ~200 bp long. One run of 2 × 50 bp paired-end sequencing was performed on the HiSeq2000 instrument, using 1 lane per tissue, to produce approx. 80 million read pairs per tissue.

Reconstructing transcriptomes
To determine splice variants within each tissue, we aligned reads to the hg19 genome using TopHat v. 1.3.3 with an intron database, which extracted from the UCSC Known Genes data set (http://genome.ucsc.edu). Aligned reads were then assembled into transcript fragments using Cufflinks software (with parameters '-F 0.05'). Then, the Cuffcompare (http://cufflinks.cbcb.umd.edu/) was used to compare these transfrags to the Ensembl annotation. After eliminating likely artifacts and assemblies not associated with Ensembl genes, the Cuffdiff was used to redistribute reads with a high-confidence transcript sets obtained. 15,16 Cuffcompare classifies assembled transcripts into multiple categories: the known coding genes, the known noncoding genes, the potentially novel transcripts, and the undefinable transcripts. The potentially novel transcripts and the undefinable transcripts were used for downstream analysis. FPKM expression level for this set was then re-estimated from the original alignments using Cuffdiff.

Discovery of alternative splicing events
To determine alternative splicing events in the LMCN, we developed a software package that can identify the average shortest (AS) profile and analyse the exons-involved in 1 transcript in the transcripts profiles. We restricted the analysis to Ensembl genes with FPKM ≥ 0.1 re-estimated using Cuffdiff as described above. We define the exon skipping event as the pair between exon-containing splice form and exon-skipping splice form, where the boundaries of flanking introns are required to match precisely. To determine the novel events, the exons and spanning introns were compared against the annotation data sets (UCSC Genes and Ensembl).

Comparison of alternative splicing events among groups
For each event, we calculated the exon inclusion ratio R = FPKMon/(FPKMon + FPKMoff) for each tissue. To account for minor differences in the annotation of splice junctions, the differentially expressed level of average shortest (AS) mRNA was calculated with limma package in R. All differentially expressed mRNA across groups data were further screened by limma package based on Bayes method (http://www.bioconductor.org/). After processing using the top table method of limma package, the differentially expressed genes were analyzed. Then, they were compared with AS genes. The AS genes and DEG were selected for further network assay. The procedure was as follows: lncRNA, miRNA and mRNA expression profiles were obtained from starBase v. 2.0 (http://starbase.sysu.edu.cn/starbase2/index. php). Interactions between miRNA-lncRNAs were identified using existing miRNA target prediction methods. Putative miRNA-lncRNA interactions were identified using the miRanda algorithm with the default parameters. The miRNA target binding sites were predicted on entire lncRNAs. The data of the miRNA-mRNA interactions was downloaded from miRNA reference databases, the TarBase (v. 6, http://diana.imis.athena-innovation.gr/DianaTools/ index.php?r=tarbasev6/index) and the mirTarBase (v. 4.5, http://mirtarbase.mbc.nctu.edu.tw/php/download.php). Both databases store manually curated collections of experimentally supported miRNA targets.

Functional enrichment analysis
The gene ontology based on gene ontology database and the pathway enrichment analysis based on the database of KEGG pathway were carried out using Fisher's test. P-value <0.05 and false discovery rate (FDR) <0.05 were considered to be statistically significant.

Integrated analysis of gene signatures
After preprocessing the expression profile of mRNA datasets, we identified 3,105 differentially expressed genes of SCLA using limma package in R language. Limma is an R/Bioconductor software package that provides an integrated solution for analyzing data from gene expression experiments. Over the past decade, limma has been a popular choice for gene discovery through differential expression analyses of microarray and high-throughput PCR data. 16 Figure 1A shows the log fold changes (log FC) in top 10 gene signatures. Figure 1B shows the volcano plot of these genes (log odds). Log FC value (Fig. 1A) indicated that the upregulated genes were SNAP25 and MANEAL; the downregulated genes were AGER, EPAS1, EMP2, RP 11-114H24.5, HBA2, NFKBIA, ABCA3, and CA4, and they were considered as suppressor genes in SCLC ( Table 1).

Identification of the cluster in LMCN
While the LMCN could provide a global view of all possible competing ceRNA interactions that could be used to investigate the regulatory properties of the lncRNAs, the partial sub-networks revealed a more detailed pattern of how  Log FC -logarithmic fold change; AveExpr -the average expression; t -two-sample t-test; B -potential scenario for the statistical test for differential expression. A the lncRNAs synergized with competing mRNAs. We derived a cluster from the LMCN by applying a Pearson's correlation coefficient threshold >0.9 (Fig. 3). The genes with a top degree were EPB41L4A-AS1, HOXA-AS2, XIST, DLEU2, FGD5-AS1, ALMS1-IT1, SNHG12, MIR17HG, MIR4720, and SCARNA10 (Table 2).

Functional enrichment analysis
Functional enrichment analysis of the obtained gene signatures was carried out, including Gene Ontology (GO) term enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway. The GO analysis was carried out in 3 categories, including biological processes (BP), molecular functions (MF) and cellular components (CC). The results of GO analysis showed that cluster genes were significantly enriched in vacuolar transport, Golgi vesicle transport, mitotic cell cycle phase transition, endosomal transport, and cell cycle phase transition. After the false discovery rate was calculated, the most significant 10 GO terms including BP, MF and CC were presented in Table 3.
For KEGG pathway enrichment, 4 terms were obtained and the most significant ones were Olfactory transduction

Identification of the hub genes with average shortest events in LMCN
The genes in the LMCN-shared high degree were HOXA-AS2, EPB41L4A-AS1, DLEU2, XIST, MIR17HG,   FGD5-AS1, SNHG12, ALMS1-IT1, ZRANB2-AS2, and MIR4720. In the obtained cluster, the average rank value based on p-value was calculated and genes with the highest value were considered as hub genes. The result was shown in Table 2: genes in a cluster with an average degree value of 6.6 were connected significant. This is consistent with our earlier AS analysis. The average shortest events (AS events) present in several genes (NEAT1, FOXP1, STX12, KLF7, ITPR1) in the cluster of cell cycle pathways further confirmed the significance of hub genes.

Pathway analysis of small-cell lung cancer pathogenesis
In the current study, based on LMCN and pathway enrichment of gene signatures, several KEGG pathways were found. Some of these pathways have already been reported by other researchers, such as cell cycle and endocytosis. [17][18][19] One of the significant pathways was endocytosis with the p-value of 3.49E-08. Endocytosis is a form of active transport and an important process for the regulation of cell signal transduction, cellular dynamics and cell-to-cell communication. Endocytosis and its counterpart, exocytosis, take place in all cells, which include pinocytosis and phagocytosis. The current study showed that in this pathway ligand-induced endocytosis requires NM IIA to induce the degradation of EGFR, causing the endocytosis process. 20 Our present study verified the endocytosis process in cancer and predicted that the interaction of ALMS1-IT1 and VPS37B as a novel mechanism in trans-Golgi network (Fig. 4). VPS37B is located in chromosome 12. This gene binds to TSG101 at several different sites, which function in vacuolar protein sorting. Several other genes are also in the pathway list, e.g., GIT2, ZFYVE9, HSPA6, RAB11FIP4, CHMP2B, TGFB2, RAB4A, RAB5A, RK5, and others. Some of them are lncRNAs, such as XIST, TRAF3IP2-AS1, DLEU2, DLEU2, MIR4720, and SNHG12. These genes were analyzed via the LMCN, indicating that they were involved in all the 4 main functions (receptor-mediated endocytosis, caveolae, macropinocytosis, and phagocytosis). They were considered as the critical genes. Therefore, endocytosis pathway as well as these genes was thought to be very important in the pathogenesis of SCLC. The map of OXPHOS pathway was shown in Fig. 4.

Critical genes in pathogenesis of SCLC
The genes in the cluster were considered as critical genes for the pathogenesis of SCLC. These results strongly support SRSF1 as a prognostic biomarker in SCLC and provide a rationale for personalized therapy in SCLC. The critical genes with AS events in this study included NEAT1, GO ID -Gene Ontology identification; Bg ratio -background ratio indicates the number of target genes regulated by the transcription factor; Count -indicates how many reads map to each feature. FOXP1, STX12, KLF7, and ITPR1. Among these genes, NEAT1 and FOXP1 were reported to be associated with cancer. Nuclear-enriched abundant transcript 1 (NEAT1, non-protein coding) produces a lncRNA, which is transcribed from the multiple endocrine neoplasia locus. This lncRNA is retained in the nucleus where it forms the core structural component of the paraspeckle sub-organelles. Its high expressions are associated with poor prognosis in different types of solid tumors, such as lung cancer, esophageal cancer, colorectal cancer, and hepatocellular carcinoma. 21 It was aberrantly upregulated and may act as a transcriptional regulator for numerous genes, including some genes involved in unfavorable prognosis in cancer patients. 22 Increased serine/arginine splicing factor 1 (SRSF1) DNA copy number and mRNA overexpression has been linked to poor prognosis for SCLC. Functional studies in vitro and in vivo demonstrate that SRSF1 may play a central role in DNA repair and chemosensitivity. 23 On the other hand, downregulated NEAT1 promotes leucocyte differentiation in acute promyelocytic leukemia. Our study examined the relationship between NEAT1 and SCLC and confirmed that the defection of this gene was closely related to the pathogenesis of SCLC.
FOXP1 gene belongs to the subfamily P of the Forkhead box (FOX) transcription factor family. Forkhead box transcription factors play an important role in the regulation of tissue-and cell type-specific gene transcription during both development and adulthood. This gene may express in several types of human malignant tumors. It also has been reportedly associated with the metastasis and patient prognosis in NSCLC. 24 The expression of FoxP1 in NSCLC was significantly higher. Since the FOXP1 protein contains both DNA-binding-and protein-protein binding-domains, it was reported to be a regulator in cancer cells. The low expression of FOXP1 was independent factors for predicting the poor prognosis for NSCLC. 24 Alternative splicing of FOXP1 results in multiple transcript variants encoding different isoforms, which is also related to many other cancers. 25 NEAT1 and FOXP1 genes connected closely in the cluster and might be another potential way for the pathogenesis of SCLC. The other genes, STX12, KLF7 and ITPR1, had not been reported to be in any relationship with SCLC and they need to be investigated further.
In this study, a comprehensive process of data in datasets of SCLC for protein-protein interactions (PPI) network was conducted. Some clusters where the genes connected closely to SCLC were selected. Then, the function and signaling pathways for SCLC were presented systematically. Finally, the potential pathways and 2 critical genes in SCLC pathogenesis were identified. These may aid in SCLC therapy. However, this study was based on the bioinformatics reports. These results need to be further proven with laboratory tests.