dataesr/scientific-paragraphs-categorization huggingface dataset:
| Statistical tests and visualization Data analysis was performed in R version 3. 6. 1 (R Core  Team, 2016)  and read counts were normalized using DESeq2  (Love et al. , 2014). Normalized read counts were filtered to contain only repeats with nonzero expression in at least one sample. Binning of repeats by distance to the closest gene was performed by summing the median repeat expression of each caste in 1 kb intervals in pre-and postgenic regions. If a repeat stretched over several bins, the affiliation was determined by the location of the middle of the repeat sequence. Expression of repeats that overlapped by at least 1 base pair with introns were summed in a bin called '0′, corresponding to 0 kb distance to a gene. For the PCA in Figure  2b , normalized expression of the top 500 repeats with the highest variance in intergenic regions (>10 kb distance to the closest gene) was used.
We thus set out to search for eGFR loci with DM/noDMdifference and for novel eGFR loci allowing for difference. For this, we gathered GWAS data on eGFR separately for 178,691 DM and 1,296,113 noDM individuals. We prioritized genes underneath identified loci using in-silico functional evidence. We also evaluated the impact of DM-/noDM-specific weights on the genetic risk score (GRS) for eGFR in data independent from the variant-and weight-identifying step  30.
Results Overview of the GWAS meta-analyses. We first analyzed GWAS data on 7,046,926 single nucleotide polymorphisms (SNPs) and their association with logarithm-transformed eGFR separately for 109,993 DM (type 2 or type 1) and 1,070,999 noDM individuals (stage 1, including 72 studies from CKDGen Consortium  31  and UK Biobank, UKB  32  ; mostly European-ancestry; Supplementary Data 1, Methods, Fig. 1a ). These DM-/noDM-stratified GWAS summary statistics allowed us to apply the difference test, the joint test and the stratified tests (Methods). Based on these tests, we searched for eGFR-associated loci with DM/noDM-difference and for novel eGFR-associated loci allowing for DM/noDMdifference. We sought replication of identified loci in independent DM-/ noDM-stratified data (stage 2, n DM = 68,698, n noDM = 225,114; Million Veterans Program, MVP  33  , Michigan Genomics Initiative, MGI  34  , Trøndelag health study, HUNT  35  ; all European-ancestry; Supplementary Data 1, Methods).
While the discovery+replication design augments confidence in identified loci, this comes at the cost of lower power  36. We also conducted GWAS searches in a combined stage design (stage 1 + 2: n DM = 178,691, n noDM = 1,296,113), to yield significant loci exploiting the full sample size, yet without independent replication. The two stages of data were also used to separate the variant identification and weight quantification for the GRS (stage 1 data) from the GRS association analyses (stage 2 data). The analysis workflow is shown in Fig. 1. Seven eGFR loci identified with DM-/noDM-differential effects. To search for eGFR loci with DM/noDM-difference, we applied two approaches  24, 37  (Fig. 1b ): (i) a genome-wide difference test (P Diff < 5 × 10 -8 , difference test approach), and (ii) a search for genome-wide significant association with overall eGFR followed by a difference test in the same data (P Overall < 5 × 10 -8 and P Diff < 0.
We also queried the Open Targets Genetics database  44  for associations of the 11 variants with other traits. We found 126 genome-wide significant associations (P < 5 × 10 -8 , Supplementary Data 9), particularly for hypertension and blood counts: for 3 of the 11 variants, the eGFR-decreasing alleles were associated with increased risk of hypertension (near UMOD/PDILT, DCDC5, and PIK3CG). The variant rs55722796 near MED1/ NEUROD2, solely associated with eGFR in noDM individuals, was not associated with hypertension, but the eGFR-decreasing allele (in noDM) was associated with decreased blood counts (i. e. , decreased red blood cells, hemoglobin, and hematocrit). GWAS for eGFR in individuals with DM identified 29 genomewide significant loci. Previous searches for SNP-effects on eGFRbased kidney function outcome in DM patients had identified eight genome-wide significant loci (i.
Methods Definition of the outcome eGFR and study-specific participant information. GFR was estimated in all study participants in all studies using serum creatinine measurements via the Chronic Kidney Disease Epidemiology Collaboration (CKD-EPI) formula  20. Study-specific information on the utilized assay and year of measurement is given in Supplementary Data 1. We used the R package nephro  63  , winsorized at 15 and 200 ml min -1 per 1. 73 m 2 , and logarithmized using a natural logarithm. For a better interpretation of effect size, we there used eGFR on the original scale in GRS analyses. Study-specific information on study design, sample size, sex and age, utilized serum creatinine assay, year of measurement, and average eGFR is given in Supplementary Data 1. Each study is conducted according to the declaration of Helsinki; local ethics committees approved research protocols and participants provided written informed consent.
1 (HRC) or 1000 Genomes Project phase 3 v5 (1000Gp3v5) ALL or phase 1 v3 (1000Gp1v3) ALL panel. Each study conducted linear regression GWAS on log(eGFR) using natural logarithm, an additive genotype model as well as adjusted for sex, age, and other study-specific covariates. Details on study-specific genotyping, imputation, and GWAS were described previously for CKDGen  20. For UKB, we utilized the fastGWA tool  65  to conduct GWAS for log(eGFR) based on linear mixed models while accounting for sex, age, age × sex, age 2 , age 2 × sex, and 20 genetic principal components and assuming an additive genetic model, which allowed to include related individuals in the GWAS  65. For stage 2 studies, MVP, MGI, and HUNT, the same analysis plan was distributed. Details on genotyping, imputation, and GWAS in UKB as well as for stage 2 studies are shown in Supplementary Data 1.
For quality control, we excluded variants with low imputation quality, Info < 0. 6, or rare variants with minor allele frequency, MAF < 0. 1%. We utilized the software packages GWAtoolbox  66  and EasyQC  67  for the quality control of study-specific GWAS results. We conducted a correction for genomic control lambda of the results stratified by diabetes status. DM/noDM-stratified GWAS meta-analyses. In stage 1 of our analysis, separately in DM (n DM = 88,829) and noDM (n noDM = 620,665) strata, we conducted fixed-effect inverse-variance weighted meta-analyses of 72 GWAS of log(eGFR) using metal  68  , and then meta-analyzed these results with DM/noDM-stratified GWAS from UKB (n DM = 21,040; n noDM = 414,628; European only). To adjust for population stratification within studies, we applied genomic control (GC) correction  69  to each study prior to the meta-analysis.
To search genome-wide for eGFR loci with DM-/noDM-differential effects, we applied two approaches  24, 37  : (i) a genome-wide difference test (P Diff < 5 × 10 -8 , difference test approach), and (ii) a search for genome-wide significant association with overall eGFR followed by a difference test in the same data (P Overall < 5 × 10 -8 and P Diff < 0. 05/k, k = number of followed SNPs, overall+difference test approach). The two approaches complement each other in terms of power to detect difference loci that depend on the magnitude and DM-/noDM specificity  24, 37. We implemented the approaches in two designs: first, we applied a discovery +replication design, where we searched for significant differences in the stage 1 meta-analysis and moved the selected SNPs to a replication stage using the stage 2 meta-analysis (applying a Bonferroni-corrected alpha-level accounting for the SNPs tested).
Second, to make full use of stage 1 and stage 2 data, we searched for differences in stage 1 + 2 meta-analysis combined (combined stage design). We utilized the R package EasyStrata 72 to apply the difference test approaches to the DM/noDM-stratified meta-analysis results. Approaches to search for eGFR loci allowing for DM/noDM-differences. By allowing for differences between DM and noDM in SNP-effects on eGFR, one can possibly detect novel loci that have been masked in analyses on overall eGFR (without DM-status stratification). Again, we used the meta-analyzed SNP-specific summary statistics to apply two approaches to search for novel eGFR loci allowing for difference: (i) a screen using a two degrees of freedom joint test (P joint < 5 × 10 -8 , joint-test approach)  38  derived from a χ 2 -test using 38 C joint ¼ βDM se DM ! 2 þ βnoDM se noDM ! 2 ð2Þ where βDM and βnoDM are the genetic effect estimates for eGFR in DM or noDM from the stratified GWAS meta-analysis, respectively, and corresponding standard errors se DM and se noDM ; and (ii) a on eGFR association in individuals with or without DM separately (stratified tests approach, P DM < 5 × 10 -8 or P noDM < 5 × 10 -8 ).
Analogously to the search for difference loci, we applied twostage designs, the discovery+replication, and the combined stage designs. We utilized the R package EasyStrata  72  to apply the joint and stratified test approaches to the DM/noDM-stratified meta-analysis results. Variant selection and region definition. In order to derive non-overlapping locus regions and locus lead variants, we clumped genome-wide significant variants from the respective test results (P Diff < 5 × 10 -8 , P Joint < 5 × 10 -8 , P DM < 5 × 10 -8 , or P noDM < 5 × 10 -8 ) as done previously  16  : the most significant variant was selected genome-wide (first locus lead variant) and the corresponding locus was defined as the smallest physical interval on the corresponding chromosome containing this variant such that there were no genome-wide significant variants within 500 kb outside the two borders.
To prioritize genes, we used the results from the Gene PrioritiSation (GPS) published previously  16  for known eGFR loci and generated the GPS de novo for novel eGFR loci accordingly: we queried each gene underneath identified loci and prioritized the gene when it mapped to a 99% credible set variant that was protein-relevant with high predicted deleteriousness (CADD  45  PHRED-Score ≥ 15) or expression-relevant in kidney tissue (eQTL, NEPTUNE  46  , or GTEx v7  76  , falsediscovery rate, FDR < 5% for the eQTL) and for splice quantitative trait loci in kidney tissue (sQTL)  76  (FDR < 5% for the sQTL). We also prioritized genes that were known as human kidney disease monogene: for this, we queried each gene at identified loci for a documented kidney phenotype in human as done previously  16  (Online Mendelian Inheritance in Man® database, OMIM  77  or Groopman et al.
78  ) with additional manual curation by expert review to focus on kidney disease monogenes. GRS analyses. GRS analyses based on the previously identified 634 eGFR variants  16  and based on the variants identified for DM-/noDM-difference on eGFR were conducted in unrelated European-ancestry individuals from HUNT (n DM = 3023, n noDM = 25,484, stage 2 study). For each individual, the GRS was computed as the weighted sum of eGFR-lowering alleles across the 634 variants applying three different weighting schemes: (i) weighted by the respective overall per-variant effect as derived previously (i. e. , estimated in CKDGen and UK Biobank, our stage 1 data 16 ) (ii) weighted by the variant's DM-/noDM-specific effects as derived here in stage 1 depending on whether the individual in HUNT had DM or noDM, respectively (and vice versa), (iii) weighted by the variant's DM-/noDMspecific effect from stage 1 for HUNT individuals with DM/noDM, respectively, when the variant was among the seven with an identified significant difference, and weighted by the overall-effect size otherwise.
Acknowledgements The  Deutsche Forschungsgemeinschaft (DFG, German Research Foundation ) supported the meta-analysis-Project-ID 387509280-SFB1350 (Subproject  C6  to I. M. H. ). A. M. H. , B. R. , and R. T. were supported by  VACSR&D MVP  grant  CX001897. This research is based on data from the  Million Veteran Program, Office of Research and Development, Veterans Health Administration , and was supported by  VACSR&D MVP  grant  CX001897  (A. M. H. ). This publication does not represent the views of the Department of Veteran Affairs or the  United States Government. We conducted this research using the UK Biobank resource under the application number 20272. We thank  Paola Bilani  for collecting author information. Extended acknowledgements are provided in Supplementary Note 4 for all studies, in Supplementary Note 5 for MVP and in Supplementary Note 6 for LifeLines.
Then, by network deconvolution, we determine their coexpressed gene modules. Moreover, we validate the conservancy of their upregulation using the Coronascape database (DB). Finally, tissue-specific regulation of the five predictive hub genes indicates that OAS1 and MX1 expression levels are lower in healthy kidney tissue. Altogether, our results suggest that these genes could play an essential role in developing severe outcomes of COVID-19 in DKD patients.
Materials and Methods For this study, we designed a bioinformatics pipeline that uses transcriptomic data to identify predictive hub genes involved in DKD complications derived from COVID-19 infection as depicted in Figure  1. Used codes are available at https://github. com/kap8416/Transcriptomics-Diabetes-Kidney-Disease. Figure  1. Bioinformatics pipeline for obtaining predictive hub genes involved in the pathophysiological impact of COVID-19 on DKD condition. Highlighted in yellow are the steps to perform the transcriptomics analysis, while those in blue are steps for the interactomics methodology. First, transcriptomic data was acquired from the NCBI GEO database. Datas were pre-processed, and outlier detection with WGCNA analysis was performed. Then, we obtained the DEGs and determined their potential functional roles in biological and molecular pathways.
Those DEGs that overlapped among COVID-19 and DKD pathways were selected, and interactome analysis was performed using protein-protein interaction (PPI) associations. Subsequently, the potential hub genes were identified using the Maximal Clique Centrality (MCC) topological algorithm. Then, their co-expressed gene modules were determined using the "Corto" algorithm. Finally, hub gene expression was validated using analysis of systems levels employing Coronascape COVID-19 datasets and gene expression value comparisons in healthy tissues using the GTEx database.
Differential gene expression analysis Raw data were processed with the "affy" package available in Bioconductor using the function robust multiarray average (RMA) to convert an "AffyBatch" object into an "ExpressionSet" object (https://www. bioconductor. org/). Then, data normalization was performed. After that, we removed the "batch effect" using "ComBat" from the sva, and Weighted Gene Co-Expression Network Analysis (WGCNA) to identify outliers. Also, the collapseRows function available in the WGCNA package for collapsing gene expression of several probes in a single measurement per gene, was used. Finally, differential gene expression analysis was performed using the limma package. Log2FC > 1, and FDR < 0. 05 values were considered to determine the differentially expressed genes (DEGs).
Functional enrichment analysis for DEGs in the DKD dataset The functional enrichment was assessed using the Gene Set Enrichment Analysis (GSEA) of DEGs. Enrichment includes all the Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome pathways. The cluster profiler R package  [40]  was employed for data analysis and SRplot server (http://www. bioinformatics. com. cn/en?keywords=pie) for data visualization. Further, we considered focusing only on 42 genes shared in both DKD and COVID-19 pathogenetic processes.
Functional enrichment analysis for the shared pathogenic DEGs Functional enrichment analysis was performed for the DEGs using KEGG and Reactome pathways, including the GO terms for biological processes. Then, the ClueGO (version 2. 5. 7) module of the Cytoscape software (version 3. 9. 0) was used to examine the inter-relational pathway analysis of the significantly enriched functions to identify the most significant genes  [41]. We identified up and downregulated pathways using the top 200 DEGs. Then, the 42 overlapped genes were considered to perform the same analysis with a kappa score of 0. 4, pvalue < 0. 05. For statistical tests, hypergeometric two-sided and Benjamini-Hochberg methods were employed  [41]. Finally, the Cluster profiler R package was employed for data analysis and visualization  [40].
PPI network and identification of hub genes The list of 42 DE-human genes previously identified in the shared pathogenetic process were analyzed to deep insight into their interatomic roles. Network analysis and visualization were performed using the Cytoscape version 3. 9. 0. STRINGDB platform (https://string-db. org/) was used to obtain physical and functional experimental validated data. For this study, high confidence scores (0. 7) and no additional interactors were filtered, keeping only experimental, co-expression and database characterized interactions. Each node in the PPI network represents a protein whereas edges are interactions and connections between them. Finally, the predictive hub genes were obtained using the CytoHubba module selecting the top five nodes employing the topological algorithm of Maximal Clique Centrality (MCC)  [42].
Deconvolution of the Gene Regulatory Network (GRN) In order to identify co-expression modules of the five hub genes, the "Corto" algorithm with default parameters was used for the GRN inference. This package is freely available on the CRAN repository. The "Corto" algorithm is implemented as a co-expression-based tool that infers GRNs using as an input a TFs list with their targets and an expression matrix data set  [43]. In brief, Corto uses a combination of Spearman correlation and Data Processing Inequality (DPI), adding bootstrapping to evaluate the significant edges by removing indirect interactions. We used the DKD expression matrix and a human-TF list as input. For this analysis, a p-value= 1 × 10-8 was used as a cut-off, and 100 bootstraps were carried out. Subsequently, an output file containing an inferred enriched GRN was obtained.
Functional enrichment analysis for co-expressed DEGs in the GRN The functional enrichment analysis was assessed using the GSEA of each of the five previously identified modules of the hub genes co-expressed DEGs. Functional enrichment considered all the GO terms and the KEGG and Reactome pathways. The cluster profiler R package was employed for data analysis  [40]. For data visualization the parent terms were depicted in the GRN using the Biorender environment (https://biorender. com/).
In silico validation of expression of the predictive hub genes in GTEx database Gene expression data was accessed at the GTEx repository (https://www. gtexportal. org/). GTEx is a resource for the analysis of tissue-specific level datasets. The five hub genes were tested for multi-gene query expression in kidney, lung, heart, pancreas, and whole blood healthy tissues. Data was visualized in a heatmap plot. Then, singlecell data from lung tissue was accessed to determine the specific cell expression of the hub genes. Data was visualized in pie charts.
Identification of Differentially Expressed Genes across Diabetic Kidney Disease First, we performed raw data normalization using the No GSE305029 Affymetrix microarray dataset. This step allowed reducing false positives, as depicted in Figure  2A. Later, principal component analysis (PCA) showed clustered distances among each renal tubule sample, including those belonging to the disease condition (DKD) and those to the control group (Figure  2B ). After this step, one batch was observed for the sample GSM757034 belonging to the control group. This sample was removed for further dataset proccesing. After that, we performed WGCNA. This approach is a widely used data mining and exploratory method. It allowed us to identify clusters and outliers to discard samples that did not cluster with their corresponding group (disease or control).
While For instance, complement activation, defense response, and humoral immune response were the top enriched terms. Moreover, other affected processes identified in the upregulated genes are clustered in angiogenesis and vasculature development events. Besides, these induced genes are participating in cellular response to viruses and to cytokine stimulus. Further, to determine their involvement in multiple molecular-biological processes, an inter-relational analysis was performed using the ClueGO module of the Cytoscape environment to obtain functionally organized pathway-terms networks. ClueGO employs kappa statistics to link the functional terms in the network. As shown in Figure  6A , the upregulated DEGs were mainly involved in the innate immune system, leukocyte differentiation, angiogenesis regulation, cytokine induction, extracellular matrix organization, defense response to other organisms, and regulation of peptidase activity.
Interactome analysis and hub genes identification The PPI network analysis is a useful method for identifying the hub genes playing a critical role in complex disease interactomes. The PPI networks for up-DEGs are shown in Figure  8a. In the predicted PPI networks for the 42 up-DEGs, an operation with a network scoring cutoff of 0. 7 was applied using the MCODE plug-in of the Cytoscape environment. It resulted in seventeen DEGs with the highest interactions in the network. Among those seventeen most relevant genes, we find the following TLR2, TLR7, CXCL10, TRAF2, MX2, JAK1, IGNAR2, and a cluster of HLA family proteins. Then, the five prominent upregulated hub genes, namely STAT1, MX1, ISG15, IRF7, and OAS1, were identified and exhibited protein-protein-validated experimental interactions with a high confidence score of 0.
Overall, our interactome analysis results show that the five hub genes identified are co-expressed in modules with other genes participating in the events associated with heart disease, virus immune response, and genetic regulation at the transcriptional and post-transcriptional levels. 3. 5 In-silico validation of hub genes across COVID-19 and healthy tissue datasets. Finally, we performed in silico validation analysis to determine the expression profiles of the hub genes obtained by the interactome analysis, Figure  9. First, we employed the Coronascape Database, which represents a compendium of large-scale studies of transcriptomics data that allows us to perform biological systems analysis. When comparing the expression of our hub genes in DKD samples against this compendium, we identified that studies that comprise the following immune cell samples: dendritic, Calu3, CD4T, CD16/CD14 monocytes, NK and A549lowMOI cells are the ones that overlapped in the expression-being ISG15 the most conserved different expressed gene (8), followed by STAT1  (7) , OAS1(6), and IRF2 (2), Figure  9a.
Then, we delve into the expression of those genes in healthy tissue linked to the potential comorbid development such as the pancreas, heart, kidney, and lung, Figure  9b. Additionally, we included whole blood comparison of the baseline expression levels of those biomolecules. Our results indicate that MX1, IRF7, ISG15, STAT1 possess lung tissue-specific expression. Single-cell expression testing identified that this specific expressio higher in macrophages and B and T cells, Figure  9c. In contrast, IFR7 is a gene with the highest expression valu whole blood. While OAS1 and MX1 are the genes with the lower expression levels in pancreas, heart, and kid healthy samples, Figure  9b. Hence, the in-silico validation of our hub genes in those expression datasets pointed that ISG15 and MX1 are the most conserved up-regulated genes across COVID-19 datasets of immunological while OAS1 and MX1 are the genes with lower gene expression levels in tissue-specific healthy samples.