Coleta de dados
A matriz de expressão gênica e os dados clínicos correspondentes dos conjuntos de dados GSE142025 e GSE96804, relacionados à DN, foram obtidos do banco de dados Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/)15. O conjunto de treinamento (GSE142025) incluiu amostras de tecido renal de 27 pacientes com DN e nove controles, sequenciados utilizando a plataforma GPL20301. O conjunto de validação (GSE96804) compreendeu dados de sequenciamento de 41 pacientes com DN e 20 controles, processados utilizando a plataforma GPL17586. Ambos os conjuntos de dados concentram-se em tecido renal; o conjunto de dados GSE96804 examina especificamente o glomérulo, a unidade primária de filtração do rim (Figura 1). O conjunto de dados GSE142025 (conjunto de treinamento) compreende amostras de tecido renal inteiro e fornece uma visão abrangente da paisagem transcritômica da DN. O conjunto de dados GSE96804 (conjunto de validação), por outro lado, foca especificamente no tecido glomerular, que é o local primário da lesão na filtração glomerular. Como esses dois conjuntos de dados não foram diretamente mesclados devido às diferenças na plataforma e no tecido, a correção de efeito de lote não foi aplicada. Em vez disso, a validação cruzada entre conjuntos de dados foi realizada de forma independente. Genes com alterações direcionais consistentes e uma área sob a curva (AUC) maior que 0,7 em ambos os conjuntos de dados foram selecionados como candidatos robustos, o que apoia a generalização entre compartimentos renais.
Um total de 1.136 genes relacionados às mitocôndrias (MRGs) foi extraído do banco de dados MitoCarta3.0 (https://www.broadinstitute.org/mitocarta). Os ingredientes ativos do LLF foram previstos utilizando o banco de dados Traditional Chinese Medicine Systems Pharmacology (TCMSP) (http://sm.nwsuaf.edu.cn/lsp/tcmsp.php), com base em um limiar de biodisponibilidade oral (OB) de ≥30% e um limiar de semelhança com fármacos (DL) de ≥0,18. Posteriormente, os genes-alvo potenciais para os constituintes ativos foram previstos utilizando o banco de dados Swiss Target Prediction (http://www.swisstargetprediction.ch/).
Análise de expressão diferencial
A análise de expressão diferencial de GSE142025 (DN versus controle) foi realizada utilizando o pacote limma (v3.54.1), com critérios de significância definidos em P.adj < 0,05 e |log2FoldChange| > 0,516. Os gráficos de vulcão e mapas de calor foram visualizados utilizando os pacotes ggplot2 (v 3.3.6) e ComplexHeatmap (v 2.14.0)17,18, respectivamente. Os genes diferencialmente expressos (DEGs), genes relacionados a resistência (MRGs) e os genes-alvo potenciais dos ingredientes ativos foram intersectados, e os genes sobrepostos foram definidos como genes candidatos. A rede que liga os ingredientes ativos aos genes candidatos foi construída utilizando o software Cytoscape (v 3.9.0)19.
Análise de enriquecimento funcional e construção da rede de interação proteína-proteína (PPI)
As análises de enriquecimento da Ontologia Genética (GO) e do Banco de Dados Kyoto de Genes e Genomas (KEGG) para os genes candidatos foram realizadas utilizando o pacote clusterProfiler (v 4.6.2) para explorar suas funções biológicas e vias de sinalização associadas (P.adjust < 0,05). Os genes candidatos foram então inseridos no banco de dados STRING (https://cn.string-db.org/) para recuperar as interações de PPI (nível de confiança ≥ 0,4), seguido pela construção de uma rede de PPI utilizando o Cytoscape (v 3.9.0)20.
Aprendizado de máquina
Quatro algoritmos de aprendizado de máquina, incluindo floresta aleatória (RF), k-vizinhos mais próximos (KNN), mínimos quadrados parciais (PLS) e máquina de vetores de suporte com kernel base radial (SVM), foram implementados usando o pacote caret (v6.0-93) com base no conjunto de dados GSE142025. Os genes candidatos identificados na análise anterior foram utilizados como variáveis preditoras, e o status da doença (DN ou controle) foi utilizado como resultado. Para o modelo KNN, foi implementada uma validação cruzada de 10 dobras usando a função trainControl, com tuneLength = 10. O modelo RF foi ajustado com 20 árvores (ntree = 20); os modelos PLS e SVM foram ajustados usando as configurações implementadas no fluxo de trabalho original do caret.
As distribuições residuais do modelo foram avaliadas utilizando o pacote DALEX (v2.4.3). Em seguida, foi realizada a avaliação de variáveis baseada em permutação usando a função variable_importance do DALEX, com o erro quadrático médio (RMSE) especificado como função de perda. O dropout_loss resultante representa o RMSE obtido após a permutação da variável correspondente21,22. No fluxo analítico original, as variáveis com dropout loss < 0,281 em todos os modelos foram mantidas, e os genes comuns aos quatro modelos foram definidos como genes de características consenso para validação subsequente.
Os quatro modelos de aprendizado de máquina foram utilizados principalmente para priorização de características, e não para a construção de um classificador clínico final. Assim, a discriminação diagnóstica foi posteriormente avaliada ao nível do gene individual por meio de análise ROC nos conjuntos de dados de descoberta e de validação.
Identificação de biomarcadores candidatos
Em GSE142025 e GSE96804, as diferenças de expressão dos genes característicos entre amostras DN e controles foram avaliadas utilizando o teste de Wilcoxon. Apenas genes com expressão diferencial significativa (P < 0,05) e tendências de expressão concordantes nos dois conjuntos de dados foram selecionados para análise da curva característica de operação do receptor (ROC). O pacote pROC (v1.18.0) foi utilizado para gerar as curvas ROC e calcular a área sob a curva (AUC), sendo os genes com AUC > 0,7 em ambos os conjuntos de dados classificados como biomarcadores candidatos23.
Análise de enriquecimento de conjuntos de genes (GSEA)
As funções biológicas e as vias de sinalização associadas aos biomarcadores candidatos foram exploradas mais a fundo utilizando a análise de enriquecimento genômico (GSEA) no conjunto de dados GSE142025. Primeiramente, foi realizada uma análise de correlação de Spearman entre os biomarcadores candidatos e todos os demais genes utilizando o pacote psych (v2.2.9).24. Os coeficientes de correlação foram calculados e ordenados (do mais alto para o mais baixo). O conjunto de genes de referência utilizado foi c2.cp.kegg.v2023.1.Hs.symbols.gmt do Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/). Em seguida, foi realizada a análise GSEA (enriquecimento de conjuntos de genes) para avaliar o enriquecimento dos genes ordenados no conjunto de genes de referência, utilizando o pacote clusterProfiler (v4.6.2). Foi aplicada a correção para testes múltiplos. via o método FDR e ajustado P-valores (denotados como P.ajustado) foram considerados significativos se < 0.05.
Análise de modificações m6A
Para investigar modificações de metilação no RNA dos biomarcadores candidatos, utilizou-se o banco de dados SRAMP (http://www.cuilab.cn/sramp/) para prever os sítios de modificação m6A nos biomarcadores candidatos, com foco em posições de alta confiança dentro de suas estruturas secundárias. Em seguida, utilizou-se o banco de dados ENCORI (https://starbase.sysu.edu.cn/) para identificar proteínas modificadas por m6A que interagem com os biomarcadores candidatos, empregando o parâmetro |HepG2 (shRNA)| > 1 para selecionar as proteínas-chave. Posteriormente, utilizou-se o banco de dados RPISeq (http://pridb.gdcb.iastate.edu/RPISeq/) para prever a probabilidade de interações entre as proteínas-chave e os biomarcadores candidatos. As sequências de RNA de ambos foram enviadas em formato de texto simples para gerar pontuações de predição pelos classificadores RF e SVM. Considerou-se significativa uma interação quando a pontuação excedeu 0,525. A análise no SRAMP foi realizada com o limiar definido como ‘Alto’ para a predição de sítios de m6A, utilizando o modo ‘Transcript’ e parâmetros padrão. A análise no ENCORI utilizou a função de interação ‘miRNA-mRNA’ com o parâmetro ‘HepG2 (shRNA)’ > 1. A análise no RPISeq utilizou o classificador RF com parâmetros padrão; pontuações > 0,5 indicaram interação positiva. Trata-se de predições computacionais, não de evidência experimental de modificação m6A ou de interações RNA-proteína em tecido renal. O critério HepG2 shRNA foi obtido a partir de conjuntos de dados pré-calculados do ENCORI e pode não refletir regulação específica do rim.
Análise da infiltração imune
O algoritmo CIBERSORT foi aplicado para estimar as proporções de 22 tipos de células imunes em amostras de controle e de DN provenientes de GSE142025, com visualização por meio de um mapa de calor gerado usando o pacote ggplot2 (v3.3.6)26. O CIBERSORT foi executado utilizando a matriz de assinatura LM22, com 1.000 permutações e normalização por quantis desativada (como recomendado para dados de microarranjo). Amostras com p < 0,05 no CIBERSORT foram mantidas para análises adicionais. O CIBERSORT estima frações de células imunes a partir da expressão em tecido renal total, o que não permite distinguir infiltração específica por compartimento (por exemplo, glomerular versus tubulointersticial) nem diferenciar leucócitos infiltrantes de células imunes residentes. Assim, as correlações relatadas são ao nível do tecido e devem ser validadas por métodos espaciais. A análise de correlação de Spearman entre células imunes diferencialmente abundantes e biomarcadores candidatos foi realizada utilizando o pacote psych.
Construção de redes e docking molecular
MicroRNAs (miRNAs) que interagem com biomarcadores candidatos foram previstos utilizando o banco de dados miRNet (https://www.mirnet.ca). Posteriormente, RNAs longos não codificantes (lncRNAs) que direcionam os miRNAs identificados foram previstos por meio dos bancos de dados TarBase (http://www.diana.pcbi.upenn.edu/tarbase) e starbase (http://starbase.sysu.edu.cn/). Os lncRNAs comuns a ambos os bancos de dados foram selecionados para a construção da rede. Uma rede regulatória lncRNA-miRNA-mRNA foi então construída utilizando o software Cytoscape. Ingredientes ativos potenciais que direcionam biomarcadores candidatos foram escolhidos para construir uma rede ingrediente ativo-biomarcador. Além disso, ingredientes ativos, biomarcadores candidatos e vias identificadas na GSEA foram incorporados ao Cytoscape para criar uma rede ingrediente ativo-biomarcador-via.
A análise de docking molecular foi realizada para avaliar a afinidade de ligação entre ingredientes ativos potenciais e biomarcadores candidatos. As estruturas tridimensionais das proteínas biomarcadoras foram obtidas do Protein Data Bank do Research Collaboratory for Structural Bioinformatics (RCSB PDB, https://www.rscb.org/pdb) no formato de arquivo PDB. As estruturas bidimensionais dos ingredientes ativos potenciais foram recuperadas no formato SDF a partir do banco de dados PubChem (http://pubchem.ncbi.nlm.nih.gov). O docking molecular foi conduzido utilizando a plataforma CB-Dock (http://clab.labshare.cn/cb-dock/php/blinddock.php). Uma energia de ligação inferior a -5 kcal/mol indicou uma afinidade de ligação forte27.
Preparação e autenticação de Ligustri Lucidi Fructus
Aqui, Ligustri Lucidi Fructus (LLF) refere-se ao fruto maduro seco de Ligustrum lucidum W. T. Aiton (Oleaceae). O material botânico foi autenticado de acordo com a Farmacopeia Chinesa, e um espécime-testemunha nº 20240506,20240911,20241103 foi depositado na Universidade de Medicina Tradicional Chinesa de Shanxi.
Para a preparação da decocção, 200 g de fatias qualificadas de LLF foram embebidas em 1.000 mL de água destilada por 30 min à temperatura ambiente. A mistura foi fervida vigorosamente e depois mantida em fervura suave por 60 min. O filtrado foi coletado, e os resíduos das plantas foram submetidos novamente à decocção com outros 1.000 mL de água destilada por 60 min. Os dois filtrados foram combinados, filtrados, centrifugados e concentrados sob pressão reduzida para obter uma concentração final de estoque de 1 g de droga bruta/mL (volume total de 100 mL). A decocção preparada foi armazenada a 4 °C para uso de curto prazo ou a −20 °C para conservação de longo prazo. A qualidade do LLF e de sua decocção foi rigorosamente identificada e verificada de acordo com os padrões da Farmacopeia Chinesa para garantir a confiabilidade e reprodutibilidade experimentais.
Para identificação qualitativa, realizou-se a cromatografia em camada delgada. Resumidamente, um volume apropriado da decocção preparada foi centrifugado e o sobrenadante foi extraído com metanol. Após a filtração, a solução da amostra e a solução padrão de referência de especnuezhenídeo foram aplicadas na mesma placa de gel de sílica G. Após a eluição, secagem e inspeção sob luz ultravioleta, a mancha da solução da amostra apresentou cor e posição consistentes com o composto de referência, confirmando a presença do componente bioativo característico do LLF.
Para controle de qualidade quantitativo, foi realizada a detecção por cromatografia líquida de alta eficiência. A análise foi realizada utilizando uma coluna C18 com metanol–água como fase móvel. O comprimento de onda de detecção foi ajustado em 224 nm. O teor de especnuezhenídeo na decocção foi determinado com base na curva padrão. Os resultados demonstraram composição química estável e uniforme da decocção preparada, garantindo qualidade constante do medicamento durante todo o experimento de intervenção em animais.
Experimentos com animais
Doze camundongos machos db/db de grau SPF (8–9 semanas de idade) e seis camundongos db/m de idade correspondente foram mantidos na instalação animal SPF da Universidade de Medicina Tradicional Chinesa de Shanxi. Antes dos experimentos, os animais foram aclimatados por 7 dias em ciclo de 12 h de luz/12 h de escuro com acesso ad libitum a alimento e água. O estudo foi aprovado pelo Comitê de Ética da Universidade de Medicina Tradicional Chinesa de Shanxi (Número de Aprovação: 2022DW167). Animais que apresentaram perda de peso corporal >20%, condição moribunda ou incapacidade de acessar alimento ou água foram eutanasiados de forma humanitária antes do término planejado do estudo. Ao final do estudo, todos os camundongos remanescentes foram eutanasiados por injeção intraperitoneal de pentobarbital sódico seguida de fratura cervical, conforme os protocolos institucionais.
Após o período de aclimatização, o estabelecimento de um modelo de DN em camundongos db/db foi confirmado por meio de nível de glicose no sangue da veia da cauda ≥ 16 mmol/L e microalbuminúria, conforme indicado por uma tira reagente positiva para microalbumina urinária. Após a confirmação bem-sucedida do modelo de DN, os camundongos db/db foram divididos aleatoriamente em dois grupos (n = 6 por grupo): o grupo modelo de DN (DN) e o grupo tratado com LLF (Tratamento). Além disso, camundongos db/m (n = 6) foram utilizados como grupo controle (Controle). A dose foi escolhida com base em estudos farmacodinâmicos prévios de LLF em ratos diabéticos e foi convertida para a dose equivalente em humanos utilizando a normalização pela área de superfície corporal28. Os grupos Controle e DN receberam água destilada, enquanto o grupo tratado recebeu 3,5 g/kg de LLF durante 8 semanas. Após 8 semanas de administração, todos os camundongos foram eutanasiados para coleta de soro, urina e tecido renal para testes subsequentes.
Indicadores sanguíneos e urinários
Os níveis séricos de glicose foram analisados utilizando um Analisador Automático de Bioquímica Sanguínea. A concentração de microalbumina urinária foi medida de acordo com as instruções do kit (Arquivo Suplementar 1).
Observação patológica de tecidos renais de camundongo
Os tecidos renais foram processados para exame histopatológico. Após fixação em paraformaldeído a 4%, os tecidos foram lavados, desidratados, emblocados em parafina e seccionados. Em seguida, foi realizada a coloração com hematoxilina e eosina (HE), e as lâminas coradas foram examinadas ao microscópio óptico para avaliação das alterações patológicas.
Reação em cadeia da polimerase quantitativa com transcriptase reversa (RT-qPCR)
A RT-qPCR foi utilizada para determinar a expressão dos biomarcadores candidatos em tecido renal de camundongo. O RNA total foi extraído de acordo com as instruções do fabricante, após o que a concentração e qualidade do RNA foram avaliadas (Tabela 1). O cDNA foi sintetizado a partir do RNA extraído utilizando o Kit de Síntese de cDNA. A amplificação foi realizada com os pares de primers listados na Tabela 1, utilizando GAPDH como gene de referência. Os níveis relativos de expressão foram calculados pelo método 2−ΔΔCt13,26.
Análise estatística
Todas as análises estatísticas foram realizadas utilizando o software R (versão 4.2.2) juntamente com o software necessário para os procedimentos experimentais correspondentes. Salvo indicação em contrário, todos os testes estatísticos foram bicaudais, e as diferenças foram consideradas estatisticamente significativas em P < 0,05. A análise de expressão diferencial transcriptômica foi realizada utilizando o pacote limma. Genes com P ajustado < 0,05 e mudança absoluta no log2 da razão maior que 0,5 foram definidos como diferencialmente expressos.
Para comparações da expressão de genes candidatos entre amostras independentes de DN e controles, foram utilizados testes não paramétricos de postos de Wilcoxon, conforme indicado no fluxo analítico original. As correlações entre biomarcadores candidatos e frações de células imunes foram avaliadas utilizando a correlação de postos de Spearman.
Os dados experimentais são apresentados como média ± DP. As comparações entre os três grupos independentes de animais foram realizadas utilizando análise de variância unidirecional quando as premissas da análise paramétrica foram atendidas. Para comparações post hoc, o teste da menor diferença significativa foi utilizado quando as variâncias eram homogêneas, enquanto o teste T3 de Dunnett foi utilizado quando as variâncias eram desiguais. As medidas farmacodinâmicas basais e na semana 8 foram analisadas e apresentadas separadamente; nenhuma inferência sobre interação grupo-por-tempo foi feita. Os scripts em R e os dados brutos utilizados para as análises de bioinformática e aprendizado de máquina estão disponíveis no Arquivo Suplementar 1.