RNA-SSNV: A Reliable Somatic Single Nucleotide Variant Identification Framework for Bulk RNA-Seq Data

The usage of expressed somatic mutations may have a unique advantage in identifying active cancer driver mutations. However, accurately calling mutations from RNA-seq data is difficult due to confounding factors such as RNA-editing, reverse transcription, and gap alignment. In the present study, we proposed a framework (named RNA-SSNV, https://github.com/pmglab/RNA-SSNV) to call somatic single nucleotide variants (SSNV) from tumor bulk RNA-seq data. Based on a comprehensive multi-filtering strategy and a machine-learning classification model trained with comprehensively curated features, RNA-SSNV achieved the best precision–recall rate (0.880–0.884) in a testing dataset and robustly retained 0.94 AUC for the precision–recall curve in three validation adult-based TCGA (The Cancer Genome Atlas) datasets. We further showed that the somatic mutations called by RNA-SSNV tended to have a higher functional impact and therapeutic power in known driver genes. Furthermore, VAF (variant allele fraction) analysis revealed that subclonal harboring expressed mutations had evolutional selection advantage and RNA had higher detection power to rescue DNA-omitted mutations. In sum, RNA-SSNV will be a useful approach to accurately call expressed somatic mutations for a more insightful analysis of cancer drive genes and carcinogenic mechanisms.

INTRODUCTION
Cancer is the leading cause of death and an important barrier to increasing life expectancy (Sung et al., 2021). According to GLOBOCAN 2020 estimates of cancer incidence and mortality, 19.3 million new cancer cases and 10.0 million cancer deaths occurred in 2020 (Sung et al., 2021). Somatic mutations are usually induced by environmental factors, and it is well known that their accumulation with aging and evolution in human cells will lead to malignant transformation and eventually cancer (Watson et al., 2013). Thus, comprehensive somatic mutation identification in cancer such as the Catalogue Of Somatic Mutations In Cancer (Tate et al., 2019) (COSMIC database) can help characterize its genomic complexities (Watson et al., 2013) and discover oncogenic mutations and driver genes which significantly influence cancer development . Furthermore, person-level somatic mutations also have their own oncogenic and therapeutic implications in multiple cancers (lung cancer (Skoulidis and Heymach,2019), bladder cancer (Cazier et al., 2014;Wen et al., 2021), and glioblastoma McDonald et al., 2015)), targeting the corresponding mutant proteins or pathways. Currently, most somatic mutation identification studies were based on DNAlevel, actionable practices in somatic mutation detection within whole-genome or whole-exome sequencing data have been developed to facilitate precision oncology (Xiao et al., 2021).
Mutations within exons are supposed to be transcribed into RNA, and be reflected in the translated protein. However, many DNA mutations within exons were not found in RNA because they were located in the non-transcribed allele or had no or low expression (O Brien et al., 2015). Yizhak et al. (2019) reported that 65% of DNA somatic mutations within 243 TCGA tumor samples were not detected in RNA. Rashid et al.(2014) found that only 27% of mutated alleles got expressed in multiple myeloma. The significant lack of DNA mutations in RNA indicated that not all DNA mutations have certain effects finally. RNA can be a reliable source to distinguish mutations that have been expressed to affect cellular functions. Although RNA-seq is mainly used for gene expression and fusion discoveries in clinical oncology , previous studies showed that calling genomic variants in expressed exons using RNA-seq data was feasible and cost-effective (Chepelev et al., 2009;Cirulli et al., 2010;Gonorazky et al., 2019;Piskol et al., 2013;Quinn et al., 2013). The advantages included making the most abundant RNA-seq data resources and discovering rare somatic mutations with the low-level DNA allele fraction at higher sequencing depths in sufficiently expressed genes (Chepelev et al., 2009;Cirulli et al., 2010;Gonorazky et al., 2019;Piskol et al., 2013;Quinn et al., 2013;Liu et al., 2014). However, calling somatic mutations within RNA-seq data was challenging compared with calling variants in WES data. The main challenge was the high false-positive rate, deriving from errors during reverse transcription, misalignment near splicing junctions (exon ends), RNA editing, and modification during post-transcriptional processing (Cirulli et al., 2010;Xu, 2018). Multiple RNA somatic mutation calling tools and pipelines have been developed to remove these false-positive calls, which can be divided into two categories: statistical filtering strategy-based (García-Nieto et al., 2019;Neums et al., 2018;Yizhak et al., 2019) and machine learning-based approaches (Muyas et al., 2020;Sheng et al., 2016). For instance, GLMVC (Sheng et al., 2016) calls RNA somatic mutations based on a bias-reduced generalized linear model trained by the characteristics of RNA-seq data. VaDiR (Neums et al., 2018) integrates results from three variant callers and produced higher precision results through consensus combination but sacrificed sensitivity. RNA-MuTect (Yizhak et al., 2019) comprehensively filtered mutations within artifact sites and achieved optimal performance. RF-RNAMut (Muyas et al., 2020) utilized a machine learning model to distinguish somatic variants from germline variants identified in RNA-seq data. Although existed tools have their advantages and highlights, they had their limitations: (1) unsatisfying precision-recall performance with the maximum reported precision-recall to be 0.87-0.72 (Yizhak et al., 2019), (2) required restricted resources such as DNA and RNA panel of normal (PoN) calls from~6500 GTEx samples to achieve a desired result (Yizhak et al., 2019), and (3) model not specifically trained to recognize excessive artifacts in RNA but to identify germline mutations as negative (Muyas et al., 2020).
Here, we introduce a framework named RNA-SSNV (https://github.com/pmglab/RNA-SSNV). It is a unified framework containing a universal pipeline to call RNA somatic single nucleotide variants from the combination of tumor RNA-seq and normal WES data, a multi-filtering strategy to remove doubtful calls with little loss of sensitivity and a supervised machine learning model to identify somatic mutations and artifacts. Our framework achieved the best overall performance for precision and recall, requiring only public reference resources. To validate the generalization performance of our framework, we utilized RNA-SSNV within TCGA lung squamous cell carcinoma (LUSC), bladder urothelial carcinoma (BLCA), and glioblastoma multiforme (GBM) independent datasets. RNA-SSNV achieved similar performance in the area under curve (AUC) for the precision-recall curve with 0.94 for all three datasets. Given its high precision-recall performance, RNA-SSNV will help exploit expressed somatic variants, further extend the range of RNA-seq applications and make full use of abundant RNA-seq data resources.

Framework Overview
Our RNA somatic single nucleotide variant identification framework (RNA-SSNV) consists of three major components, including a RNA somatic mutation calling step, a multi-step filtering process and a machine-learning based prediction (Figure 1). The underlying hypothesis of RNA-SSNV is that RNA-specific mutations have unique biological and technique features; thus, a comprehensive filtration process and a machine learning model based on these features can substantially improve the accuracy of RNA somatic mutation calling.

Datasets
Our datasets were retrieved from GDC, which had harmonized pipelines (https://docs.gdc.cancer.gov/Data/Introduction/) to generate RNA-seq and DNA-seq data. All RNA-seq datasets were aligned to GRCh38 build using a two-pass method with STAR, which required preprocessing before mutation calling. All DNA-seq datasets were aligned to the GRCh38 reference using bwa (Li and Durbin, 2009) and co-cleaned using the GATK toolkit (McKenna et al., 2010), which can directly be utilized in mutation calling.
We chose the TCGA lung adenocarcinoma (LUAD) cohort as the training dataset that contained the largest patient scale (511) compared with other available cancer cohorts. Our training dataset comprised paired tumor RNA-seq and tumor/normal WES data derived from 511 LUAD patients, which simultaneously generated DNA and RNA somatic mutations. Our independent validation datasets comprised paired tumor RNA-seq and normal WES data derived from 498 LUSC, 441 BLCA, and 198 GBM patients, for which we called RNA somatic mutations to get validating records.

Mutation Calling
Theoretically, calling somatic mutations within RNA-seq data can be easily conducted using callers designed for DNA. Haplotype-based callers (GATK Mutect2 (Benjamin et al., 2019;Cibulskis et al., 2013), TNscope) had been proven to outperform position-based variant callers due to their inherent technical advantage in complex variants and high mutation loading regions (Pei et al., 2020;Xu, 2018). In addition, we queried the TCGA helpdesk and learned that our RNA-seq data (TCGA LUAD, LUSC, GBM, and BLCA projects) were sequenced by the UNC center using poly-T mRNA enriching strategy, which indicated that only transcribed exon regions (GENCODE v22 annotated exon regions) within mature mRNA can be sequenced (Kukurba and Montgomery, 2015) and our paired normal targeted capture exome sequencing (WES) data had a canonical target region (Agilent SureSelect TargetInterval). Thus, we chose to utilize Mutect2 to perform somatic variant calling and only retain mutations within targeted coding regions (overlap of exons and WES targets).
Normally, our STAR-2-pass aligned RNA-seq data required a co-cleaning process to conduct mutation calling. Following GATK recommended procedures (RNAseq Best Practice), our aligned RNA-seq bam was passed to the MarkDuplicates tool to identify duplicate reads and help remove PCR-related artifacts. Next, SplitNCigarReads hard-clipped and reformat some alignments which span introns causing large-scale mistaken indels. Finally, it shall undergo base quality recalibration conducted by BaseRecalibrator and ApplyBQSR to detect and correct patterns of systematic errors in the base quality scores.
After obtaining analysis-ready bam files, we utilized Mutect2 to call RNA somatic mutations from paired tumor RNA-seq and normal WES data, DNA somatic mutations from tumor and normal WES data. For the TCGA LUAD training set, we called RNA and DNA somatic mutations to help construct the training dataset. For TCGA LUSC, GBM, and BLCA validation sets, calling RNA somatic mutations were sufficient to validate our framework's performance. For DNA somatic mutations omitted in RNA which required verification, we applied the force-calling mode in Mutect2 to retrieve their RNA mutational status. Finally, FIGURE 1 | Schematic overview of the framework for RNA somatic mutation identification. RNA calling: RNA-seq and WES data were aligned and co-cleaned accordingly. Mutect2 was used to conduct RNA somatic calling with paired tumor RNA-seq and normal WES data. Features were extracted from outputs of FilterMutectCalls and Funcotator. Multi-filtering: multi-filtering strategy was conducted in Mutect2 called mutations by removing multiallelic, RNA-editing, immunoglobin, and HLA sites. Model prediction: using the trained model, mutations with extracted features were predicted as positive or negative, only positives were regarded as reliable mutations. Result analysis: pairwise analysis can be conducted when DNA evidence was available. RNA-SSNV will output a generic entry

Multi-Filtering Strategy
Before model training or predicting, RNA somatic mutations shall be comprehensively filtered to remove known possible artifacts (García-Nieto et al., 2019;Yizhak et al., 2019). Our multi-filtering strategy included removing multi-allelic mutations, RNA-editing sites, IgG, and HLA regions. For multi-allelic mutations, we removed mutations containing three or more allele types to avoid misaligning artifacts. For RNA editing events, we combined A-to-I RNA editing information from the REDIportal (Mansi et al., 2021) database and further editing information from the DARNED (Kiran et al., 2013) database. We removed all mutations which located in the union set of RNA editing events to prevent these false-positive calls. For IgG regions, we removed mutations falling into IgG genes to avoid noisy alignments (Ye et al., 2013). For HLA regions, we removed the HLA mutations in chromosome 6 which contained a high density of germline variants (Buhler and Sanchez-Mazas, 2011).

Construct a Training Dataset
For all TCGA projects involved in our study, the GDC Data Portal (https://portal.gdc.cancer.gov/) already provided openaccess DNA somatic mutations detected by four different callers MuSE, MuTect2, SomaticSniper, and VarScan (Ellrott et al., 2018) with stringent thresholds. Using the GDC MAF Concatenation Tool (https://github.com/wwysoc2/gdc-maftool), we combined the curated mutations from four callers, and constructed a union set of all available DNA somatic mutations for each cancer type (LUAD, LUSC, BLCA, and GBM) to maximize the sensitivity. In addition, given that GDC somatic variant calling pipeline had strict criteria leading to the loss of some true positive somatic mutations, we called our own DNA somatic mutations using raw sequencing data and retrieved GDC-omitted DNA somatic mutations.
Normally, variations in DNA will be passed and presented in RNA through transcription. Reciprocally, any RNA somatic mutations presented in DNA should be true positive since they have got evidence from DNA. Moreover, other RNA somatic mutations lacking support from DNA will be regarded as true negative. To construct a reliable training dataset for model training, we split our RNA somatic mutations into three categories ( Figure 2) based on evidence from the GDC database and GDC-omitted DNA somatic mutations. Finally, based on the information from FilterMutectCalls output and annotation information of Funcotator, we systematically extracted features for each training record with three categories: variant, genotype, and annotation levels (Supplementary Table S1).

Performance Metrics
Due to the extreme distribution bias for true positive and true negative classes (TP : TN = 1:8), our main purpose was to identify true positive RNA somatic mutations correctly. We chose precision, recall, F1 scores, and areas under the precision-recall curve (PR-AUC) as major performance metrics in our study because they are insensitive to class imbalance. Other metrics derived from the confusion matrix (

Model Training and Validation
Records within the training dataset were split into training and testing subsets (9:1). We utilized the training subset for model parameter tuning, feature selection, and model training. For the testing subset, we utilized them for testing the model's generalization performance.
To handle the imbalanced distribution for TP and TN classes, we chose a weighted random forest classifier (RandomForestClassifier, scikit-learn 0.24.2) to reduce the bias  by assigning inversely proportioned weights to different classes (Zhu and Pierskalla, 2016). First, we utilized recursive feature elimination with 10-fold cross-validation (RFECV, scikit-learn 0.24.2) to select optimal features. Second, we utilized a 10-fold cross-validated grid-search over a parameter grid (GridSearchCV, scikit-learn 0.24.2) to fine-tune optimal parameters (max_depth, min_samples_split, min_samples_leaf, max_features, etc.). Finally, we constructed a machine learning model for RNA somatic mutation identification with optimal features and parameters, and applied it in testing subset to assess its generalization performance. Following the procedures mentioned earlier, we conducted somatic single nucleotide variants calling in LUSC, BLCA, and GBM cohorts, utilized a multi-filtering strategy and built validation datasets based on extracted features. We applied our discriminant model in these validation datasets and retrieved assessing metrics to further demonstrate the generalization performance.
Also, we validated the necessity of introducing a new training dataset from another cancer type. We added the GBM dataset into the initial training dataset and constructed a new randomforest classifier. After retrieving and assessing metrics for the new random-forest classifier within LUSC and BLCA independent validation datasets, we compared them with our initial model's performance.

Model Interpretation and Visualization
We utilized impurity-based feature importance for tree-based machine learning models to help interpret features' contributions within our model. The higher its contribution, the more important the feature. Impurity-based feature importance (Gini importance) is computed as the total reduction of the criterion brought by that feature and retrieved through our model's attribute feature_importances_. Because traditional feature importance mainly focused on overall model interpretation, we also introduced the SHAP (SHapley Additive exPlanations, https://github.com/slundberg/shap) (Lundberg and Lee, 2017) python package to help visualize prediction (Lundberg et al., 2018) and provide local explanations (Lundberg et al., 2020). We provided feature contributions calculated by SHAP for predicted probability and conducted a single prediction's visualization by invoking the force_plot function. We also investigated the feature contributions of the training dataset. We calculated and visualized the sum of SHAP value magnitudes by summary_ plot function in SHAP to show the distribution of each feature's impacts on the model output (lift or lower prediction probability).

Whole Framework Implementation
We built our whole framework using Snakemake (Köster and Rahmann, 2012) and class-oriented python scripts. Snakemake (https://github.com/snakemake/snakemake) was applied to manage standard bioinformatic workflows involved in this study (co-cleaning, calling, and annotation) and conduct task auto-management without complicating shell scripts. Functionoriented python scripts contained feature extraction, model training and testing, and model utilizing function. Both Snakemake-based workflows and python scripts were available within our project repository (https://github.com/pmglab/RNA-SSNV), which helped create reproducible analysis.

Analyze RNA Mutations With DNA Evidence
We integrated predicted RNA somatic mutations with known DNA mutations to analyze the relevance between RNA and DNA. We examined their intersectionality and split them into three parts (RNA-DNA overlap, DNA-only, and RNA-only) and two sub-categories (positive and negative class, Figure 3). Each part and sub-category had its biological implication and interpretation requiring further investigation. The RNA-DNA overlap part stood for RNA mutations with DNA evidence support. DNA-only part stood for DNA mutations not detected in RNA, and we utilized the Mutect2 force-call mode to inspect their coverage status in RNA. RNA-only part stood for RNA mutations not detected in DNA, and most of them were artifacts due to lack of DNA evidence or low sequence qualities.
Cancer driver genes were under positive selection during tumorigenesis (Martinez-Jimenez et al., 2020). Here, we focused on cancer-specific driver genes (https://www.intogen. org/) to explore their enrichment patterns (number distribution, functional impact, and therapeutic power) between expressed (RNA-DNA overlap part) and unexpressed (DNA-only part) somatic mutation panels.  (Frederic et al., 2009), and Sorting Intolerant From Tolerant (SIFT) (Ng, 2003) were top-performing prediction tools on somatic variants Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 865313 5 (Suybeng et al., 2020). Thus, we used the dbNSFP v4.1a  database to annotate missense mutations with the aforementioned prediction scores. The chi-squared test was used to calculate the significance (p-value) of enriched distribution and odds ratio (OR). A two-sided independent t-test was used to determine the significance (p-value) of the difference between the means of two prediction groups.
We also conducted an analysis of transcriptome-wide allelespecific expression (ASE) to identify ASE events in somatic mutations and their impacts on gene expression which affected carcinogenesis. We chose cases containing both tumor and paired-normal RNA-seq data from LUSC and BLCA cohorts (LUSC: 49 cases, BLCA: 19 cases), and curated their gene expression profiles from the UCSC Xena database (https:// xena.ucsc.edu/). Then, we chose only heterozygous SNVs in both tumor RNA-seq and WES data (RNA-DNA overlap part), and implemented chi-squared tests on the RNA and DNA allelic depths with a significance cutoff of p-value 0.01 to identify somatic SNV ASEs (Heap et al., 2010;Liu et al., 2016). Finally, we compared the TPM value of tumor and pairednormal samples of cases harboring the somatic SNV ASEs to examine the alteration of total gene expression, and defined the TPM fold change (FC) of 2 and 1/2 as the thresholds of upregulated and downregulated genes .

General Performance of the Framework
After the initial RNA somatic mutation calling and multi-filtering step, we collected 467,654 mutations in the LUAD training dataset and 721,234, 323,323, and 126,449 mutations in LUSC, BLCA, and GBM independent validation datasets, respectively. To evaluate the effectiveness of multi-filtering strategy, we validated the loss of GDC mutations in the LUAD training dataset ( Figure 4A) and LUSC, BLCA, and GBM independent validation datasets (Supplementary Figure S1). We found that the loss was negligible (0.1%), whereas the reduction of possible artifact calls was rather significant (70%); such preprocessing guaranteed a relatively pure mutation set for training and predicting. Furthermore, our framework's built-in machine learning model was trained and fine-tuned by 10-fold crossvalidation. In total, 37 features from three categories were kept for model training after feature selection conducted in the initial 40 features ( Figure 4B). Finally, our framework achieved 88.0% precision and 88.4% recall rate within the testing dataset ( Figure 4C), and other assessing metrics ( Table 2) were also satisfying. For example, the false-positive rate was 0.014, the falsenegative rate was 0.013, and the true-negative rate was 0.987. Moreover, most RNA somatic mutations were at the upper or lower ends of the bay plan plot according to the predicted probability distribution of the testing dataset ( Figure 4D), which suggested a clear classification result.
For case-level performance, as expected, LUSC and BLCA retained a median precision of 0.885 and 0.876 across cases, but GBM only reached 0.739 median precision ( Figure 4F), contradicting its general precision of 0.902. Such contradiction was caused by four high-mutation-rate (harbored more than 100 DNA mutations) cases having high precision (>0.950). In contrast, most GBM cases had extremely low somatic mutation rates with less than 30 DNA mutations transcribed in RNA. Thus, some less identifiable RNA editing events and novel mutations rescued by RNA can easily twist GBM's caselevel precision but are hard to affect GBM's general precision. In addition, LUSC, BLCA, and GBM reached a median recall of 0.905, 0.880, and 0.857, concordant with their general recall. Also, RNA somatic mutation counts were highly correlated with DNA (Pearson correlation coefficient: LUSC: 0.905, BLCA: 0.937, and GBM: 0.607, Supplementary Figure S2) after excluding outlier cases with extreme mutation counts, suggesting the high accuracy of our framework.
We investigated the contributions of 37 features using an importance plot based on Gini impurity ( Figure 4G) which showed that STRANDQ was the most important feature for discriminating RNA somatic mutations, followed by AF_tumor, TLOD, ROQ, and ECNT with nontrivial feature importance scores. In addition, features containing other sequencing qualities and population allele frequencies also played a role in prediction because they represented mutations' reliability and germline evidence. We found that the prevalent RNA editing allelic changes "A>G" came at the bottom of the importance list, which indicated that our multifiltering strategy adequately removed these editing sites and reduced their influence. Furthermore, we, in detail, illustrated the effects of 37 features on the prediction model by SHAP (Muyas et al., 2020) and ascertained whether their variations lowered or lifted the predicted probability (Supplementary Figure S3). After feature selection, we excluded "A>C," "A>T," and "MMQ_alt" features. Among all allelic change features, "A>G," "C>A," "C>G," and "G>A" were retained. Out of which, "A>G" and "G>A" represented A-to-I (Wang et al., 2021) and C-to-U (Lerner et al., 2019) RNA editing events, and their existence had negative impacts on the model output. On the contrary, "C>A" and "C>G" represented RNA-editing exclusive allelic changes that exhibited positive impacts. Interestingly, we also found that high tumor allele depth for reference base and alternative base had opposite impacts, which indicated that RNA somatic mutations with high reference allele depth or low alternative allele depth in the tumor sample tended to be artifacts.

Evaluation With Known DNA Evidence
We compared RNA-level somatic mutations with DNA-level to investigate the biological mechanisms for their intersection and uniqueness. As a result, we made a tabular overview ( Table 3) and Venn diagrams (Supplementary Figure S4) to illustrate detailed distribution for the combination of RNA and DNA-level somatic mutations. Here, our framework successfully identified authentic mutations from the RNA-only part (which got ignored/not covered in WES data) to increase information gain and improve diagnostic yield. For all three parts, the RNA-only part had the largest mutation counts. The vast majority were labeled as negative (97.7-99.2%), indicating that our framework had successfully identified most artifacts in RNA because all these negative calls shall be filtered in the final output. Interestingly, when comparing mutation counts of the DNA-only part with the RNA-DNA overlap part, we found that less than 1/3 DNA somatic mutations got expressed in RNA. Such phenomenon was concordant with another study, mainly due to insufficient sequence coverage in low-expression or un-expression genes (Yizhak et al., 2019). Further analysis was listed in the following section for elaborate explanations.

Variably RNA-Expressed Mutations Harbored a Special Enrichment Pattern
In detail, we explored the RNA expression ratios (number of expressed DNA somatic mutations/number of all DNA somatic mutations) for each case of three cancer types ( Figure 5A), median expression ratios for LUSC, BLCA, and GBM were 0.312, 0.349, and 0.256, respectively. Highly variable expression ratios (0.000-0.632) in three types of cancer suggested that different DNA somatic mutations had various expression statuses in RNA. Notably, although the brain has a high number of expressed genes than other human tissues (Naumova et al., 2013), expression ratios of GBM were still significantly lower than those of LUSC or BLCA. These results indicated that DNA somatic mutations might be variably expressed or not expressed at all, and RNA somatic mutations were important to evaluate possible expression status.
To investigate whether the RNA-expressed somatic mutations tended to have larger functional impacts than those that only existed in DNA, probably resulting from the positive selection of cancer subclonal, we compared the impact scores of mutations within cancer-specific driver genes (Martinez-Jimenez et al., 2020) between RNA-DNA overlap and DNA-only parts. Interestingly, the cancer driver genes' mutations were enriched in the RNA-DNA overlap part (LUSC: OR = 2.01, p = 1.14 e-68, BLCA: OR = 2.57, p = 9.89 e-119, GBM: OR = 2.70, p = 1.73 e-16. Supplementary Table S2), even though the DNA-only part had excessive mutation counts than the RNA-DNA overlap part (DNA-only/RNA-DNA overlap:~2/1). Moreover, we compared the predicted pathogenicity scores for missense mutations located within cancer driver genes between RNA-DNA overlap and DNA-only parts, and found that all RNA-DNA overlap parts had significantly higher pathogenicity scores across three cancer types and seven prediction tools (p-value < 1 e-5, Figure 5B). The significantly higher prediction scores implied that predicted damaging mutations tended to be selectively expressed in driving tumorigenesis, and our RNA-level somatic mutation identification framework effectively enriched the functional mutations.
Furthermore, we want to explore whether actionable mutations tend to get expressed in RNA and exhibit clinical effects. Thus, we assessed the therapeutic power for mutations in cancer driver genes between RNA-DNA overlap and DNA-only variants using the OncoKB database (Chakravarty et al., 2017) (https://www.oncokb.org/, Supplementary Table S3). Therapeutic sites within the RNA-DNA overlap part were far more than DNA-only across three cancer types (LUSC: OR = 13.34, p = 8.35e-19, BLCA: OR = 3.27, p = 4.26e-16, GBM: OR = 4.26, p = 3. 66e-4, Table 4), indicating that the RNA-level somatic mutations calling can enrich clinical therapeutic variants. Notably, we observed that some therapeutic mutations from the OncoKB database also occurred in the DNA-only part. For example, except for 52 RNA-DNA overlap somatic mutations in BLCA, PIK3CA also had 12 DNA-only somatic mutations with "Level_3B″ OncoKB annotation (Chakravarty et al., 2017). We found that even if the 12 TCGA BLCA cases containing 12 DNA-only somatic mutations had sufficient expression level for the PIK3CA gene (TPM: 23.7-51.7, curated from UCSC Xena (Goldman et al., 2020) dataset), the 12 mutations' alternative allele still got un-expressed (median alt allele-depth: 0) leading to unlikely benefit from certain targeted therapies. Therefore, although PIK3CA is a valuable therapeutic target for inhibitors of PI3K/AKT/mTOR pathways in advanced bladder cancer (Ross et al., 2016;Willis et al., 2020), the detailed expression status of the mutations should be carefully evaluated when the targeted therapy is considered. Such phenomenon was opposed to the assumption that mutations located within sufficiently RNA total-Counts of mutations observed in the total RNA somatic mutation set. P-R-Precision-recall metric for RNA somatic mutations with GDC mutations as a golden standard dataset.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 865313 8 expressed genes had undoubtful effects making them potential therapeutic targets, and RNA-level mutations were required to validate these targets' transcription status.
Given that RNAs were enriched with mutations of higher functional impact and therapeutic value, we assessed the performance of RNA-level somatic mutations for discovering cancer driver genes by other statistical methods. Here, WITER (Jiang et al., 2019) was adopted to test the enrichment of somatic mutations due to positive selection in tumorigenesis (Jiang et al., 2019;Martinez-Jimenez et al., 2020). We compared the significant genes based on RNA-level somatic mutations to those based on the DNA-level somatic mutations in three cancer datasets. Among all significant genes (FDR < 0.1), the RNA somatic mutations led to a higher proportion of known cancer driver genes from the Intogen database (Martinez-Jimenez et al., 2020) in two of the three datasets than DNA (LUSC: 6/7 vs. 6/9 and GBM: 5/5 vs. 5/18, see details in Supplementary Tables  S4, S5) with identical cancer-driver genes, another cancer type (BLCA: 12/18 vs. 15/19) also had a similar proportion. This result suggested that the RNA-level may lead to fewer false-positive estimations for driver genes than DNA-level.
In addition to known cancer driver genes, other significant genes based on the RNA-level somatic mutations, though un-registered in the Intogen database, were also functionally important to cancer development. CTNNB1, for example, had a significant q-value of 0.081 in BLCA. CTNNB1's mutations have been found to cause aberrant WNT/CTNNB1 signaling and are associated with the susceptibility and prognosis of breast, endometrial, and gastric cancers (Kurnit et al., 2017;van Schie and van Amerongen, 2020;Wang et al., 2012). CHEK2 (q = 0.087 in LUSC, q = 0.052 in BLCA) played an important role in the repair of DNA damage, and its heterozygous mutations had been found to be causing genetic susceptibility to lung cancer (Wang et al., 2014) and bladder cancer (Złowocka et al., 2008). Although our detected CHEK2 somatic mutations were not inherited or passed on, their heterozygosity was similar and induced cancer risk. In a word, RNA can also prioritize potential cancer driver genes.

RNA Increased Mutation Detection Power.
VAF (variant allele fraction) was the fraction of sequencing reads harboring the mutation when performing NGS (Friedlaender et al., 2021), measuring the subclonal prevalence of specific mutations (Benard et al., 2021). We compared the DNA VAF distribution for DNA-only and RNA-DNA overlap parts within three cancer types to examine the subclonal selection advantage for expressed mutations. Higher DNA VAF was observed in expressed DNA somatic mutations ( Figure 5C left comparison, p < 1 e-5), indicating the trend of cancer evolution for subclonal harboring RNA somatic mutations. Interestingly, RNA VAF was significantly higher than DNA VAF within expressed mutations of RNA-DNA overlap part ( Figure 5C right comparison, p < 1 e-5), suggesting an expression tendency for the mutant allele. The common cancer WES study has a mutation limit of detection (LoD) at 5% VAF, and reporting these subclonal mutations incurs the risk of sequencing error-induced false positives (Yan et al., 2021). For these low-VAF (<0.05) DNA somatic mutations, their RNA VAFs were much higher, with median values of 0.374 in LUSC, 0.342 in BLCA, and 0.241 in GBM. Therefore, RNA somatic mutations exhibited subclonal selection superiority and increased the power for low-VAF mutation detection.
Here, we, in detail, demonstrated the recovery of DNAomitted mutations for our framework. For the RNA-only part, we found that our framework helped rescue~10% of mutations (Table 3) which were missed based on DNA sequencing data. Most of the rescued mutations had low alternative allele depth (median: 0-1) or alternative allele fraction (median: 0-0.03) in WES data but opposite situations (median alt allele depth: 8-10, median alt allele fraction: 0.31-0.67) in RNA-seq data. There were also 102, 120, and 8 mutations located within cancer driver genes out of 6,997, 6,233, and 969 positive mutations from LUSC, BLCA, and GBM, respectively (Supplementary Table S6). Furthermore, we discovered biologically important cancer variants within these overlooked "driver" mutations using the DoCM database (Ainscough et al., 2016) (http://docm.info). We found that 17 out of 102, 14 out of 120, and 2 out of 8 DNAoverlooked "driver" mutations in LUSC, BLCA, and GBM had literature support from one or more publications (Supplementary Table S7). For example, TCGA-FD-A5BS had TP53 p.R282W mutation rescued by RNA with its reference-alternative allele depth in DNA: 19-1, RNA: 17-14. The R282W mutant had been found to cause the gain of novel oncogenic functions (GOF) in p53 proteins and associate with  Level_1 FDA-approved drug  0  0  43  1  0  0  Level_2 standard care  0  0  0  0  1  0  Level_3 clinical evidence  58  3  140  40  13  1  Level_4 biological evidence  85  10  96  22  28  7  Counts sum  143  13  279  63  42  8  Total  1, poorer cancer outcomes with a more prominent GOF effect (Zhang et al., 2016). Low tumor purity can bias somatic mutation detection with the positive correlation between mutation numbers and tumor purities (Cheng et al., 2020). For example, TCGA-90-6837 in LUSC with its CPE (Aran et al., 2015) (consensus measurement of purity estimations) lower than average (0.56 vs. 0.68) had no official DNA mutation (WES failed to detect), we investigated its RNA somatic mutations identified by our framework to confirm its mutational status. We found that out of its 192 RNA somatic mutations, six mutations fell within cancer driver genes, and their existence had been ignored by WES (Table 5). Among these mutations, KMT2D is a lung tumor suppressor gene (Alam et al., 2020), and its mutation was one of the most significant prognostic factors in LUSC (Ardeshir-Larijani et al., 2018). We found that KMT2D p.E869p mutation could cause its truncation leading to tumor progression. In addition, TP53 p.A276G mutation had been found to locate within the DNA binding domain of the TP53 protein and presumably have deleterious impacts on protein functions (Chang et al., 2021) with pathogenic ClinVar database (Landrum et al., 2020) interpretation (Accession: VCV000185319.3). These findings confirmed that RNA-seq data could provide valuable supplementary information useful for clinical decisions and improve diagnostic yield in extreme cases when DNA failed to detect actionable mutations.

Transcriptome-Wide Allele-Specific Expression Analysis
We calculated the TPM fold change (FC) to measure gene differential expression status. After excluding infinite FC values, we found that the median gene FC for RNA-expressed mutations was significantly higher than unexpressed mutations ( Figure 5D). Thus, genes harboring RNA-expressed somatic mutations tended to have higher expression level in tumor samples than in paired normal samples. We detected somatic SNV-level ASEs, and found that 24.8% of 3876 and 23.2% of 1700 somatic mutations exhibited ASE events in LUSC and BLCA RNA-DNA overlap parts. As expected, most (~90%) ASE somatic mutations had over-expressed mutant alleles. The results showed that certain expressed somatic mutations had higher expression superiority in the mutant allele than the wild allele, which further enhanced the mutation detection power in RNA. Furthermore, we curated gene lists for 10 signaling pathways in cancer (Sanchez-Vega et al., 2018) and explored the functional alteration on signaling pathways for ASE somatic mutations. Ideally, if the ASE somatic mutation is functional, the direction of ASE event for the mutant allele should be the same as the direction of gene expression alteration for tumor vs. paired-normal samples . Thus, we mapped ASE somatic mutations to genes involving cancer signaling pathways with identical expression change direction. Finally, we identified several pathways (cell cycle, HIPPO, RTK RAS, TGF-Beta, and WNT) containing heavily altered genes with ASE events (Supplementary Table S8). Interestingly, seemly "benign" synonymous mutations also contained ASE events and altered gene expression level. For example, NF1 is a tumor suppressor that negatively regulates RAS signaling (Redig et al., 2016). NF1 p.L43L mutation in TCGA-39-5040 had an over-expressing mutant allele (DNA VAF: 0.32, RNA VAF: 0.63) and showed an upregulated gene expression (tumor/paired-normal fold change: 2.53), which activated NF1 function to under-regulate the RAS signaling pathway and suppressed carcinogenesis.

DISCUSSION
Although common somatic mutation detection practices come with WES, important and actionable mutations are often conserved in RNA-seq. Therefore, we developed RNA-SSNV, an integrative framework to identify RNA somatic single nucleotide variants called within tumor RNA-seq and paired-normal WES data. To maximize performance, we combined multi-filtering strategies and a machine-learning model. For the multi-filtering strategy, we found that it removed massive artifacts (~70%) while omitting few true positive calls (~0.1%). Before constructing the classification model, we also evaluated the performance of the GATK-recommended filtering tool (FilterMutectCalls) for the LUAD training dataset and LUSC, BLCA, and GBM validating datasets using precision-recall metrics. The result showed that FilterMutectCalls achieved a satisfying recall but a low precision rate (LUAD P-R: 0.380-0.865, LUSC P-R: 0.399-0.871, BLCA P-R: 0.442-0.886, and GBM P-R: 0.540-0.881), which may lead to large false-positive calls. Because FilterMutectCalls was originally designed based on DNA somatic mutation filtering strategy, which may not be fully compatible with RNA, we adopted a machine learning model with comprehensive features to conduct classification. For model training, we adopted various techniques to ensure its reliability. To construct a high-quality training dataset, we used GDC DNA mutations as the golden standard and self-called DNA mutations as important supplementary When utilized in independent validation datasets (TCGA LUSC, BLCA, and GBM), RNA-SSNV achieved similar performance as in the testing dataset, which had 0.871-0.895, 0.876-0.871, and 0.902-0.830 precision-recall rate, respectively. Not only can our framework reliably detect RNA somatic mutations, but it also can conduct pairwise analysis with provided DNA mutations. Although our framework achieved satisfying performance within somatic RNA single-nucleotide variants' identification, limited scenarios in which only RNA somatic mutations can be retrieved such as the GTEx project (Lonsdale et al., 2013) (contained RNA-seq data from~6700 samples across 29 normal tissues). Common RNA-seq practices involving research always included DNA-seq data which generated somatic DNA mutations simultaneously; thus, the investigation for the relationship between DNA-level and RNA-level somatic mutations was essential. Multiple studies have found that combining DNA-level and RNA-level somatic mutation can achieve maximum performance for mutational investigation (Krug et al., 2018;Newman et al., 2021;Wilkerson et al., 2014;Zhang et al., 2020). Thus, we split DNA and RNA somatic mutations into three parts: DNA-RNA overlap part, DNA-only part, and RNA-only part; and each part had positive and negative sub-parts representing our model's classifications. The DNA-RNA overlap part represented orthogonal validated DNA and RNA mutations; its positive sub-part contained reliable cancer somatic mutations with clinical usage, but its negative sub-part contained false-negative calls misclassified by our model. When using SHAP to analyze these falsenegative calls (Supplementary Figure S5), we found that G>A mutant status had significant impacts, which implicated that APOBEC-mediated C-to-U RNA editing events (Lerner et al., 2019) contributed to misclassification and current RNA editing resources were insufficient to filter C>U editing sites. DNA-only part represented DNA mutations omitted in RNA somatic mutation calling, and we found that some DNA mutations' reference allele got selectively expressed while their alternative allele got silenced. To explore how many DNA-only somatic mutations got selectively expressed, we calculated the selective expression ratios (number of mutations with reference allelic depth>10/number of DNA somatic mutations not identified in RNA) for DNA-only parts across three cancer types (Supplementary Figure S6). The median mutation selective expression ratios for LUSC, BLCA, and GBM were 0.134, 0.120, and 0.154, respectively, confirming that DNA somatic mutations within GBM had higher selective expression tendency than LUSC (p = 0.003) and BLCA (p = 5.63 e-6), possibly due to innate upregulation of DNA repair mechanisms (Ferri et al., 2020). We retrieved their information in RNA using Mutect2's force-calling mode and utilized our model to classify them. Most of them were predicted negative as expected, but a small portion (1.9%) was predicted as positive, suggesting that our selected caller (Mutect2) might have a little neglection. We also observed that mutations' reference allele-specific expression within driver genes leads to doubtful translation effects. In addition, most mutations located within collagen-related genes (COL11A1, COL6A3, COL5A2, etc.) were found silenced while these genes got sufficiently expressed in RNA (Supplementary Table S9). Interestingly, the proteome database (Human Cancer Proteome Variation Database) also contains nearly no evidence for mutant collagen proteins across three cancer types which were abnormal because massive DNA somatic mutations had been found in these genes. The RNA-only part represented RNA mutations without DNA evidence support. Its negative sub-part was artifacts, but its positive sub-part included RNA-rescued mutations missing in DNA that contained mutations within cancer driver genes (1.4%) to provide more therapeutic targets and help with clinical decisions. A major shortcoming of WES is uneven coverage of sequence reads over the exome targets contributing to many low-coverage regions Xiao et al., 2021), and substantial inter-individual variation in coverage of medically implicated genes caused false-negative mutation calls due to low coverage (Barbitoff et al., 2020;Kong et al., 2018). Although using replicate exome-sequencing can improve WES coverage by 4.3-12.7% (Cherukuri et al., 2015), improve variant calling accuracy (Zhang et al., 2014), and enhance clinical interpretation, information redundancy and excess costs limited its usage. Compared with replicate exomesequencing, RNA-seq has improved somatic single nucleotide variants, and clinically actionable mutations are often conserved in RNA.
We also examined the potential of improving our model's performance by introducing additional training data from different cancer types. After adding GBM cancer-type data into the training dataset, we only observed a slight improvement within the testing dataset (recall rate increased 1.3%) and the AUC for P-R curves for TCGA LUSC, BLCA-independent validation datasets remained stable at 0.94 (Supplementary Figure S7). The unchanged performance suggested that our model trained with LUAD datasets probably has already contained key features of RNA somatic mutation in cancer cells and is applicable for other cancers. Although the general performance for our model was identical across three validation datasets, performances under default threshold (0.5) slightly differed and a dynamic shift of threshold according to different aims (prefer higher precision or recall) was required. In addition, due to insufficient C-to-U RNA editing database resources, the current model sacrificed high recall to ensure removing editing events for the G>A mutation type. The high distribution of G>A mutations (52.3%) in false-negative sets of TCGA LUSC-independent validation dataset reflected this imperfection. Therefore, we recommended that users manually review predicted-negative G>A mutations within known driver genes to improve diagnosis. To facilitate user to inspect predictions, we provided codes to visualize the contribution of important features using SHAP library and a canonical table to exhibit all useful information for user-specified records. A major limitation of our framework was that it was designed to identify RNA somatic mutations only from tumor RNA-seq and paired-normal WES data. Future works will include extending RNA somatic mutation identification scope into other sequencing data types (single-cell RNAseq or whole-genome DNA-seq).
For cancer research involving both WES and RNA-seq data, the conventional analysis strategy uses WES data to call somatic mutations and then validates whether somatic mutations exist in RNA-seq data. However, the conventional strategy may still omit some somatic mutations in RNA-seq data. Our study significantly improved the capability to call RNA somatic mutations and further revealed the association between somatic mutations derived from RNA and DNA, providing valuable supplementary information for conventional cancer somatic mutation analysis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.


s2fieldsofstudy: Biology Medicine

text:
High Expression of DEPDC1B Predicts Poor Prognosis in Lung Adenocarcinoma

Introduction Lung adenocarcinoma (LUAD) is the most common type of lung cancer. DEP domain-containing 1 B (DEPDC1B) is involved in the development of several cancers; however, its role in LUAD is unknown. Therefore, we aimed to determine the biological function and prognostic value of DEPDC1B in LUAD. Material and Methods We analyzed the correlation between DEPDC1B expression and the clinical features of LUAD and lung squamous cell carcinoma (LUSC). Survival was evaluated by generating Kaplan–Meier curves, which were used to analyze the relationship between DEPDC1B expression and prognosis in LUAD and LUSC. DEPDC1B expression in tumor and normal tissues from patients with LUAD and LUSC was determined using immunohistochemistry, and its clinical significance was analyzed. Finally, the correlation between the expression and biological function of DEPDC1B in LUAD was examined. Results Our findings revealed that DEPDC1B expression was higher in tumor tissues than that in normal tissues from patients with LUAD and LUSC (P < 0.001). These results were confirmed in clinical samples from patients using immunohistochemistry. Analysis of a dataset from The Cancer Genome Atlas (TCGA) showed that high DEPDC1B expression was associated with poor prognosis only in patients with LUAD (P < 0.001). Similarly, high DEPDC1B expression was related to shorter overall survival (OS) and progression-free interval (PFI) in patients with LUAD. These associations were not observed in LUSC. Functional enrichment analysis suggested that DEPDC1B promoted tumor development in LUAD by regulating the cell cycle. Conclusion High DEPDC1B expression predicts poor prognosis in patients with LUAD. Thus, DEPDC1B has potential as a therapeutic target for LUAD.

Introduction
Lung cancer is a leading cause of cancer and cancer-related death worldwide. [1][2][3] The majority of lung cancers (about 85%) are non-small cell lung cancer (NSCLC), which includes lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC). Adenocarcinoma is the most common pathological subtype. Prior to the 1990s, LUSC was the most common pathological subtype, especially in males. However, the incidence of adenocarcinoma has increased and now exceeds that of squamous cell carcinoma. 4 Lung cancer is one of the deadliest cancers, accounting for approximately 28% of cancer-related deaths. Most lung cancers are diagnosed at advanced stages, and the five-year survival rate is poor. 5 Despite recent advances in the timely diagnosis and treatment of NSCLC, the prognosis of advanced NSCLC remains poor. The diagnosis and treatment of NSCLC are multidisciplinary and require comprehensive clinical, radiological, molecular biological, and pathological data to accurately assess the prognosis and identify therapeutic targets. Studies have shown that molecular biological detection plays a significant role in improving the early diagnosis rate and prognosis of NSCLC. [6][7][8] Therefore, identifying oncogenes associated with prognosis in NSCLC is important for developing molecular biomarkers.
DEP domain-containing 1 B (DEPDC1B) is involved in various signaling pathways. DEPDC1B, located on chromosome 5q12, encodes a 61 kDa protein of 529 amino acids. 9,10 DEP domains play significant roles in several processes. DEPDC1B recognizes G-protein-coupled receptors and regulates signaling pathways through effector and regulatory factors. 11 DEPDC1B interacts with diverse signaling molecules, including splicing regulatory molecules and transmembrane proteins. Studies have shown that DEPDC1B participates in cell adhesion, cell proliferation, and cell cycle regulation. 12,13 DEPDC1B is overexpressed in various cancers and is a potential biomarker and therapeutic target. 9,[14][15][16][17][18] However, few studies have explored the role of DEPDC1B in NSCLC. Therefore, the objective of the present study was to explore the correlation between DEPDC1B expression and the clinicopathological features and prognosis of NSCLC. We also aimed to clarify the biological processes involving DEPDC1B in LUAD.
We analyzed and compared DEPDC1B expression between tumor and normal tissues in LUAD and LUSC in dataset from The Cancer Genome Atlas (TCGA) database. We then assessed DEPDC1B expression in clinical LUAD and LUSC samples using immunohistochemistry. We separately validated the association between DEPDC1B expression and prognosis in patients with LUAD and LUSC. Then, we used gene set enrichment analysis (GSEA) to determine the possible functions of DEPDC1B in LUAD. Our results strongly support DEPDC1B as a biomarker for predicting prognosis and treatment outcomes in LUAD. GSEA and immune-associated infiltration analysis revealed the biological impact of DEPDC1B in LUAD and the underlying mechanism. Our findings indicate that DEPDC1B may be an useful biomarker for predicting the prognosis of LUAD, but not LUSC.

Bioinformatics Data Source
We collected 594 LUAD cases (535 tumors and 59 normal tissues) and 551 LUSC cases (502 tumors and 49 normal tissues) from TCGA (https://cancergenome.nih.gov; Table S1). Next, we divided the cases into two groups, low and high DEPDC1B expression groups, using the median expression level as a cutoff.

Pathological Sample Collection
We collected tumor tissue samples and paired adjacent normal tissues (normal lung tissue more than 5 cm from the tumor margin) from 100 patients diagnosed with LUAD and 60 patients diagnosed with LUSC between December 2017 and December 2020. All patients were admitted to the First Affiliated Hospital of Anhui Medical University and provided written informed consent. The follow-up date was December 31, 2021 or the date of death. We collected the following patient information: age, gender, smoking status, pathological stage, TNM stage, tumor differentiation, follow-up status, progression free interval (PFI), and overall survival (OS). The clinicopathological information for enrolled patients is presented in Table S2. Two pathologists diagnosed lung cancer based on pathological findings and determined histologic types and stages according to the eighth edition of the tumor, node, metastasis (TNM) staging system for lung cancer. 19 This study was approved by the Institutional Review Board of First Affiliated Hospital of Anhui Medical University.

Immunohistochemical Analysis of DEPDC1B Expression
DEPDC1B expression in LUAD and LUSC was evaluated through immunohistochemical analysis of tumor and normal tissue samples collected from 100 patients diagnosed with LUAD and 60 patients diagnosed with LUSC. The tissues were fixed with formalin and embedded in paraffin, and DEPDC1B was detected using a rabbit anti-DEPDC1B polyclonal antibody (ab237542; Abcam). The paraffin embedded slides were deparaffinized in xylene, dehydrated using an ethanol gradient, and rinsed with distilled water and phosphate-buffered saline. Antigens were extracted by incubation in EDTA buffer for 30 min. Sections were incubated with an anti-DEPDC1B antibody (1:100 dilution) overnight at 4°C and then incubated with HRP-conjugated goat anti-rabbit IgG (ab111909; Abcam) at room temperature for 2 h. Finally, the slides were developed with DAB, stained with hematoxylin, dehydrated, mounted, examined, and evaluated. Two independent pathologists examined the prepared slides under a light microscope (400×) along with positive and negative controls. DEPDC1B overexpression was calculated based on the percentage of stained cells and staining intensity using a semi-quantitative scoring system. 20 Scores were assigned according to the percentage of stained cells as follows: 0 (0-5%), 1 (5-25%), 2 (25-50%), 3 (50-75%), and 4 (75-100%), and staining intensity was scored according to the color of the cell as follows: 0 (no staining), 1 (light staining), 2 (moderate staining), or 3 (intense staining). The scores for the percentage of stained cells and staining intensity were multiplied to obtain the final score (0-12 points). We defined 0-2, 3-6, and 7-12 points as no expression, low DEPDC1B expression, and high DEPDC1B expression, respectively.

Differential Gene Expression Analysis
We divided patients with LUAD into low and high DEPDC1B expression groups using the median z-score. Differentially expressed genes (DEGs) were analyzed using R software with a |log fold change| (log FC) >1 and adjusted P value <0.05 as threshold values. 21 Volcano and heat maps were created to visualize the results of the analyses.

Gene Set Enrichment Analysis
The biological function of DEPDC1B in LUAD was explored using Gene Set Enrichment Analysis (GSEA) of Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GSEA has the advantage of gene sets grouped based on common biological function, chromosomal regulatory mechanism, or location. 22 Enrichment results were considered statistically significant when they showed a false discovery rate (FDR) <0.25 and an adjusted P-value <0.05.

Immune Infiltrate Analysis
Immune infiltration in LUAD samples was analyzed by single-sample GSEA, and infiltration of 24 immune cell types was analyzed using the GSVA package in R (http://www.biocondutor.org/package/release/bioc/html/GSVA.html). 23

Statistical Analysis
We analyzed DEPDC1B expression in unpaired and paired samples using the Wilcoxon rank-sum and Wilcoxon ranksum sign tests, respectively. A nomogram of the clinical variables associated with DEPDC1B in LUAD was generated using the RMS R package. 24 The prognostic value of DEPDC1B in LUAD and LUSC was assessed using Kaplan-Meier analysis. Univariate and multivariate Cox regression analyses were employed to determine the prognostic factors in LUAD. All statistical analyses were conducted and plots were generated using R software (version 3.6.3). A P-value less than 0.05 was considered statistically significant.

Elevated Expression of DEPDC1B in LUAD and LUSC
To examine DEPDC1B expression in patients with LUAD and LUSC, we compared the expression levels in tumor and normal tissues in a dataset from TCGA. In LUAD, DEPDC1B levels were significantly higher in tumor tissues than in normal tissues (P < 0.001; Figure 1A) and were higher in tumor tissues than in paired normal tissues (P < 0.001; Figure 1B). Similar results were observed in LUSC (P < 0.001; Figure S1A and B). To verify this difference in DEPDC1B expression, immunohistochemical analysis was performed using paired normal and tumor tissues collected from 100 patients with LUAD and 60 patients with LUSC. The results of our analysis confirmed elevated DEPDC1B levels in LUAD tumor tissue (P < 0.001; Figure 1C), and representative immunohistochemical images are shown in Figure 1D and E. Similar results were observed in LUSC (P < 0.001; Figure S1C-E).

Relationship Between DEPDC1B Expression and the Clinicopathologic Features of LUAD and LUSC Obtained from a TCGA Dataset
Kruskal-Wallis analysis showed that DEPDC1B expression levels in patients with LUAD were significantly related to T stage (P < 0.001), N stage (P = 0.046), primary therapy outcome (P = 0.004), and pathologic stage (P = 0.028; Figure 2A-D). DEPDC1B expression levels in patients with LUSC were significantly related to T stage (P = 0.042), N stage (P < 0.001), and pathologic stage (P< 0.001) ( Figure S2A-C). Results of the Mann-Whitney U-test revealed that DEPDC1B levels were significantly related to gender (P = 0.003), age (P = 0.025), and smoker (P = 0.048; Figure 2E-G) in LUAD, but only significantly related to age (P = 0.003; Figure S2D) in LUSC.

High Expression of DEPDC1B is Associated with Poor Prognosis in LUAD
Kaplan-Meier survival analysis of TCGA dataset showed that patients with high DEPDC1B-expressing tumors had shorter survival than those with low DEPDC1B-expressing tumors (OS, Figure 3A and PFI, Figure 3B; P < 0.001). These results were validated in 100 patients with LUAD (OS, Figure 3C and PFI, Figure 3D; P < 0.001). Interestingly, no significant difference was observed in patients with LUSC in TCGA dataset ( Figure S3A and B) or our 60 LUSC patient cohort ( Figure S3C and D). Therefore, we only conducted subgroup analysis of patients with LUAD in TCGA dataset, which showed that high DEPDC1B expression was correlated with OS for different TNM stages: T1 and T2 (P = 0.002; Figure S4A), T3 and T4 (P = 0.016, Figure  S4B); N0 (P = 0.021; Figure S4C), N1, N2, and N3 (P = 0.065; Figure S4D); M0 (P < 0.001; Figure S4E) and M1 (P = 0.75; Figure S4F). Similar results were observed for clinical stages (I and II, P = 0.002, Figure S4G; III and IV, P = 0.116, Figure  S4H). These results suggest that DEPDC1B expression levels impact the prognosis of patients with LUAD ( Figure 4A).

Univariate and Multivariate Cox Analyses of Prognosis in Patients with LUAD
We examined the correlation between DEPDC1B expression and clinical characteristics of patients with LUAD in TCGA dataset using univariate Cox analysis. Clinicopathological variables associated with shorter survival included TNM stage, pathologic stage, primary therapy outcome, and high DEPDC1B expression. To comprehensively explore the factors related to high DEPDC1B expression and survival, multivariate Cox regression analysis was performed; high DEPDC1B expression remained an independent factor correlated with OS (Hazard ratio [HR]: 1.949, 95% confidence interval [CI]: 1.292-2.940; P = 0.001) and was correlated with N stage and primary therapy outcome (Table 1). We also conducted univariate and multivariate Cox analyses of our 100-patient cohort with LUAD. Univariate analysis showed that gender,  Table 2. Multivariate analysis was performed using variables with P values <0.05 in the univariate analysis, and the results showed that male sex, terminal pathologic stage, poor differentiation, and high DEPDC1B expression were independent risk factors for OS and PFI of patients with LUAD, and the detailed results are shown in Table 2.

Construction of a Nomogram Based on DEPDC1B Expression in LUAD
To better predict the outcomes of patients with LUAD, we constructed a nomogram using the multivariate analysis of TCGA dataset. This line chart incorporates the clinical features independently associated with survival in the multivariate analysis (ie, TNM stage, primary therapy outcome, and DEPDC1B expression; Figure 4B). High DEPDC1B expression was associated with poor prognosis, in agreement with the results of the multivariate analysis. The c-index of the line chart was 0.719 (95% CI: 0.693-0.745). The calibration plot ( Figure S5) showed that the predicted values are in good agreement with the observed values.

Identification of DEGs Between LUAD Patients with Low and High DEPDC1B Expression
We used R (with |logFC| >1, and an adjusted P value <0.05) to analyze TCGA dataset ( Figure 5). We identified 3533 DEGs (2325 upregulated and 1208 downregulated) between the low and high DEPDC1B expression groups. Volcano and heat maps were created to visualize the results ( Figure 5A and B).

DEPDC1B-Related Signaling Pathways and Enrichment Analyses in LUAD
To elucidate the role of DEPDC1B in LUAD, we conducted GO enrichment analysis of the 1400 DEGs most closely related to DEPDC1B using the ClusterProfile R package. We detected 253 enriched terms in GO biological process categories, including "organelle fission", "nuclear division", and "chromosome segregation" (Figure 5C), which indicates that DEPDC1B expression is involved in diverse biological processes. Sixty-four GO terms were related to the cellular component category "chromosomal region" ( Figure 5D). Furthermore, the analysis suggested enrichment in the molecular function category "receptor-ligand activities" ( Figure 5E).
Related signaling pathways were identified in the low and high DEPDC1B expression groups using GSEA, and results with an adjusted P value <0.05 and a FDR <0.25 were regarded as significant. The following eight pathways were significantly different between the two groups: retinoblastoma gene in cancer, cell cycle checkpoints, PLK1 pathway, mitotic spindle checkpoint, deposition of new CENP-A-containing nucleosomes at the centromere, meiotic recombination, DNA replication, and resolution of sister chromatid cohesion ( Figure 6A-H).

Correlations Between Immune Cell Infiltration and DEPDC1B Expression
We next explored the correlation between DEPDC1B expression levels and infiltration of 24 immune cell types in LUAD ( Figure 7A). DEPDC1B expression was negatively associated with mast cells, eosinophils, immature dendritic cells (iDCs), follicular helper T cells (TFHs), dendritic cells (DCs), CD8 T cells, and natural killer cells ( Figure 7B-E) and was positively associated with type-2 T helper (Th2) cells and gamma delta T (Tgd) cells ( Figure 7F and G).

Discussion
Lung cancer is one of the deadliest cancers worldwide, and LUAD is the most common histological subtype. 2 Developments in chemotherapy and advances in targeted therapy and immunotherapy have led considerable progress in the treatment of LUAD. However, the 5-year survival rate of LUAD remains unsatisfactory. 25,26 Owing to the insidious symptoms of early lung cancers, many patients are diagnosed at an advanced stage. Therefore, identifying molecular biomarkers of LUAD is critical to develop novel therapeutic targets and establish optical treatment strategies. The DEPDC protein family comprises seven members, which are characterized by the presence of a DEP domain. Accumulating evidence indicates that the DEPDC protein family is involved in the development and progression of several cancers, especially hepatocellular carcinoma, and is closely related to poor prognosis. 27 DEPDC1B, a member of the DEPDC1 family, is involved in the growth and progression of numerous tumors. 9,[15][16][17][27][28][29][30][31] We performed bioinformatics analyses using RNA-seq data from TCGA to assess the prognostic value of DEPDC1B for LUAD and LUSC. The results indicated that high DEPDC1B expression was related to poor OS and PFI in patients with LUAD, but this association was not observed in LUSC. This result was also confirmed using immunohistochemistry. Overexpression of DEPDC1B in LUAD was associated with poor clinicopathological factors, indicating that DEPDC1B may function as an oncogene in LUAD. Typically, LUAD patients with high DEPDC1B expression had poorer OS than patients with low DEPDC1B expression, while no significant difference related to DEPDC1B expression was observed in LUSC. This finding differs from a previous study showing that DEPDC1B acts on the Wnt/β-catenin signaling pathway to accelerate NSCLC migration and invasion, suggesting that elevated DEPDC1B expression may be associated with poor prognosis in both types of NSCLC. 32 Previous investigations have indicated that most tumors are related to activation of the Wnt/β-catenin signaling pathway, which is involved in major developmental processes and may promote tumor development by influencing immune exclusion and immunosurveillance. 33,34 Therefore, the Wnt/β- catenin signaling may also play an important role in NSCLC. Although DEPDC1B was previously shown to have an oncogenic effect on the progression of NSCLC, the sample size in the study was relatively small, and the possible biological functions of DEPDC1B in LUAD and LUSC have not been comprehensively explored. The results of our study indicate that DEPDC1B plays a role in LUAD but has no significant impact on the prognosis of LUSC. Accordingly, DEPDC1B is more valuable for the prognosis of LUAD.
Given that DEPDC1B expression did not significantly influence the prognosis of patients with LUSC in the present study, functional enrichment and immune infiltration analyses were performed for LUAD only. Targeted therapies for LUAD and LUSC differ substantially; 35 targeted therapies for LUAD tend to be more sensitive than those for LUSC. Enrichment analysis suggests that DEPDC1B expression is related to various biological processes, such as mitosis, and plays a significant role in the cell cycle. Studies have also shown that DEPDC1B acts at G2 phase of the cell cycle to control the progression of mitosis through adhesion-dependent signaling mechanisms. 12 Sustained proliferative signaling is one of the most important features of tumors. 36 Cancer cells are unstable, and our results suggest that high DEPDC1B 4178 expression may accelerate mitosis and promote cell migration and proliferation, resulting in tumor development, consistent with previous reports on other cancers. 16,31 Evasion of immune destruction is a well-known characteristic of tumors, and immune cell infiltration is an important feature of the tumor microenvironment. The composition and distribution of infiltrating immune cells can impact prognosis. We investigated the relationship between DEPDC1B expression and immune cell infiltration in patients with LUAD. Our results showed that DEPDC1B expression was positively correlated with the levels of Th2 and Tgd cell infiltration. Studies have shown that T-cell stalling dominates the lung cancer landscape, along with immune dysregulation, such as Th1/Th2 cytokine ratio imbalance and T-cell infiltration, leading to tumor invasion. [37][38][39] A previous report showed that a Th2 response and tumor immune evasion lead to poor prognosis; however, this study was conducted in mice, and no significant progress has been made in human studies. Conversely, the development of lung cancer was associated with the expansion and phenotypic changes of resident Tgd cells in the tumor microenvironment. 40 In different cancers, Tgd cells regulate  tumor activity through distinct mechanisms. 41,42 Our findings revealed that DEPDC1B expression in LUAD was negatively associated with infiltration of mast cells and eosinophils. However, infiltration of eosinophils and mast cells has been correlated with a more favorable prognosis in several malignancies. [43][44][45] These results suggest that DEPDC1B might influence immune infiltration to alter the prognosis of patients with LUAD. Our study has several limitations. First, the study included only one dataset and 160 clinical samples. Thus, the sample size was relatively small, and studies of larger samples are needed to verify the results. Second, although the relationship between DEPDC1B expression and immune cell infiltration was examined using bioinformatics, the specific mechanism requires experimental verification.

Conclusion
Our results showed that the DEPDC1B expression was significantly higher in tumor tissues than in normal tissues in both LUAD and LUSC. We examined the prognostic value of DEPDC1B in LUAD and LUSC, and the results suggest that high DEPDC1B expression predicts poor prognosis in patients with LUAD but not in patients with LUSC. DEPDC1B has potential as a biomarker for predicting the occurrence, progression, and prognosis of LUAD. However, additional experiments are needed to clarify the biological function of DEPDC1B and the underlying mechanism.

Ethics Approval and Consent to Participate
This research had been approved by the Institutional Review Boards of First Affiliated Hospital of Anhui Medical University. This project was following the principles of the Helsinki Declaration.

Consent for Publication
All patients signed an informed consent approved by the Institutional Review Board.


Coxsackievirus A6 strains causing an outbreak of hand-foot-and-mouth disease in Northeastern Brazil in 2018

ABSTRACT Hand-foot-and-mouth disease (HFMD) is a highly contagious viral disease commonly associated to Enteroviruses (EV). During 2018, Brazil faced massive HFMD outbreaks spread across the country. This study aimed to characterize the EV responsible for the HFMD outbreak that occurred in Paraiba State, Brazilian Northeastern region, in 2018, followed by a phylogenetic analysis to detail information on its genetic diversity. A total of 49 serum samples (one from each patient) collected from children ≤ 15 years old, clinically diagnosed with HFMD were tested for EV using conventional RT-PCR and RT-qPCR. EV infection was confirmed in 71.4% (35/49) of samples. The mean and median ages were 1.83 years and one year old, respectively. Twenty-two EV-positive samples were successfully sequenced and classified as EV-A species; 13 samples were also identified with the CV-A6 genotype. The phylogenetic analysis (VP1 region) of three samples revealed that the detected CV-A6 strains belonged to sub-lineage D3. The CV-A6 strains detected here clustered with strains from South America, Europe and West Asia strains that were also involved in HFMD cases during the 2017-2018 seasons, in addition to the previously detected Brazilian CV-A6 strains from 2012 to 2017, suggesting a global co-circulation of a set of different CV-A6 strains introduced in the country at different times. The growing circulation of the emerging CV-A6 associated with HFMD, together with the detection of more severe cases worldwide, suggests the need for a more intense surveillance system of HFMD in Brazil. In addition, this investigation was performed exclusively on serum samples, and the analysis of whole blood samples should be considered and could have shown advantages when employed in the diagnosis of enteroviral HFMD outbreaks.

INTRODUCTION
Hand-foot-and-mouth disease (HFMD) is a highly contagious acute viral disease that commonly affects children < 5 years old. HFMD is characterized by a mild febrile illness with multiple oral ulcers, and vesiculo-papular eruptions on the hands, feet, mouth and buttocks, which normally undergo spontaneous involution without major complications 1 .
The most commonly associated viruses leading to HFMD infections are enteroviruses (EV), non-enveloped single-stranded RNA viruses of the genus Enterovirus (family Picornaviridae; realm Riboviria) 1 . EV genome is composed of one open reading frame encoding four structural viral proteins (VP1 to VP4 ) and seven nonstructural proteins (2A-2C and 3A-3D). VP1 is the main coding region for neutralizing epitopes containing virulence determinants, being used for serotyping and molecular epidemiology studies. Seven EVs species have been associated with human disease: EV-A to -D and rhinovirus A to C 1,2 . EV-A is responsible for more than 90% of HFMD cases worldwide 3,4 . Enterovirus A 71 (EV-A71), coxsackievirus A6 (CVA6) and coxsackievirus A16 (CVA16) are the leading pathogens associated with HFMD outbreaks globally [2][3][4] . Although EVs have long been recognized as the most common etiological agents of HFMD1 1 , they have been poorly investigated in Brazil [5][6][7] . HFMD is not a notifiable disease in Brazil, and as a consequence, associated epidemiological data are scarce probably due to the lack of a continuous and systematic laboratory surveillance 6 . On the other hand, HFMD outbreaks must be notified 8 . During 2018, Brazil faced HFMD outbreaks spread across the country 9,10 , but the EV type (or types) involved in those outbreaks are still unknown.
The aim of the present study was to identify the EV genotype implicated in Brazilian HFMD outbreaks that occurred in Paraiba State, Northeastern region, during 2018. The phylogenetic analysis of VP1 sequences was also performed to obtain more detailed information on the genetic relationships with previously reported strains and to identify the lineage implicated in the outbreaks of 2018. In addition, this investigation was performed exclusively on serum samples, and the analysis of whole blood samples should be considered and could have shown advantages when employed in the diagnosis of enteroviral HFMD outbreaks.

Study population
From February 26 th to March 11 th 2018, a total of 49 children (≤ 15 years old) sought the outpatient clinic of Arlinda Marques Children's Hospital with suspected HFMD, were included in this study. Arlinda Marques Children's Hospital is the largest pediatric hospital of the Joao Pessoa city, the capital of Paraiba State. Paraiba is the third most densely populated State of the Brazilian Northeastern region. It is a tourist and industrial hub, and the city of Joao Pessoa is the sea-bordered State capital, with an estimated population of 817,511 inhabitants 11 . The definitive HFMD diagnosis was performed at the Arlinda Marques Children's Hospital based on clinical and epidemiological criteria ( Table 1). The HFMD cases were defined based on the presence of oral ulcers, mainly on the tongue and oral mucosa, as well as on the hard and the soft palate, accompanied by typical vesicular eruptions most commonly observed on the extensor surfaces of the hands, feet, buttocks and/or knees 1 . Basic demographic data were collected including age and gender, the date of the onset of disease, the date of blood sampling and clinical manifestations (cough, diarrhea, vomiting, fever and skin lesions). At the time of clinical onset, all 49 children were considered immunocompetent and were defined as mild HFMD, since no neurological and cardiopulmonary complications such as aseptic meningitis, encephalitis, encephalomyelitis, acute flaccid paralysis, as well as pulmonary hemorrhage or cardiorespiratory failure were observed 4,12 . Some patients presented with respiratory or gastrointestinal symptoms and both findings are recognized as part of EV infections 13 . The HFMD cases were reported to the Epidemiological Surveillance Center of Paraiba State (CIEVS/PB), which initiated an investigation in order to identify the EV genotype associated with the HFMD cases. A total of 49 serum samples (one from each patient) collected during the HFMD acute phase were initially sent to the Public Health Laboratory of Paraiba (LACEN-PB) and stored at -20 °C. Oropharyngeal swabs were not available. Frozen serum specimens were then sent to the Institute of Tropical Medicine, University of Sao Paulo (IMT-USP) in order to confirm the EV infection and identify the involved genotype.

EV detection
Viral RNA was extracted from serum samples using the PureLink RNA Mini kit (Invitrogen, Life Technologies, Carlsbad, CA, USA) according to the manufacturer's instructions. EV was screened by both real-time RT-PCR (RT-qPCR) 14 and conventional semi-nested RT-PCR targeting the 5'UTR region 15 . As a control for PCR inhibitors and to monitor RNA extraction efficiency, each sample was tested by qRT-PCR for the presence of the human ribonuclease (RNase) P gene following standards protocols 16 . Conventional semi-nested RT-PCR assays were monitored for false-positive or false-negative results using reactions controls containing a positive control, a negative control and sterile water.

Amplification of partial VP1
All EV-positive samples,obtained by RT-qPCR and/or conventional semi-nested RT-PCR, were submitted to VP1   amplification. Reverse transcription was carried out with random hexamer primers, and the region encoding the VP1 region of EV was used in a semi-nested PCR following the protocol described by Iturriza-Gómara et al. 17 . The primers designed to EV genogroups A, B and C amplified fragments ranging from 355 to 414 base pairs (bp) 17 .

Sequencing and genotyping
PCR amplicons were sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) with the same primers used in the semi-nested PCRs targeting the 5' UTR and/or the VP1 regions 15,17 . Dye-labeled products were sequenced in an ABI 3130 DNA Analyzer (Applied Biosystems). Sequences were edited with the Sequencher 4.7 (Gene Codes Corporation, Ann Harbor, MI). The Enterovirus Genotyping Tool (National Institute for Public Health and the Environment -RMVI) 18 was used to assign the genotype of the studied EV strains.

Phylogenetic analysis
The phylogenetic analysis was carried out using VP1 partial coding-region sequences (310 bp) in order to infer the lineages and phylogenetic relationships among relevant strains. The VP1 region sequences generated here and a set of prototype EV sequences available in GenBank were aligned using the CLUSTAL W algorithm in the BioEdit Sequence Alignment Editor software, version 7.0.5.2 (Ibis Therapeutics, Carlsbad, CA, USA). The selection of the best model considered the likelihood ratio test (LRT) implemented in the jModeltest software, version 2.1.10 (University of Vigo, Vigo, Spain), and the Generalized time reversible (GTR) model using gamma distribution with invariants sites (G + I) was applied. A maximum likelihood tree was constructed using the FastTree 2.1 (Lawrence Berkeley National Laboratory, Berkeley, CA, USA). To obtain reproducible data and provide a major tree clustering reliability pattern, the statistical support of branches was evaluated by the approximate likelihood ratio test (aLRT).

Accession numbers
The sequences of the present study were deposited in GenBank under the accession numbers: MZ457015, MZ464040-MZ464053 and MZ476967-MZ476973 for the 5´UTR region and MZ779206-MZ779208 for the VP1 region.

RESULTS
EV infections were laboratory confirmed in 71.4% of samples (35/49). RT-qPCR and conventional seminested RT-PCR were positive for EV in serum samples of 22 (44.9%) and 32 (65.3%) patients, respectively. Both methods resulted positive in 19 (38.8%) patients. The day of blood sampling ranged from the first day after the onset of symptoms to the sixth days, most samples were obtained on the first day (42.8%; 21/49). The percentage of confirmed viremia was the highest on the first day of illness (40%; 14/35) and viremia was lower in patients sampled after day 3 (14.3%; 5/35). Twenty-six of the cases presenting with viremia were detected within the first three days of disease (74.3%; 26/35). Table 1 summarizes the clinical findings of HFMD patients. The age of the 49 patients ranged from four months to 15 years old, being 35 (71.4%) infants under the age of one year old, 11 (22.4%) were toddlers between 2-5 years old and three (6.1%) children were ≥ 6 years old. No significant difference in the frequency of EV infection was observed when children were sub-divided into these three age groups. The mean and median ages were 1.83 years and one year old, respectively. A total of 32 (65.3%) samples were from males, and 17 (34.7%) were from females (χ 2 =3.6234, p=0.056). Thirty-three (67.3%) patients showed common clinical signs and symptoms of HFDM, with fever and vesicles, mostly on hands, feet and oral mucosa. Eight cases (16.3 %) exhibited vesicular eruptions on the buttocks and perianal area. Twenty-two patients (44.9%) showed severe skin manifestations all over the body, and one had seizures associated with the severity of manifestations. There were no reported cases of meningitis, encephalitis or fatalities.
A 450 bp fragment of the 5'UTR region was sequenced to classify the detected EVs into species for all 32 conventional RT-PCR-positive samples. Twenty-two of the 32 (68.75%) positive samples could be successfully sequenced. All 22 sequenced samples were genotyped as members of EV-A. Ten of them were genotyped as CV-A6 according to the data obtained by the Enterovirus Genotyping Tool, and in twelve strains the genotype could not be assigned. CV-A6 could also be obtained for three extra samples by direct VP1 sequencing of RNA from serum samples, confirming the genotype previously obtained by the 5'UTR sequencing (Table 1).
To investigate the genetic relationships among contemporaneous CV-A6 strains, the three Brazilian VP1 partial sequences (nucleotides at positions 2588-2888; 310 bp comprising 102 codons) were characterized here (BRA-38B, BRA-41B and BRA-48B strains) and compared to 42 homologous sequences from worldwide countries detected between 1949 and 2019, including nine CV-A6 sequences previously isolated in Brazil. The CV-A6 lineages (A-D) were determined as formerly proposed 19 . The three Brazilian CV-A6 samples detected in the city of Joao Pessoa were classified as belonging to clade D, sub-lineage D3. In addition, CV-A6 strains BRA-38B, BRA-41B and BRA-48B grouped together with strain 52403 (MT74440) detected in Brazil in 2018, but were separated from the other Brazilian CV-A6 strains sub-lineage D3 that had circulated earlier in the country, from 2012 to 2017 6,13 (Figure 1).
The phylogenetic tree has also shown that the Brazilian CV-A6 strains appears to belong to a distinct phylloclade (herein designated D3-VS clade) and this finding was highly supported by the statistical test (hLRT=0.941). Clade D3-VS includes strains that were simultaneous circulating in Latin America countries, including Argentina, Uruguay and Venezuela 3,20 , as well as in European countries, such as France and the United Kingdom, in addition to West Asian countries like Turkey (Figure 1). All strains within this clade, with the exception of MK754263/G/URY/2018, share the same amino acid substitutions at positions I671V and A702S. In order to confirm that the clade D3-VS was not an artifact of convergent amino acid substitutions, a novel tree was inferred excluding codons 671 and 702. The resulting tree displayed the same topology of the original one, maintaining an elevated statistical support (hLRT=0.935) (data not shown).

DISCUSSION
During 2018, major HFMD outbreaks occurred throughout Brazil 9,10 . We tested 49 samples collected at the Pediatric Complex Arlinda Marques, Paraiba State, Brazilian Northeastern region, during the outbreak period, and CV-A6 was identified. Approximately three-fourth of the CV-A6 patients had viremia detected during the acute phase of disease. Cheng et al. 21 in a multi-center cohort study conducted in Taiwan showed that viremia associated with EV-71 occurred more frequently in children of one year of age in HFMD patients, agreeing with the results found in the present study. Cheng et al. 21 also demonstrated that after the third day of disease onset, the proportion of patients with viremia decreased to around one-fourth and no viremia was detected after the seventh day of disease onset, even if patients had persistent clinical symptoms, indicating that the resolution of EV viremia occur within seven days. Our study found a similar trend.
It is well known that RT-qPCR provides rapid diagnosis and increases the sensitivity for EV detection 22 ; however, our results have shown that conventional semi-nested RT-PCR was more sensitive than RT-qPCR. In fact, this is not an entirely surprise, as conventional PCR systems have often shown to be more sensitive than real-time PCR systems in the detection of different pathogens, including viruses 23,24 . Both methods targeted the EV 5'UTR region, but the sensitivity of the conventional PCR could be increased by adding a second round of amplification. On the other hand, conventional PCR methods are more prone to errors, including the lack of amplification products due to higher detection limits depending on the amount of cDNA per reaction 25 . Traditional detection of amplified DNA relies upon electrophoresis of nucleic acids in the presence of ethidium bromide (or similar) and visual analysis of the resulting bands in an ultraviolet light transillumination device 26 . Consequently, some fragments are not clearly detectable on agarose gels since they are visualized as faint bands or are not seen when the template cDNA concentrations was insufficient 25 , including the three EV samples (9B, 27B and 28B) recognized exclusively by RT-qPCR even when other samples with higher Ct values (i.e., samples 4B, 13B, 16B and 17B) were positive by both methods.
CV-A6 caused a massive HFMD outbreak in Finland in 2008, and then spread to other European countries and to North America [27][28][29][30] . Most of the largest CV-A6 HFMDassociated outbreaks have been reported in Asian Pacific countries, and CV-A6 has replaced EV-71 and CA-16 as the main enterovirus causing HFMD in the region since 2009 [31][32][33][34] . In Latin America, CV-A6 associated with HFMD cases began to circulate between 2011 and 2012, with reports from Cuba 35 and Brazil 6 , reaching Peru in 2016 36 . Since 2018, CV-A6 has been more frequently detected in HFMD outbreaks, as well as in sporadic cases in the region 3,5,7,37 . The present study reinforces the importance of the emergence of CV-A6 causing HFMD in Latin America. CV-A6 is additionally characterized by an atypical HFMD clinical presentation with lesions in other anatomical areas, such as the trunk, limbs and buttocks 34 , supporting the clinical data observed here.
EV typing is crucial for assessing the enteroviral biodiversity associated with HFMD infections and also to improve our understanding on their potential public health burden. A limitation of the present report was the failure to genotype all EV-positive samples. Simultaneous detection and identification of all EV genotypes in clinical samples is very challenging due to its high genetic diversity 20 . EV genotyping mainly depends on Sanger sequencing targeting the VP1 region 38 which clearly impaired our capacity to identify all EV-positive samples.
Brazilian CV-A6 strains detected here are genetically linked to the strains concurrently circulating in Latin America and European countries 3,37 . Brazilian CV-A6 strains BRA-38B, BRA-41B and BRA-48B are also distantly related to those CV-A6 strains that circulated earlier in Brazil, from 2012 to 2017 6,13 . These data suggest that the Brazilian CV-A6 strains detected in 2018 in Paraiba State and the CV-A6-associated HFMD strains recently reported in Latin America countries, may have a common origin and did not seem to have emerged through the evolution of local CV-A6 strains that circulated previously in Brazil, therefore suggesting a global co-circulation of a set of different CV-A6 strains introduced in the country at different times. Moreover, the sub-lineage D3 became the main CV-A6 sub-lineage since 2008 associated not only to HFMD outbreaks 6,7,34 , but also to meningitis 13 and severe acute respiratory illness 19 . The sub-lineage D3 may have a stronger ability of transmission and it has been suggested that this sub-lineage can be responsible for the growing circulation of CV-A6 worldwide 39 .
Of outstanding significance in the present study was the detection and molecular characterization of EV strains in serum samples. It is well known that the detection of EV exclusively in one serum sample must be carefully interpreted. Single serum samples alone are not suitable for diagnosing EV-associated HFMD. In general, for the EV-associated HFMD diagnosis, it is recommended that oropharyngeal swabs or stool samples are accompanied by paired blood samples collected in the acute-phase and in the convalescent-phase of the disease. The diagnosis of EV is only confirmed when serological tests (neutralization test) present with a significant increase in the antibody titers (three times or more) in the paired serum samples, and these antibodies are raised to the same EV serotype isolated from feces or oropharyngeal swabs 40 . During the HFMD outbreak in Paraiba State reported here, only the acute-phase serum samples were collected, therefore they were the only samples available for the EV diagnosis. The lack of oropharyngeal swabs impaired the testing of EV in both samples in order to evaluate if blood EV PCR testing in patients with suspected HFMD is additionally favored when compared with PCR testing performed on oropharyngeal specimen. Despite this limitation, we have shown that serum samples collected in the acute phase of HFMD are suitable for the diagnosis of EV-associated HFMD by rapid molecular methods, especially during outbreaks emergencies.

CONCLUSION
In conclusion, our findings showed that CV-A6 was involved in the Brazilian HFMD outbreaks that occurred during 2018. The increasing circulation of the emerging CV-A6 associated with HFMD, along with the detection of more severe cases worldwide, suggest the need for a more intense surveillance system of HFMD in Brazil. Our data also support the use of the RT-qPCR method as a powerful tool for a rapid EV diagnosis in blood samples during HFMD outbreaks. The early and timely EV detection in serum samples could positively affect the patients' management by reducing unnecessary antibiotic and antiviral therapies. In addition, the serum EV genotyping has the potential to assist health care planners to improve the health outcomes in EV-associated HFMD disease.