Acquisition of human MG risk genes and miRNAs
To ensure a high quality, the data of MG risk genes was obtained following previously described criteria [16, 17]. First, we obtained all publications from PubMed database using the keyword “myasthenia gravis” to collect experimentally supported MG-related genes. We then manually curated MG-related genes following the criteria we set: (1) the gene was notably differentially expressed in at least five human MG samples and was detected using reliable low-throughput experiment methods, such as RT-PCR and northern blot; or (2) the gene was contained in several current public databases we mentioned previously [16, 17]. In addition, we updated the list of MG risk miRNAs referred to our previous studies published online [18].
Analysis of PPI and module
Based on the Search Tool for the Retrieval of Interacting Genes (STRING) database [19], PPI network of MG risk genes was established with a confidence score ≥ 0.4 and visualized by means of Cytoscape software. Furthermore, Plug-in Molecular Complex Detection (MCODE) [20] was utilized to identify functionally related and highly interconnected modules from the PPI network with a degree cutoff = 2, node score cutoff = 0.2, k-core = 2 and maximum depth = 100.
Functional enrichment analysis
To screen the potential functions of modules in the PPI network, gene ontology (GO) functions and KEGG pathway enrichment analysis were performed by clusterProfiler package [21], which is an R package of Bioconductor and can perform statistical analysis and visualization of functional clustering of gene sets or gene clusters. P-value < 0.05 was considered to be significantly enriched function annotations.
Identification of miRNA-gene and miRNA-lncRNA interactions
First, the interactions of miRNA-gene were acquired from miRTarBase (release 8.0) [22], which is a database that curates experimentally validated miRNA-gene interactions. miRNA-gene pairs were selected when they were supported by western blot and/or reporter assay showing high-confidence functional. Further, the gene and miRNA in each miRNA-gene pair were filtered by aforementioned human MG risk genes of modules and MG risk miRNAs. Ultimately, we only retained the intersections among human MG risk genes in modules, human MG risk miRNAs and experimentally supported miRNA-gene interactions as the candidate miRNA-gene pairs. Then, the interactions of MG risk miRNA and lncRNA were extracted from starBase (v3.0) [23], which is a database decoding the interactions of miRNAs, lncRNAs, mRNAs from large-scale CLIP-Seq data, such as HITS-CLIP, iCLIP, PAR-CLIP, and CLASH.
Identification of potential gene-lncRNA competing pairs
Hypergeometric test and the Pearson correlation coefficient (PCC) were performed to identify the potential competing gene-lncRNA pairs. First, hypergeometric test was applied for assessing the significance of the shared miRNAs between each gene and lncRNA. The formula is as follows:
$$P=F\left(x|M,K,N\right)=\sum_{i=0}^{x}\frac{\left(\genfrac{}{}{0pt}{}{K}{i}\right)\left(\genfrac{}{}{0pt}{}{M-K}{N-i}\right)}{\left(\genfrac{}{}{0pt}{}{M}{N}\right)}$$
For each gene-lncRNA pair, M was the total number of MG risk miRNAs, K was the number of miRNAs interacting with one gene, N was the number of miRNAs interacting with one lncRNA, and x was the number of common shared miRNAs between the gene and the lncRNA. The gene-lncRNA pairs with P-value < 0.01 were considered as statistical significance.
We further evaluated the co-expression association of gene-lncRNA pairs obtained above by performing the PCC. We downloaded the expression data of the genes and lncRNAs from the Genotype-Tissue Expression (GTEx, v8 release) project (https://gtexportal.org/home/datasets) [24]. We adjusted the P-value of co-expression analysis using the false discovery rate (FDR). The gene-lncRNA pairs with PCC > 0.7 and FDR < 0.01 were finally considered significant gene-lncRNA pairs.
Construction of the lncRNA-mediated module-associated ceRNA network
Based on the theory that lncRNA shares common miRNA-binding sites with gene and functions as miRNA’s sponge to regulate gene, LMMAC network was constructed. For a given gene-miRNA-lncRNA interaction, both gene and lncRNA that shared common miRNAs were co-expression for merging into a competing triplet. After integrating all of the gene-miRNA-lncRNA competing triplets, we constructed the LMMAC network and visualized using Cytoscape software. In the network, nodes denoted genes, miRNAs, and lncRNAs, and edges denoted their interactions.
Random walk with restart to prioritize lncRNAs for MG
We used random walk with restart (RWR) [25, 26] to prioritize potential lncRNAs in the LMMAC network for MG by simulating a random walker and randomly moving from a set of source nodes to its network neighbors. It is defined as follows:
$${P}^{t+1}=\left(1-r\right){WP}^{t}+{rP}^{0}$$
where P0 is the initial probability vector, which is constructed so that a value of 1 was assigned to nodes representing genes known to be associated with disease while 0 was assigned to other nodes; W is the column-normalized adjacency matrix of the LMMAC network; r is the restart probability at each step of the random walk at source nodes; and Pt is a vector in which the ith element has the probability of being at node i at the time point t.
Furthermore, we analyzed the statistical significance of scores of each candidate lncRNA. Through comparing the scores of lncRNAs in the network following n iterations of that known MG risk genes shuffling, the statistical significance for rejection of the null hypothesis was determined. In order to strictly maintain the network topological properties, random sampling without replacement was performed when doing random disturbance, and the degree distribution was guaranteed the same between selection seed node and the real. In iterations, the times of each lncRNA score that was higher than the actual value was recorded as m. The statistical significance P value for each lncRNA was calculated by the ratio of m and n; in this study, n was set at 1000 times.
Clinical samples
A total of 31 MG patients were recruited for the sampling of peripheral blood from the Second Affiliated Hospital of Harbin Medical University. All the patients met the diagnostic criteria for MG [27]. None of the patients had undergone hormone treatment for the previous 6 months. Meanwhile, 31 sex- and age-matched healthy subjects with no autoimmune disease were enrolled as the control. Peripheral blood samples obtained from participants were put in tubes containing ethylenediaminetetraacetic acid. Peripheral blood mononuclear cells (PBMCs) were isolated using lymphocyte separation medium. Written informed consent was obtained from all individuals and the study was approved by the Ethics Committee of The Second Affiliated Hospital of Harbin Medical University.
Real-time PCR analysis
TRIzol reagent (BioTeke, Beijing, China) was applied to isolate total RNA from PBMCs or Jurkat cells. Total RNA was reverse transcribed into cDNA using a RevertAid First Strand cDNA Synthesis kit (Sangon, Shanghai, China) according to the manufacturer’s instructions. Quantitative real-time PCR analysis was performed using Power SYBR Green PCR Master Mix (Solarbio, Beijing, China). The 2−ΔΔCt method was used for the calculation of relative expression level normalized by GAPDH and U6.
Cell culture
The human embryonic kidney 293T (HEK293T) cells were purchased from Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd (Shanghai, China) and were cultured in Dulbecco’s modified Eagle’s medium (DMEM; Servicebio, China). T cell leukemia line (Jurkat cells) were purchased Procell Life Science & Technology Co., Ltd (Wuhan, China) and were cultured in Roswell Park Memorial Institute 1640 medium (Gibco, USA). All the culture media were supplemented with fetal bovine serum (FBS; Tianhang, Zhejiang, China), penicillin, and streptomycin and incubated at 37 °C in a humidified atmosphere of 5% CO2.
Cell transfection
lncRNA short hairpin RNA (shRNA) of human HCG18, miR-145-5p mimics, miR-145-5p inhibitor, and negative control (NC) were transfected into cells using Lipofectamine 3000 (Invitrogen, USA) following the manufacturer’s protocol. The pool of HCG18 shRNA (shHCG18) contained three shRNAs and three antisense oligonucleotides, which was used to knock down HCG18 expression. The shRNA sequences were as follows: shHCG18-1 target sequence, 5′-AGCTGAAAGTCGACGAAGA-3′; shHCG18-2 target sequence, 5′-CAGCAACTCCTGATGAACA-3′; shHCG18-3 target sequence, 5′-GGTGTAGACAAGACAGCAA-3′; antisense oligonucleotides target sequence-1, 5′-TCTTCGTCGACTTTCAGCT-3′; antisense oligonucleotides targetsequence-2, 5′-TGTTCATCAGGAGTTGCTG-3′; and antisense oligonucleotides target sequence-3, 5′-TTGCTGTCTTGTCTACACC-3′. Total RNA was extracted from transfected cells and subjected to RT-PCR analysis to evaluate transfection efficiency.
Western blot analysis
Following total protein extraction from Jurkat cells using radioimmunoprecipitation assay solution, protein concentration was determined by bicinchoninic acid protein assay kit (Wanleibio, China). Protein separation was performed using sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE), and the separated proteins then were transferred onto polyvinylidene fluoride (PVDF) membranes (Millipore, Billerica, MA, USA). PBS containing 5% non-fat milk was used for membrane blocking at room temperature for 1 h. Subsequently, the membrane was incubated with the following primary antibodies overnight at 4 °C: CD28 and β-actin antibodies (Wanleibio), followed by incubation with HRP Goat Anti-Rabbit secondary antibody (IgG, Wanleibio) at room temperature for another 45 min. Signal development was performed using enhanced chemiluminescence (ECL) (Wanleibio) and data was quantified using Gel-Pro-Analyzer software. β-actin was used as the internal control. The experiments were independently repeated in triplicate.
Luciferase reporter assay
The HCG18 fragment containing miR-145-4p binding sites were synthesized to generate wild type (HCG18-WT) or mutant type (HCG18-MUT). The fragments of HCG18-WT or HCG18-MUT were cloned into the Renilla luciferase gene pmirGLO Luciferase Reporter Vectors (Promega, USA) to generate pmirGLO-HCG18-WT and pmirGLO-HCG18-MUT vectors, respectively. Then, the above vectors were co-transfected together with miR-145-5p mimics or negative control into HEK293T cells using Lipofectamine 3000 (Invitrogen). Luciferase activity was evaluated by a dual-luciferase reporter assay system (Promega, Madison, WI, USA) after 48 h of transfection according to the manufacturer’s protocol.
Cell proliferation assay
The proliferation of Jurkat cells was determined following the protocol of Cell Counting Kit-8 assay (CCK-8) assay (Wanleibio, China). Each well of a 96-well plate was filled with 0.1 ml suspension containing 4 × 103 cells. Cell culture was performed at 37 °C in a humidified atmosphere of 5% CO2. CCK-8 solution was added at 24, 48, 72 and 96 h after the beginning of incubation. The absorbance of each well after incubating at 37 °C for 2 h was measured at 450 nm by using a microplate reader (BIOTEK, USA). The experiments were independently repeated in triplicate.
Apoptosis analysis
The apoptosis rate of cells was estimated with an annexin V-fluorescein isothiocyanate (FITC)/propidium iodide (PI) apoptosis detection kit (Wanleibio, China) in accordance with the manufacturer’s protocol. Jurkat cells were transferred and cultured in a six-well plate for 48 h and, subsequently, digested with 0.25% trypsin (Gibco). Cells were stained with annexin V-FITC and PI for 15 min in dark, followed by flow cytometry using a NovoCyte Flow Cytometer (Agilent Technologies, Inc, Shanghai, China). Data was analyzed using NovoExpress Flow Cytometry Software (Agilent Technologies, Inc, Shanghai, China). The experiments were independently repeated in triplicate.
Statistical analysis
GraphPad Prism version 6.0 was used to perform the statistical analysis. Data were presented as the mean ± standard deviation (SD). Student’s t test was performed to analyze the differences between the two groups, while one-way analysis of variance (ANOVA) to analyze that among multiple groups. Correlation analysis was done using Pearson correlation. A p-value < 0.05 (*) was considered statistically significant, and **indicates p < 0.01, and ***p < 0.001, and ****p < 0.0001.