Em formação

21: Redes regulatórias - inferência, análise, aplicação - biologia


21: Redes regulatórias - inferência, análise, aplicação

Inferir progressão temporal latente e redes regulatórias a partir de dados transcriptômicos transversais de amostras de câncer

Desvendar redes regulatórias moleculares subjacentes à progressão da doença é extremamente importante para a compreensão dos mecanismos da doença e identificação de alvos de drogas. Os métodos existentes para inferir redes regulatórias de genes (GRNs) dependem principalmente de dados de expressão de genes de curso no tempo. No entanto, a maioria dos dados ômicos disponíveis de estudos transversais de pacientes com câncer muitas vezes carecem de informações temporais suficientes, levando a um desafio importante para a inferência de GRN. Por meio da quantificação da progressão latente usando distâncias múltiplas baseadas em caminhadas aleatórias, propomos um método Bayesiano baseado em progressão latente-temporal, PROB, para inferir GRNs a partir de dados transcriptômicos transversais de amostras de tumor. A robustez do PROB às variabilidades de medição nos dados é matematicamente comprovada e verificada numericamente. A avaliação de desempenho em dados reais indica que PROB supera outros métodos em inferência de pseudotempo e inferência de GRN. Aplicações para câncer de bexiga e câncer de mama demonstram que nosso método é eficaz para identificar os principais reguladores da progressão do câncer ou alvos de drogas. O ACSS1 identificado é validado experimentalmente para promover a transição epitelial para mesenquimal de células de câncer de bexiga, e as interações de alvos FOXM1 previstas são verificadas e são preditivas de recidiva no câncer de mama. Nosso estudo sugere novas maneiras eficazes de modelagem de dados transcriptômicos clínicos para caracterizar a progressão do câncer e facilita a tradução de abordagens baseadas em redes regulatórias em medicina de precisão.

Declaração de conflito de interesse

Os autores declararam que não existem interesses conflitantes.

Bonecos

Fig 1. Ilustração da estrutura PROB ...

Fig 1. Ilustração da estrutura PROB para inferir a rede regulatória de genes causais de ...

Fig 2. Demonstrando a robustez do PROB usando ...

Fig 2. Demonstrando a robustez do PROB usando conjuntos de dados sintéticos em diferentes níveis de variabilidades.

Fig 3. Comparação de PROB com outro ...

Fig 3. Comparação de PROB com outros métodos de inferência de pseudotempo e métodos de inferência GRN existentes ...

Fig 4. Reconstruindo redes regulatórias de EMT durante ...

Fig 4. Reconstruindo redes regulatórias de EMT durante a progressão do câncer de bexiga.

( uma ) Padrões de expressão ...

Fig 5. Validação experimental do previsto ...

Fig 5. Validação experimental do papel previsto de ACSS1 em EMT de câncer de bexiga.

Fig 6. FOXM1 foi revelado como um ...

Fig 6. FOXM1 foi revelado como um gene-chave subjacente à progressão do câncer de mama por PROB.


CÊNICO: inferência e clustering de rede regulatória de célula única

Apresentamos o SCENIC, um método computacional para reconstrução simultânea da rede reguladora de genes e identificação do estado da célula a partir de dados de RNA-seq de uma única célula (http://scenic.aertslab.org). Em um compêndio de dados de uma única célula de tumores e cérebro, demonstramos que a análise cis-regulatória pode ser explorada para orientar a identificação de fatores de transcrição e estados celulares. SCENIC fornece percepções biológicas críticas sobre os mecanismos que conduzem a heterogeneidade celular.

Declaração de conflito de interesse

Concorrência de interesses financeiros

Os autores declaram não haver interesses financeiros concorrentes.

Bonecos

Figura 1. O fluxo de trabalho SCENIC e seu ...

Figura 1. O fluxo de trabalho SCENIC e sua aplicação ao cérebro do mouse.

Figura 2. Comparação entre espécies de redes neuronais ...

Figura 2. Comparação entre espécies de redes neuronais e tipos de células.

Figura 3. SCENIC supera os efeitos do tumor e ...

Figura 3. SCENIC supera os efeitos do tumor e desvenda os estados celulares relevantes e GRNs no câncer.


Resultados

Formulação do método proposto

LASSO é uma regressão linear que penaliza a soma do valor absoluto dos coeficientes de regressão 49. Baseia-se na combinação de eu2-norma, ou seja, a soma residual dos quadrados, com o eu1-norma dos coeficientes de regressão, que equivale a esparsidade devido à redução dos coeficientes para zero. Além disso, a proposta de considerar a fusão na formulação clássica LASSO pretendia abordar problemas com uma ordem significativa nas características consideradas (ou seja, regressores). No LASSO fundido, minimização do eu1A norma é imposta não apenas aos coeficientes de regressão (ordenados), como em LASSO, mas também às diferenças consecutivas dos coeficientes de regressão com base na ordem assumida dos regressores correspondentes 50.

Nesta seção, formulamos a abordagem que estende o conceito de LASSO fundido, incorporando informações sobre a similaridade do comportamento diferencial para os genes de resposta e regressor. Aqui, cada gene é tomado como resposta e os genes que codificam fatores de transcrição (conhecidos como genes TF) são usados ​​como regressores. A fusão impõe a restrição de similaridade nas redes obtidas de cada conjunto de dados. Dessa forma, a abordagem garante a dispersão, comumente observada em redes regulatórias de genes, com evidências baseadas em dados biologicamente significativas para as relações inferidas. Além disso, apresentamos o cenário experimental e o tipo de dados nos quais a abordagem é prontamente aplicável.

Cenário experimental e matrizes de peso

Para consistência ao longo do texto, todos os sobrescritos referem-se aos experimentos / condições e os subscritos denotam genes. Vamos definir a abordagem para P genes usados ​​como regressores e um único gene usado como resposta. A abordagem pressupõe que k conjuntos de dados de expressão gênica, denotados por X eu , estão reunidos sob k diferentes condições ao lado da referência correspondente (ao controle) conjuntos de dados X c,eu , 1 ≤ euk, todos eles contendo os níveis de expressão de P genes acabados N eu , 1 ≤ euk, pontos no tempo ou perturbações (Fig. 1, etapas 1 e 2). Observe que os conjuntos de dados do k diferentes condições não precisam estar no mesmo domínio de tempo ou com a mesma frequência de amostragem, o único requisito é que cada conjunto de dados permita obter com segurança a matriz de covariância correspondente para o P genes. Além disso, os perfis Y eu para o único gene de resposta sobre o k condições juntamente com os perfis de controle correspondentes Y c,eu , 1 ≤ euk são dados.

Representação do fluxo de análise, conjuntos de dados e sua transformação para uso no modelo.

Representação de conjuntos de dados de condição e controle X eu e X ao controle , respectivamente, contendo os níveis de expressão de P genes, como regressores, sobre N eu pontos de tempo, 1 ≤ euk. Matriz de bloco X contém os valores de expressão do gene X eu , 1 ≤ euk, para os genes regressores (fatores de transcrição) na diagonal. Y inclui os perfis Y eu , 1 ≤ euk, para os genes de resposta. Matrizes de peso C eu contêm a similaridade de comportamento diferencial entre o gene de resposta e cada um dos regressores com base em conjuntos de dados de k experimentos de perturbação e os controles correspondentes. Valores do coeficiente de regressão β eu , 1 ≤ euk, consiste na relação regulatória entre um FT e seus alvos de k experimentos de perturbação.

Uma vez que os conjuntos de dados de referência X c,eu , 1 ≤ euk, estão disponíveis, primeiro determinamos o comportamento diferencial para cada um dos regressores, bem como para a resposta em cada condição. Para este fim, contamos com a estatística B 11 específica do gene correspondente às probabilidades logarítmicas de que o gene é expresso diferencialmente em um determinado ponto de tempo de uma condição particular em comparação com o controle. Let, 1 ≤ jP, 1 ≤ euk, 1 ≤ tN eu seja a probabilidade de que o gene j no ponto do tempo t sob condição eu é expresso diferencialmente. Então a matriz Pr eu de dimensão reúne as probabilidades para o comportamento diferencial dependente do tempo dos genes considerados, estimados pelas estatísticas B correspondentes. A estatística B foi estimada para cada gene e em cada ponto de tempo, comparando o conjunto de dados X eu do tratamento e do respectivo controle X c,eu , 1 ≤ euk.

As probabilidades derivadas podem ser usadas para definir matrizes de peso para cada condição eu (1 ≤ euk), denotado por, que captura as informações sobre as semelhanças entre o gene de resposta e cada um dos genes regressores com base em seu comportamento diferencial (Fig. 1, etapa 3), usando o seguinte:

onde e são as probabilidades de comportamento diferencial para o gene de resposta Y e o gene regressor j (1 ≤ jP), respectivamente, no ponto de tempo / perturbação t.

Se o valor de for próximo de zero, o j º O gene regressor tem, em média, um comportamento diferencial semelhante à resposta em todos os pontos de tempo / perturbações em um conjunto de dados de uma condição. Mais especificamente, está próximo de zero se ambos os genes forem expressos diferencialmente ou se ambos não forem afetados pela perturbação. Devido à natureza simétrica do valor absoluto na Eq. (1), a matriz de peso C eu é simétrico com zeros ao longo da diagonal.

Modelo baseado em regressão

Tendo k conjuntos de dados de expressão gênica resolvidos no tempo, incluindo P + 1 genes e N eu , 1 ≤ euk, pontos de tempo juntamente com o (s) conjunto (s) de dados de controle correspondente (s), pretendemos formular um modelo que capture os seguintes três critérios: (1) os níveis de expressão de cada gene devem ser explicados pelos níveis de expressão de um pequeno número de fator de transcrição (TF) genes codificadores, ou seja, a regressão correspondente deve ser esparsa para reduzir o número de falsos positivos e provavelmente aumentar a detecção de relações diretas (2) as redes regulatórias devem ser inferidas simultaneamente sobre o dado k conjuntos de dados para explicar simultaneamente todas as condições analisadas e (3) uma ligação direta devem ser preferidos para genes que exibem comportamento diferencial semelhante sobre as condições consideradas, uma vez que um gene de comportamento diferencial (TF) é susceptível de alterar o comportamento de um alvo direto. Aqui, com o propósito de construir os modelos de regressão, nos concentramos em usar os conjuntos de dados do k condições apenas, uma vez que são potencialmente mais informativos sobre os genes responsivos em comparação com o estado de referência. Além disso, os conjuntos de dados de controle entram na modelagem por meio das probabilidades de comportamento diferencial e são, portanto, também empregados na reconstrução da rede regulatória gênica.

O primeiro critério pode ser capturado pelo LASSO clássico, de forma que os regressores com coeficientes de regressão diferentes de zero são considerados envolvidos em uma relação direta com o gene de resposta, gerando a rede regulatória do gene. Para abordar o segundo critério de inferir simultaneamente k redes do k conjuntos de dados diferentes, deve-se garantir que: (eu) a integração e transformação dos dados devem ser realizadas de tal forma que a penalidade LASSO seja imposta simultaneamente a todos os conjuntos de dados e (ii) a k redes reconstruídas devem ser o mais próximo possível em termos de bordas (ou seja, relações) sua força dada pelos coeficientes e pelo sinal da relação (ativadora ou repressora).

A integração e transformação dos conjuntos de dados é realizada da seguinte forma: o k os conjuntos de dados transcriptômicos são combinados em uma matriz de bloco único X de dimensões kN × kP, Onde , Y é uma coleção de perfis para o gene de resposta de diferentes condições, resultando em um vetor de dimensão kN × 1. Como resultado, os coeficientes estimados formam um kP × 1-vetor correspondendo às bordas entre o gene de resposta Y e genes regressores no k rede. O vetor de coeficientes de regressão estimados é representado por β.

Para garantir que o k redes reconstruídas são o mais próximo possível, adicionamos o termo de fusão LASSO, onde β ′ = [β 1, β 2,…, β k − 1 ] T , β ′ = [β 2, β 3,…, β k ] T e a ordem do k conjuntos de dados são selecionados arbitrariamente. Em uma configuração de regressão, este termo de fusão impõe a restrição de que a soma das diferenças absolutas entre os coeficientes estimados do mesmo regressor ao longo dos conjuntos de dados consecutivos, com a ordem arbitrária assumida, seja minimizada (Fig. Suplementar S3 e Tabela Suplementar S8). Essa ideia difere da abordagem mais recente, em que redes reconstruídas a partir de diferentes conjuntos de dados são obtidas individualmente e posteriormente combinadas em uma rede de consenso pelas técnicas existentes 14,51.

O terceiro critério sobre a similaridade do comportamento diferencial entre a resposta e os regressores implica a inclusão da matriz de ponderação C eu , 1 ≤ euk, de modo que o regressor de maior poder explicativo, que está associado a coeficientes de regressão diferentes de zero, sobre os conjuntos de dados múltiplos são menos penalizados na regressão LASSO fundida. Isso pode ser alcançado multiplicando os coeficientes de regressão com a matriz de peso C (do tamanho kP × kP), reunindo a similaridade de comportamento diferencial entre cada um dos P regressores e o único gene de resposta sob o k condições. Como resultado, a expressão é incluída como uma penalidade modificada na regressão.

Portanto, o modelo final para reconstruir as interações regulatórias de genes é dado pela seguinte formulação de fusão LASSO sobre o k determinados conjuntos de dados (Fig. 1, etapa 4):

Onde Y é o gene de resposta que é regulado pelos regressores com coeficientes de regressão diferentes de zero.

Conjuntos de dados de expressão gênica

Escherichia coli

As respostas de Escherichia coli às condições de estresse já foram bem investigadas, resultando na caracterização dos componentes gerais e específicos da condição que regulam as mudanças transcricionais subjacentes ao ajuste a ambientes em mudança 52,53. Os conjuntos de dados coletados fornecem um excelente caso de teste com o qual o desempenho do método proposto e as alternativas concorrentes podem ser facilmente comparados.

Os conjuntos de dados de transcriptômica resolvidos no tempo coletados com a tecnologia de microarray foram obtidos a partir de 54, onde as mudanças na expressão de genes de E. coli cepa MG1655 foi monitorada sob quatro condições de estresse, incluindo: mudanças não letais de temperatura, ou seja, tratamento de calor e frio, estresse oxidativo (pela adição de peróxido de hidrogênio), deslocamento lactose-diauxic (ou seja, mudança da fonte primária de carbono) em relação às culturas cultivadas em condições ideais, conhecidas como controle. Todas as culturas foram cultivadas simultaneamente nas mesmas condições e diferentes perturbações foram aplicadas na fase intermediária inicial (DO 0,6). A amostragem foi realizada a partir de pontos de tempo 10–50 min pós-perturbação (em intervalos de 10 min) e dois pontos de tempo de controle antes de cada perturbação para todas as condições consideradas (os dados estão disponíveis em http: //www.ncbi.nlm.nih. gov / geo / query / acc.cgi? acc = GSE20305).

As razões para selecionar e usar esses conjuntos de dados são que as previsões podem ser prontamente validadas contra o padrão ouro de referência para E. coli fornecido pelo desafio DREAM5 14, bem como as interações de rede regulatória verificadas experimentalmente do RegulonDB 53. No RegulonDB, as informações sobre o efeito dos reguladores (em termos de ativação e repressão de interações) nos genes alvo também são fornecidas. Além disso, esses conjuntos de dados satisfazem os requisitos da abordagem proposta, uma vez que o controle único pode ser usado para estabelecer o comportamento diferencial do gene entre os conjuntos de dados mediante a aplicação dos estresses. Além disso, os dados são coletados no mesmo laboratório seguindo o mesmo protocolo, reduzindo assim o nível de ruído. Finalmente, e mais importante, buscamos usar conjuntos de dados do mundo real para estimar o desempenho real do método proposto em um ambiente realista, em vez de instâncias simuladas.

Mycobacterium tuberculosis

Mycobacterium tuberculosis (MTB) é uma bactéria patogênica cuja rede reguladora de genes é mal compreendida. No entanto, recentemente Galagan et al. 55 deu o primeiro passo na reconstrução da rede regulatória completa de MTB com base nos dados de expressão gênica ChIP-seq e microarray coletados sob condições de hipóxia e reaeração. Eles realizaram experimentos de transcriptômica resolvidos no tempo em que o nível de expressão dos genes foi medido em 1, 2, 3, 5 e 7 dias após a cultura MTB cepa H37Rv em condições de hipóxia bacteriostática. As amostras foram então colocadas de volta na cultura aeróbia de rolamento e os níveis de expressão gênica foram medidos após 1, 2, 4, 5 e 7 dias de re-aeração. Os níveis de expressão também foram medidos no ponto de tempo 0, que é considerado o conjunto de dados de controle para a análise diferencial (os dados estão disponíveis em http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE43466).

As razões para selecionar e usar esses conjuntos de dados são que: (i) as previsões podem ser validadas com relação a conjuntos de dados do mundo real recentes de espécies não-modelo, (ii) satisfaz os requisitos de nossa abordagem para inferir o gene regulador rede a partir de vários conjuntos de dados, já que as séries temporais reunidas podem ser consideradas como dois conjuntos separados de dados resolvidos no tempo (condições): hipóxia e reaeração, incluindo o ponto de tempo 0 como um controle e (iii) as previsões podem ser prontamente comparadas com a pequena sub-rede de MTB o que foi parcialmente verificado por meio de experimentos em 55 (Fig. 2 em 55).

Mus musculus

Mus musculus representa um organismo modelo para a compreensão da biologia humana e das doenças, portanto, estudar sua rede reguladora de genes em grande escala é uma tarefa desafiadora importante. Muitos experimentos foram realizados em células-tronco embrionárias (ES) de camundongos para explorar os detalhes de sua pluripotência e capacidade de se autopropagar e renovar. Para este fim, Sene et al. 56, mediu o nível de expressão de genes de três linhas de células ES de camundongo geneticamente distintas (R1, J1 e V6.5) durante a diferenciação em 11 pontos de tempo: 0 h, 6 h, 12 h, 18 h, 24 h, 36 h, 48 h, 4 d, 7 d, 9 d e 14 d. Os níveis de expressão no ponto de tempo 0 são considerados como o conjunto de dados de controle para a análise diferencial (os dados processados ​​foram baixados de http://www.maayanlab.net/ESCAPE/browse.php).

Este conjunto de dados fornece um bom caso de teste para comparar o desempenho das abordagens comparadas para a extração da rede reguladora de genes em um organismo superior. Além disso, a rede de células-tronco embrionárias baseada na literatura é usada como uma rede padrão global para validar as previsões obtidas a partir das abordagens comparadas (baixado do iScMiD (Integrated Stem Cell Molecular Interactions Database) http://amp.pharm.mssm.edu/ iscmid / Literatura / index.htm).

Análise comparativa

Para a análise comparativa, consideramos as abordagens mais recentes, bem como as do estado da arte, que incluem dois critérios, dispersão da rede e remoção de relacionamentos indiretos. Para tanto, utilizamos as seguintes abordagens: o silenciamento global 33, a deconvolução da rede 34, os Modelos Gráficos Gaussianos (GGM) 19, a informação mútua (ARACNE 25 e CLR 26), a Rede Bayesiana (catnet) 57, a abordagem GENIE3 21 e diferentes modelos baseados em regularização 45 (mais detalhes na seção Métodos).

Escherichia coli

As redes resultantes dessas abordagens foram comparadas à rede, incluindo relações regulatórias verificadas experimentalmente do RegulonDB e os padrões ouro do desafio DREAM5. As descobertas foram resumidas em termos das taxas de verdadeiro positivo e falso positivo resultantes, que resultaram nas curvas ROC (características de operação do receptor) correspondentes ilustradas na Fig. 2 (com base no valor limite usado para reter uma borda ponderada na rede). Usamos o pacote R mineto 58 para traçar as curvas ROC e as estatísticas obtidas usando este pacote estão resumidas na Tabela 1. Conforme ilustrado na Fig. 2, a abordagem proposta permitiu a previsão de relações regulatórias em uma porcentagem maior de verdadeiros positivos em comparação com o GENIE3, deconvolução de rede e os métodos de silenciamento global, enquanto GGM, ARACNE, CLR e modelos baseados em regularização tiveram um desempenho ruim para inferência precisa da rede regulatória. Além disso, é evidente que, nos conjuntos de dados utilizados, o método proposto com base na extensão do LASSO fundido teve desempenho semelhante às abordagens concorrentes consideradas quando apenas a conectividade (ou seja, a presença de uma borda / relacionamento) é considerada. Para quantificar ainda mais o desempenho dos métodos comparados, calculamos a área sob as curvas ROC (AUROCs), a área sob as curvas PR (precisão-recall) (AUPRs) e as taxas de verdadeiros positivos (TPR) em baixa taxa de falsos positivos (a maioria dos métodos comparados apresentou diferença máxima entre TPR e FPR em FPR = 0,03), apresentado na Tabela 1. Utilizou-se o pacote R pROC 59 para estimar as curvas AUROC e seus respectivos intervalos de confiança (IC) (Tabela Complementar S1). As outras estatísticas resultantes da comparação das abordagens consideradas estão resumidas na Tabela Suplementar S1. Ao todo, nossos resultados demonstram que o desempenho da abordagem proposta é estatisticamente semelhante ao desempenho do GENIE3 e da deconvolução da rede, enquanto supera o restante dos métodos concorrentes.

Curvas ROC para os métodos considerados na análise comparativa (E. coli conjuntos de dados).

As curvas ROC são mostradas para as redes regulatórias de genes previstas com base nos conjuntos de dados de (UMA) frio, (B) calor e (C) estresse oxidativo, bem como (D) deslocamento lactose-diauxic além do (E) combinação de todos os quatro conjuntos de dados usando os seguintes métodos: GGM, ARACNE, CLR, eu1/2, GENIE3, Silenciamento global, desconvolução da rede e a abordagem proposta com base no LASSO fundido que considera simultaneamente todos os quatro conjuntos de dados. Devido à alta similaridade no desempenho dos modelos baseados em regularização, as curvas ROC de eu1/2 modelo é ilustrado como o representante para evitar sobreposição de curvas. A cor das linhas tracejadas representa os métodos. TPR e FPR representam taxa de verdadeiro positivo e taxa de falso positivo, respectivamente.

Em nossa abordagem, a rede de todos os conjuntos de dados foi obtida tomando o máximo ou a média dos coeficientes de regressão correspondentes (consulte a Figura Suplementar S1 (a) para ilustração). Uma vez que a distribuição das diferenças entre os coeficientes de regressão obtidos com base em cada um dos quatro conjuntos de dados para a mesma resposta e regressores são quase idênticas (consulte a Figura Suplementar S1 (b)), não se espera que as duas redes sejam diferentes. Na verdade, a comparação entre as duas redes obtidas tomando a média e o máximo dos coeficientes de regressão não mostra uma diferença significativa em relação aos valores AUROC (p-valor = 0,6257, consulte a Tabela Suplementar S1).

O desempenho de cada método ao longo dos quatro conjuntos de dados e sua combinação foi resumido pela média geométrica para os cinco AUROCs e AUPRs correspondentes, enquanto a Pontuação Geral é a média dos AUROCscores e AUPRscores 14 calculados:

A pontuação geral na Eq. (3) pode ser usado para classificar os métodos comparados, de modo que um valor maior corresponda a um método de melhor desempenho. Conforme mostrado na Tabela 1, o desempenho do nosso método foi superior ao de outros métodos concorrentes com base no GlobalScore. Além disso, o desvio padrão das curvas AUROC para o método proposto é menor em comparação com os contendores (ver Tabela Suplementar S1). A estabilidade das redes obtidas suporta ainda o uso do termo fusão na abordagem proposta, que foi considerada com o propósito de produzir redes semelhantes em diferentes conjuntos de dados.

Em seguida, determinamos o menor peso normalizado no qual a primeira borda positiva verdadeira pode ser detectada - denominado valor do seletor os valores do seletor mais próximos de 1 indicam que as bordas positivas verdadeiras provavelmente têm pesos mais altos atribuídos. Conforme mostrado na Tabela 2, o método proposto apresentou o melhor desempenho para redes obtidas em todos os conjuntos de dados. Além disso, nossa abordagem está nos dois métodos de melhor desempenho ao considerar conjuntos de dados individuais. A discrepância pode ser explicada pela inclusão do termo fusão que insiste em relações regulatórias independentes de contexto. Como os pesos de borda das redes obtidas por ARACNE com cada um dos quatro conjuntos de dados eram idênticos, os valores do seletor correspondentes não são informativos e não podem ser usados ​​na comparação. Curiosamente, a deconvolução da rede e o silenciamento global não tiveram um bom desempenho com relação ao valor do seletor, apesar das reivindicações recentes sobre dados sintéticos. O razoavelmente alto valores do seletor obtidos a partir de modelos baseados em regularização confirmam ainda mais o poder das abordagens baseadas em regularização na atribuição de pontuações mais altas às verdadeiras arestas positivas.

Além disso, a precisão das redes regulatórias de genes inferidos depende da capacidade de prever o tipo de relações regulatórias - ativação ou repressão. Portanto, comparamos o desempenho dos métodos em relação ao tipo de relações regulatórias previstas. Para tanto, utilizamos a rede do RegulonDB que inclui no total 4.566 relacionamentos. A Tabela 3 resume as porcentagens das verdadeiras interações regulatórias de ativação e repressão previstas para os quatro conjuntos de dados e os métodos considerados. Além disso, inclui as percentagens das relações verdadeiras positivas e verdadeiras negativas previstas por diferentes métodos, independentemente do tipo de regulação.

ARACNE e CLR, bem como GENIE3, não são capazes de inferir o tipo de relações regulatórias. Conforme mencionado na Seção 3.2, os algoritmos originais para o silenciamento global e a deconvolução da rede não fornecem o sinal das interações, uma vez que visam classificar as interações regulatórias. Conforme mostrado na Tabela 3, nossos achados indicam que o desempenho do método proposto no que diz respeito à previsão do tipo de relações regulatórias é promissor. Além disso, o GGM e os modelos baseados em regularização resultaram na maior fração de arestas negativas verdadeiras, bem como na menor fração de arestas positivas verdadeiras. O método de deconvolução, no entanto, superou o resto das abordagens no que diz respeito à fração de bordas positivas verdadeiras para estresse de frio e oxidativo, bem como deslocamento lactose-diauxic (embora com valores de seletor muito menores, ver Tabela 2).

Além disso, para inspecionar especificamente a eficiência das redes inferidas, selecionamos uma sub-rede do RegulonDB incluindo quatro fatores sigma, RopD, RopE, RopH e RopS, incluídos nos conjuntos de dados analisados. Fatores sigma são proteínas necessárias para o início da síntese de RNA 60 e suas atividades dependem das condições ambientais. Enquanto RopD é o fator sigma primário que transcreve a maioria dos genes, as atividades dos outros três são específicas do ambiente, por exemplo, RopE é o fator sigma de estresse térmico extremo, RopH é o fator sigma de estresse térmico ativado mediante exposição ao calor e RopS é o fator sigma de fase estacionária. Curiosamente, para os quatro fatores sigma, nenhuma aresta foi prevista pelo GGM. Com relação às seis outras abordagens, bem como catnet (aplicada apenas em genes incluídos nas redes ilustradas), inspecionando a Fig. 3 (por uma questão de legibilidade, as redes resultantes para CLR, GENIE3, e modelos de regularização são mostrados na Fig. S2 suplementar), torna-se evidente que a abordagem proposta previu um número consideravelmente menor de falsos positivos, devido às restrições de esparsidade e um número comparável de positivos verdadeiros, devido à formulação baseada em regressão. Além disso, como resultado das restrições impostas na formulação LASSO fundida, as redes extraídas para os quatro conjuntos de dados diferentes com base em nosso método são exatamente as mesmas (exceto por pequenas diferenças nos pesos - ver Fig. Suplementar S1), que é biologicamente esperado. No entanto, este não foi o caso para as redes reconstruídas por ARACNE, CLR, GENIE3, catnet, a deconvolução de rede e os métodos de silenciamento global, além disso, as redes de GENIE3, catnet e o método de deconvolução foram as mais densas em todos os conjuntos de dados. Ao todo, os modelos baseados em regularização subestimaram amplamente as verdadeiras bordas positivas, embora fossem consistentes na previsão das bordas para cada conjunto de dados único. A abordagem proposta teve um desempenho particularmente bom no que diz respeito à previsão do tipo e não apenas da presença de uma relação regulatória, evidente nos casos dos fatores sigma RpoH e RpoD.

Sub-redes incluindo fatores sigma.

A rede reguladora de genes para quatro fatores sigma, RopD, RopE, RopH e RopS, juntamente com suas interações verificadas experimentalmente obtidas no RegulonDB 53. As bordas coloridas pertencem à sub-rede recuperada do RegulonDB, onde bordas vermelhas denotam ativação, enquanto bordas azuis indicam relações regulatórias reprimidas. As bordas marcadas em verde são de tipo regulatório não especificado. Se uma aresta foi prevista por um método, mas não está incluída na rede do RegulonDB, ela será colorida em preto. As bordas previstas para ARANCE, catnet, deconvolução de rede e métodos de silenciamento global são marcadas por ‘A’, ‘C’, ‘D’ e ‘S’, respectivamente, próximo às bordas correspondentes. As letras são codificadas por cores - fontes vermelhas, azuis ou verdes representam relacionamentos ativadores, reprimidos ou não especificados, respectivamente. As bordas pontilhadas denotam os relacionamentos previstos pela abordagem proposta. Ilustrados são os relacionamentos regulatórios previstos e seus tipos com base em dados de (UMA) frio, (B) aquecer, (C) estresse oxidativo e (D) experimentos de série temporal com deslocamento lactose-diaux.

Por fim, obtivemos o tempo de execução de cada método de inferência que é apresentado na Tabela Suplementar S2. Como o método proposto inclui muitos modelos de regressão independentes, paralelizamos a abordagem (conforme descrito em Métodos). Portanto, o tempo de execução do método proposto é fornecido apenas para realizar uma única regressão. Claramente, os algoritmos que contam com inversões de matriz na presença de parâmetros ocultos (por exemplo, deconvolução de rede) são mais rápidos em comparação com o método proposto, que requer a resolução de regressões múltiplas. Embora se espere que a validação cruzada aumente a demanda computacional, se um único valor for usado em todos os genes (ou seja, modelos) a interação regulatória na vizinhança de um gene pode ser inferida em um tempo drasticamente reduzido da ordem de alguns segundos.

Mycobacterium tuberculosis

Motivado pelas previsões da aplicação da abordagem proposta na combinação de todos E. coli conjuntos de dados, em seguida, investigamos a análise comparativa apenas na combinação de ambos os conjuntos de dados de MTB (ou seja, hipóxia e reaeração). Para obter redes regulatórias de genes, aplicamos todos os métodos comparados aos valores de expressão gênica pré-processados, bem como as mudanças de dobra da expressão gênica transformada em log 2 entre o controle (ponto de tempo 0) e as amostras de séries temporais de hipóxia e reaeração. A análise ROC para os métodos comparados foi então obtida usando o pacote R mineto 58 para ambos: os 31 primeiros (como a rede padrão ouro inclui 31 arestas) e 100 arestas altamente classificadas e as estatísticas correspondentes estão resumidas na Tabela Suplementar S3. É evidente que a abordagem proposta permitiu a previsão de relações regulatórias em uma porcentagem maior de verdadeiros positivos em ambos os casos de uso de valores de expressão gênica e alterações de dobra de expressão gênica transformada em log 2 quando comparados ao desempenho de GENIE3, deconvolução de rede e silenciamento global métodos. Da mesma forma com E. coli, GGM, ARACNE, CLR e modelos baseados em regularização tiveram um desempenho fraco para inferência precisa da rede regulatória, no entanto, o desempenho da maioria dos métodos comparados foi aumentado após a aplicação às alterações de dobra de expressão gênica transformada em log 2. Além disso, os valores do seletor para todos os métodos comparados foram iguais e iguais a um.

Furthermore, to specifically inspect the efficiency of the inferred networks, we investigated on three experimentally verified TF regulatory interactions that are highlighted in the Result section from Galagan et al. (page 180, 55 ): TF Rv0081 negatively regulates TFs Rv3597c and Rv3416 (whiB3) (Rv0081 → Rv3597c (Lsr2), Rv0081 → Rv3416 (whiB3)) and TF Rv3133c and Rv2034 (DosR) negatively regulate each other (Rv3133c ↔ Rv2034 (DosR)). All interactions and their regulatory types were successfully predicted by the proposed approach with sufficiently high ranks (Supplementary Table S4). Global silencing as well as network deconvolution methods were not successful with respect to the type of regulatory interactions and obtained rather low ranks, while GENIE3 predicted all interactions with high ranks but is not able to infer the type of regulatory relationships. Likewise, GGM, ARACNE and CLR were not successful in predicting the interactions, while regularization-based models performed very inconsistent with respect to the sign and predicting true edges.

Moreover, we counted the number of predicted true interactions in the top 100 highly ranked edges (Supplementary Tables S4, S5 and S6) considering the gold standard sub-network which includes 31 interactions (Fig. 2 in 55 ). It is evident that the proposed approach obtained the highest overlap with the gold standard while the minimum rank of the intersected interactions is remarkably high (above 0.6). Finally, the networks predicted by the compared approaches were filtered by the corresponding median edge ranks and the ROC analysis has been performed to the resulting sub-networks. Here, too, our proposed approach resulted in the highest AUROC and AUPR values while the edges in the filtered sub-network have considerably higher ranks (above 0.6075 and 0.379 when applied to gene expression levels and log 2-transformed gene expression fold changes, respectively) in comparison to the median ranks of the other contending methods.

Mus musculus

Finally, to verify the performance of the proposed approach in a higher organism, we performed a comparative analysis on the combination of tissue-specific time-series data sets from three genetically distinct mouse ES cell lines during differentiation. We applied all compared methods to the log 2-transformed gene expression fold changes between control (time point 0) and the rest of time-series from the corresponding ES cell lines, motivated by the resulting improvement in the performance of the inference approaches applied to the MTB data sets. The ROC analysis for the compared methods was then carried out on the top 248 (as the gold standard network includes 248 edges) and 500 highly ranked edges and the corresponding statistics are summarized in Supplementary Table S7.

It is evident that the proposed approach allowed the prediction of regulatory relationships at a higher percentage of true positives when compared to the performance of GENIE3, network deconvolution and global silencing methods. Likewise with other data sets CLR as well as regularization models performed poorly for accurate inference of the regulatory network however, the performance of the GGM is increased in comparison to its performance on the other data sets, while ARACNE failed to predict any true positive edge. In addition, the selector values for most of the compared methods were high and close (or equal) to one.

Moreover, we counted the number of predicted true interactions in the top 248 highly ranked edges (Supplementary Tables S7 and S9) considering the gold standard sub-network (which includes 248 interactions). It is evident that the proposed approach obtained the highest overlap with the gold standard, while the minimum rank of the intersected interactions is remarkably high (above 0.5). Finally, the networks predicted by the compared approaches were filtered by the corresponding median edge ranks and the ROC analysis has been performed to the resulting sub-networks. Once again the proposed approach resulted in the highest AUROC and AUPR values, while the edges in the filtered sub-network have similar ranks as the other predicted sub-networks by the other contending methods.


Discussão

TRNs provide an important and popular framework for better understanding a cell’s regulatory mechanisms, leading to phenotypic conditions. However, to the best of our knowledge TRN reconstruction methods today do not incorporate phenotypic information adequately or at all. As such, the reconstructed networks may be limited in pinpointing regulatory mechanisms most related to a phenotype under investigation, and often necessitate a follow-up step that filters for phenotype relevance. For example, a recent study of gene expression changes underlying Huntington’s disease (HD) 73 reconstructed a TRN specific to the mouse striatum and then short-listed TFs whose predicted targets were enriched in genes differentially expressed in HD mouse models. In another study, gene expression profiles of TFs and putative target genes were used to reconstruct a context-restricted TRN for breast cancer (using only breast cancer samples), and then a list of breast cancer-relevant TFs (called “risk-TFs”) whose regulons were enriched in risk loci were short-listed 57 . In the aforementioned study 57 , GWAS and eQTL analyses were used to define risk loci and relate them to the regulon of each TF. Such previous attempts to augment TRN reconstruction with phenotypic data motivated us to develop a systematic approach to incorporate information about the phenotype directly into TRN reconstruction.

In this study, we developed InPheRNo to reconstruct phenotype-relevant TRNs and utilized it to identify regulatory interactions that differentiate one cancer type from others while correcting for the confounding effect of tissues of origin. InPheRNo is based on a carefully designed PGM, which is key to combining TF–gene expression correlations with gene–phenotype associations. The conditional distributions of the PGM model the summary statistics of gene–phenotype and TF–gene associations, providing a succinct and efficient approach for data integration to identify phenotype-relevant regulatory relationships. The method is broadly applicable since it learns regulatory relationships from expression data alone and does not impose any restriction on the type of phenotype under investigation—the phenotype may be binary, categorical or even continuous-valued, and any appropriate statistical method for testing its association with a gene’s expression may be used in InPheRNo. Unlike several other methods that rely on the regulatory relationship of one TF–gene pair at a time, InPheRNo considers the effect of multiple TFs on each gene in the reconstruction procedure, at the time of selecting candidate TFs as well as in training the PGM. Finally, using posterior probabilities obtained from the PGM, InPheRNo provides a score representing the confidence for the identified phenotype-relevant regulatory edges.

In designing InPheRNo’s pipeline, we made the choice to first perform a feature selection step (using Elastic Net) and only use the selected TFs in the PGM. First and foremost, this was done to reduce the computational complexity, both by reducing the number of candidate TFs and also by summarizing the expression profiles of genes and TFs using summary statistics. Several previous studies have successfully used summary statistics (and particularly p values) for similar reasons 23,28,29,30,31,32 . Second, modeling summary statistics instead of the full gene expression data enables integration of other regulatory evidence (captured through data types other than transcriptomic, if available) in the PGM with a relative ease.

One important consideration when using InPheRNo, is the number of samples. As InPheRNo is based on modeling of summary statistics obtained from gene–phenotype and gene–TF associations, similar requirements on the minimum number of samples for those analyses should be also considered here 74,75,76 . However, two features of InPheRNo enable it to handle a small number of samples better than traditional co-expression analysis. First, it utilizes Elastic Net (as part of the pipeline), whose regularization terms can overcome some limitations of the small sample size by imposing sparsity criterion. Second, as its PGM models the distribution of the p values instead of relying on whether such p value are significant or not (i.e., instead of thresholding them) it is more robust towards small samples sizes.

As there are no rigorously validated metazoan TRNs to benchmark against, we evaluated the predicted TRNs indirectly through key TFs and gene expression signatures derived from them, and showed clear improvement over several related strategies. Our results showed that the TFs with many cancer type-relevant targets are potential cancer driver TFs and may suggest novel drug targets or provide new insights, regarding the development and progress of cancer. Our results also suggest a powerful approach for subtyping of cancer patients using gene expression signatures: while most approaches developed for this task do not take into account the regulatory interactions among genes, our survival analysis suggests that cancer type-relevant TRNs can improve the predicting power of gene expression signatures.

In spite of the success of the InPheRNo-based gene signatures in differentiating between patients with poor and good prognosis for the majority of cancer types, in some cases, e.g., BRCA, this method did not result in groups with significantly different survival probability, despite the existence of BRCA-driver TFs in the signature. This lack of success may partially be owing to the fact that we clustered samples of each cancer type into two clusters, whereas these cancer types may include more than two subtypes, as is the case in BRCA 26 . However, since in most cancer types a definite number for the cancer subtypes is not yet established, we preferred to keep the number of clusters equal to two. A more in-depth analysis of subtype discovery and survival analysis using InPheRNo-derived TRNs is left for future work.

We would like to emphasize that in this study, we focused only on transcriptomic data, owing to the availability of this data type in many domains, including domains outside of cancer research, and lack of other important data types such as ChIP-seq data in these domains. Even in the area of cancer research, in which large databases of ChIP-seq tracks (such as ENCODE) corresponding to various cancer cell lines are available, the datasets are extremely biased toward a small fraction of well-studied TFs (for example only

10% of all TFs are studied in ENCODE). As a result, including these data sets may significantly bias the analysis towards this small fraction of TFs. In addition, matched gene expression and ChIP-seq data for tumor samples are rarely available and combining these data types from different sources and different samples, in itself a significant challenge, will require substantial effort in the future.

We believe that including additional types of regulatory evidence (especially those representing “cis” mechanisms such as TF motifs and chromatin state changes) in the phenotype-relevant TRN reconstruction procedure is an important and essential future direction for improving InPheRNo. This is especially true considering that many efforts are under way to generate large datasets containing matching transcriptomic, genomic, epigenomic and phenotypic profiles of patients 77,78,79 . One way to achieve this goal might be to include different regulatory evidence as new observed variables in the PGM used in InPheRNo. Another alternative is to use cis-regulatory evidences to construct an initial network that is used as a “prior” for Bayesian analysis of expression data, as has been demonstrated before 80 . Future investigations should focus on these avenues of integrating multi-omics data into the InPheRNo model.


Welcome!

This is one of over 2,400 courses on OCW. Explore materials for this course in the pages linked along the left.

MIT OpenCourseWare is a free & open publication of material from thousands of MIT courses, covering the entire MIT curriculum.

No enrollment or registration. Freely browse and use OCW materials at your own pace. There's no signup, and no start or end dates.

Knowledge is your reward. Use OCW to guide your own life-long learning, or to teach others. We don't offer credit or certification for using OCW.

Made for sharing. Download files for later. Send to friends and colleagues. Modify, remix, and reuse (just remember to cite OCW as the source.)


Mapping gene regulatory networks from single-cell omics data

Single-cell techniques are advancing rapidly and are yielding unprecedented insight into cellular heterogeneity. Mapping the gene regulatory networks (GRNs) underlying cell states provides attractive opportunities to mechanistically understand this heterogeneity. In this review, we discuss recently emerging methods to map GRNs from single-cell transcriptomics data, tackling the challenge of increased noise levels and data sparsity compared with bulk data, alongside increasing data volumes. Next, we discuss how new techniques for single-cell epigenomics, such as single-cell ATAC-seq and single-cell DNA methylation profiling, can be used to decipher gene regulatory programmes. We finally look forward to the application of single-cell multi-omics and perturbation techniques that will likely play important roles for GRN inference in the future.

Bonecos

Single-cell GRNs. The goal of…

Single-cell GRNs. The goal of many single-cell studies is to understand which cell…


6. How Many Gene Regulatory Networks Exist?

It is generally acknowledged that a phenotype is an emergent property of genotype-environment interactions. Specifically, a phenotype results from molecular and cellular activity patterns from genotype-environment interactions. This implies that each observable phenotype is associated with phenotype-specific gene networks, because without changing molecular interactions a phenotype cannot change this concept is illustrated in Figure 1. In this figure, gene networks can be seen as a bottleneck between the genotype and the phenotype with respect to their coupling. That means every change on the genotype level that will result in a change of the phenotype will also inevitably lead to a change in the gene network structure as mediator between both levels.

Figure 1. Schematic overview of the general role gene networks play in understanding phenotypes.

However, since gene networks refer to all possible types of molecular networks, including the transcriptional regulatory network, protein interaction network, metabolic network, gene regulatory network and interactions between these networks, it is less clear which of these networks, or all of them, are actually changed. Moreover, because a gene regulatory network can potentially represent many types of physical biochemical interactions among genes and gene products (de Matos Simoes et al., 2013a) it can be expected that gene regulatory networks are highly phenotype specific (Schadt, 2009 Emmert-Streib and Glazko, 2011). Establishing such relationships will therefore be a complex task, but also provides an opportunity to catalog phenotypes quantitatively. An example for the analysis of tissue-specific networks can be found in Guan et al. (2012) where 107 tissue specific network have been studied. Currently, the number of GRNs is difficult to estimate but based on these preliminary results one can hypothesize that there are more than 200 different GRNs for Human alone, because this corresponds about to the number of different cell types. However, also pathological cells manifesting tumors have their own characteristic networks (Emmert-Streib et al., 2014) implying that there are probably thousands of different gene networks in Human.


RMaNI: Regulatory Module Network Inference framework

Fundo: Cell survival and development are orchestrated by complex interlocking programs of gene activation and repression. Understanding how this gene regulatory network (GRN) functions in normal states, and is altered in cancers subtypes, offers fundamental insight into oncogenesis and disease progression, and holds great promise for guiding clinical decisions. Inferring a GRN from empirical microarray gene expression data is a challenging task in cancer systems biology. In recent years, module-based approaches for GRN inference have been proposed to address this challenge. Despite the demonstrated success of module-based approaches in uncovering biologically meaningful regulatory interactions, their application remains limited a single condition, without supporting the comparison of multiple disease subtypes/conditions. Also, their use remains unnecessarily restricted to computational biologists, as accurate inference of modules and their regulators requires integration of diverse tools and heterogeneous data sources, which in turn requires scripting skills, data infrastructure and powerful computational facilities. New analytical frameworks are required to make module-based GRN inference approach more generally useful to the research community.

Resultados: We present the RMaNI (Regulatory Module Network Inference) framework, which supports cancer subtype-specific or condition specific GRN inference and differential network analysis. It combines both transcriptomic as well as genomic data sources, and integrates heterogeneous knowledge resources and a set of complementary bioinformatic methods for automated inference of modules, their condition specific regulators and facilitates downstream network analyses and data visualization. To demonstrate its utility, we applied RMaNI to a hepatocellular microarray data containing normal and three disease conditions. We demonstrate that how RMaNI can be employed to understand the genetic architecture underlying three disease conditions. RMaNI is freely available at http://inspect.braembl.org.au/bi/inspect/rmani

Conclusão: RMaNI makes available a workflow with comprehensive set of tools that would otherwise be challenging for non-expert users to install and apply. The framework presented in this paper is flexible and can be easily extended to analyse any dataset with multiple disease conditions.


4. Results and discussion

4.1. Simulated BNps with 7 genes

We first evaluate different inference algorithms on synthetically generated random networks. We generate 1000 random BNps with n = 7 genes, maximum input degree K = 3, and perturbation probability p = 0.01. For each node, we uniformly assign 1 to K regulators. Hence the average connectivity in this set of random networks is 2. After determining the regulatory relationships among nodes, the regulatory functions for each node are determined by randomly filling in the corresponding truth tables with Bernoulli random numbers with the bias following a Beta distribution with mean 0.5 and standard deviation 0.01. For each random BNp, we simulate time series of different numbers of state transitions based on its underlying Markov chain. The number of “observed” state transitions m ranges from 10 to 60 to reflect the difficulty level of network inference. For control, we choose the first node as the marker gene and define the undesirable states as these network states with the first node down-regulated. In the binary representation of network states,

= <x|x1 = 0>. As the networks are randomly generated, without loss of generality, we allow intervention on the last node as the control gene, which we can either knock up or down to derive control policies. In our simulated random BNps, we have the original average undesirable steady state mass π org = 0.5071 with standard deviation 0.3575, with π org ≈ 0.5 because we set the bias to 0.5. When we apply the MSSA algorithm to derive the optimal stationary control policies for these random BNps, the average controlled undesirable steady state mass is π = 0.3703 with the standard deviation 0.3749.

Based on these simulated time series, we have implemented REVEAL, BIC, MDL, uMDL, and Best-Fit inference algorithms and modified accordingly to reconstruct BNps, including regulatory relationships and regulatory functions represented as general truth tables. For BIC and MDL, we set the regularization coefficients to values previously reported to have good performance in Zhao et al. (2006), λ = 0.5 for BIC and λ = 0.3 for MDL.

Table ​ Table1 1 provides the network inferential validity measurements: normalized Hamming distance μham (Hamming distance over the total number of edges in true networks), the steady-state distance μss, and the controllability distance μctrl for different network inference algorithms given different numbers of state transitions. As discussed in (Zhao et al., 2006), BIC and MDL perform similarly. Regarding the accurate recovery of regulatory relationships, it is interesting to see that Best-Fit appears to achieve the best performance with respect to μham while REVEAL does not perform very well. One explanation could be that REVEAL introduces many false positives, hopefully to best fit the data by using the functions with more regulators. This is in fact what we observe from our experiments. All the other inference algorithms choose the functions with the smallest number of regulators either by complexity regularization in BIC, MDL, and uMDL or choosing the “parsimonious” functions with the minimum prediction errors in Best-Fit. For uMDL, we note that μham improves quickly with the increasing sample size compared to other complexity regularization algorithms BIC and MDL. Based on our experiments, uMDL consistently generates very low false positive edges (close to zero), even with a very limited number of samples, which is the main advantage of the uMDL algorithms. This has also been shown in the original paper (Dougherty et al., 2008). For μss, both REVEAL and Best-Fit perform consistently better than BIC, MDL, and uMDL, since both REVEAL and Best-Fit aim to find the network models that best fit the observed state transitions. With regularization on model complexity by BIC, MDL, and uMDL, the steady-state distances are greater. As mentioned earlier, REVEAL and Best-Fit, especially REVEAL, reconstruct networks with more edges to explain the observed data, which leads to smaller μss.

Tabela 1

The comparison of network inference algorithms (REVEAL, BIC, MDL, uMDL, and Best-Fit) with M different number of observed state transitions.

Validityμhamμssμctrl
M103050103050103050
REVEAL0.77740.61110.65110.67430.46570.42160.10670.02750.0049
BIC0.69660.41960.33040.86790.70890.54920.07390.03000.0126
MDL0.72040.42600.32940.94140.72250.54350.07750.03110.0121
uMDL0.80000.37280.24711.19570.69730.49350.10580.03520.0093
Best-Fit0.73110.39190.29130.63780.42440.40980.10270.02500.0045

When we investigate the inferential validity with respect to controllability, μctrl, we see interesting changes of tendency between the five algorithms. Especially with very few state transitions, M = 10, BIC, MDL, and uMDL algorithms perform better than REVEAL and Best-Fit, which indicates that the regularization on model complexity with a limited number of observations helps reconstruct network models that yield better controllers. With more observations, REVEAL and Best-Fit gradually perform better than BIC, MDL, and uMDL due to introduced bias by model complexity regularization.

Figure ​ Figure2 2 plots μham, μss, and the average undesirable steady-state mass using the control policy designed on the inferred network via the MSSA algorithm. For comparison purposes, the latter average is compared to the average original undesirable mass and the average undesirable mass following application of the MMSA control policy designed on the original network. As m increases from 10 to 60, all algorithms improve. In fact, with more than 50 observed state transitions for these generated random BNps, the derived stationary control policies achieve almost the same performance compared to the optimal control policies with complete knowledge of the network models. The average performances from inferred networks are in fact within 5% for all five inference algorithms when M = 60.

Performance comparison of five network inference algorithms by different validity indices based on simulated BNps with 7 genes and K = 3. (A) Average normalized Hamming distance μham (B) μss (C) average undesirable steady-state mass π after applying derived stationary control policies based on inferred networks to the original ground truth BNps, compared to the average undesirable mass obtained by the optimal control policy (OPT) based on the complete knowledge of original BNps and the average undesirable mass before intervention (Original).

We further evaluate inference algorithms on a similar set of 1000 random BNps with n = 7 genes with the same settings but change the maximum input degree K = 5, which increases the average connectivity to 3. For this set of random BNps, we have the average undesirable original steady state mass π org /> = 0.4841 with standard deviation 0.3171. When we apply the MSSA algorithm to derive the optimal stationary control policies for these random BNps, the average controlled undesirable steady state mass is π /> = 0.2529 with the standard deviation 0.3144. The average shift of undesirable masses is higher compared to the previous set of random networks, which is expected as the network sensitivity monotonically increases with the average network connectivity (Kauffman, 1993 Shmulevich and Dougherty, 2007 Qian and Dougherty, 2009a). With higher sensitivity, networks can be more effectively controlled. We again compare the inferential validity as in the previous experiment. Figure ​ Figure3 3 shows plots analogous to Figure ​ Figure2. 2 . Especially, we note that in this set of experiments, we can achieve close-to-optimal intervention with fairly small sample size as illustrated in Figure ​ Figure3C. 3C . It is clear that the performance of different inference algorithms depends on the characteristics of the networks, especially the network sensitivity. More specifically, all three indices become worse for all the inference algorithms, illustrating that with increasing network sensitivity, the inference problem becomes more difficult. It is also clear that the performance improves at slower rates with the increasing sample size when we have higher network sensitivity. Another important difference is that for this set of random networks, both REVEAL and Best-Fit have higher μham when the number of samples increase above 40. The reason may be due to the tendency of random perturbations forcing both algorithms to bias toward more complex Boolean functions with more input variables as regulators.

Performance comparison of five network inference algorithms by different average validity indices based on BNps with 7 genes and K = 5. (A) Average normalized Hamming distance μham (B) μss (C) average undesirable steady-state mass π after applying derived stationary control policies based on inferred networks to the original ground truth BNps, compared to the average undesirable mass obtained by the optimal control policy (OPT) based on the complete knowledge of original BNps and the average undesirable mass before intervention (Original).

4.2. Simulated BNps with 9 genes

For simulations with 9 genes, owing to run time, we generate 200 BNps with n = 9 genes and perturbation probability p = 0.01. We again make uniformly random assignments of 1 to K regulators, with K = 3 so that the average connectivity is 2. The bias for the corresponding truth tables follows the same Beta distribution with mean 0.5 and stand deviation 0.01. The number of “observed” state transitions m range from 10 to 60. The derivation of control policies is still based on the definition of the undesirable states

= <x|x1 = 0> and the last node is the control gene. In the simulated random BNps, the average undesirable steady state mass is π org /> = 0.4886 with the standard deviation 0.3764. When we apply the MMSA algorithm to derive the optimal stationary control policies for these random BNps, the average controlled undesirable steady state mass is π /> = 0.3668 with the standard deviation 0.3863. Figure ​ Figure4 4 shows plots analogous to Figure ​ Figure2 2 with the trends similar as those observed in the previous experiments with corresponding random BNps with 7 genes and K = 3.

Performance comparison of five network inference algorithms by different average validity indices based on BNps with 9 genes and K = 3. (UMA) Average normalized Hamming distance μham (B) μss (C) average undesirable steady-state mass π after applying derived stationary control policies based on inferred networks to the original ground truth BNps, compared to the average undesirable mass obtained by the optimal control policy (OPT) based on the complete knowledge of original BNps and the average undesirable mass before intervention (Original).

In the second set of simulated random BNps with 9 genes, the settings are the same except that K = 5. In these random networks, the average undesirable steady state mass is π org /> = 0.4895 with standard deviation 0.3269. When we apply the MSSA algorithm to derive the optimal stationary control policies for these random BNps, the average controlled undesirable steady state mass is π /> = 0.2781 with standard deviation 0.3268. Figure ​ Figure5 5 is analogous to Figure ​ Figure3 3 .

Performance comparison of five network inference algorithms by different average validity indices based on BNps with 9 genes and K = 5. (UMA) Average normalized Hamming distance μham (B) μss (C) average undesirable steady-state mass π after applying derived stationary control policies based on inferred networks to the original ground truth BNps, compared to the average undesirable mass obtained by the optimal control policy (OPT) based on the complete knowledge of original BNps and the average undesirable mass before intervention (Original).

In summary, when we evaluate different inference procedures with respect to different inferential validity criteria, different inference procedures show different trends with their increasing sample size. Their performance overall depends on network characteristics as well as available samples. Finally, when effective intervention is our final operational objectivel, it is promising that we can achieve effective intervention based on inferred networks, even with fairly small sample size as illustrated in Figures ​ Figures2C, 2C , ​ ,3C, 3C , ​ ,4C, 4C , ​ ,5C 5C .

4.3. A metastatic melanoma network

Finally, we evaluate different inference algorithms based on a metastatic melanoma network used in previous studies on network intervention (Qian and Dougherty, 2008 Qian et al., 2009 Yousefi and Dougherty, 2013). The network has 10 genes listed in the order from the most to the least significant bit: WNT5A, PIR, S100P, RET1, MMP3, PLCG1, MART1, HADHB, SNCA, and STC2. The order does not affect our analysis. We note here that this network was derived from gene expression data (Kim et al., 2002) collected in studies of metastatic melanoma (Bittner et al., 2000 Weeraratna et al., 2002). Table ​ Table2 2 and Figure ​ Figure6 6 together illustrate the regulatory relationships among these selected 10 genes from 587 genes profiled in Bittner et al. (2000), Weeraratna et al. (2002), which were derived based on gene expression data rather than curated regulatory relationships among genes in literature. We believe that the model is appropriate for the purpose of illustrating the effectiveness of objective inferential validity on quantifying the performance of inference procedures in this work. Based on these information, we construct a BNp with the perturbation probability p = 0.01. As in the previous studies, the control objective is based on the fact that up-regulation of WNT5A is associated with increased metastasis. Assim,

= <x|x1 = 1>. For this network, the undesirable steady-state mass is π = 0.2073 in the original network, which can be reduced as illustrated in Table ​ Table3 3 with different genes as potential targets using the MSSA algorithm on the original network. Based on this model, we simulate 20, 60, and 80 state transitions and infer the network based on these time series data using all five algorithms. As the primary objective here is to reduce the undesirable steady-state mass with WNT5A up-regulated, we focus on its shift derived by the MSSA algorithm based on the inferred networks using different inference algorithms.

Mesa 2

Regulatory functions in the metastatic melanoma network [Modified from Table ​ Table1 1 in Yousefi and Dougherty (2013)].


Assista o vídeo: Nuevos métodos de dopaje - El Dopaje Genético (Janeiro 2022).