Em formação

Valor de corte E para pesquisa RNASeq


Estou tentando restringir os possíveis critérios de pesquisa de um experimento RNASeq. Eu já tenho as sequências de cds, então as peguei e usei o programa Blast + para explodir os dados contra si mesmo para encontrar sequências que eram iguais.

Pelo que li, o valor E é baseado no tamanho do banco de dados. Meu banco de dados tem apenas cerca de 27.000 sequências. Isso significa que apenas valores E excepcionalmente baixos serão relevantes para mim? Há algum outro critério que devo pesquisar antes de escrever o programa para eliminar sequências semelhantes?


Indo por um valor e ou corte de bitscore fornecerá possíveis homólogos, mas parece que você deseja remover sequências redundantes. Se você deseja agrupar e combinar sequências semelhantes para fazer um banco de dados menor, você pode simplesmente ir por meio de identidade de sequência usando algo como CD-HIT. Isso, por exemplo, é o que é feito para produzir o conjunto UniRef do banco de dados UniProt.


Quais métodos estão disponíveis para encontrar um valor de corte para genes não expressos em RNA-seq?

Eu tenho uma matriz de contagem de expressão gênica produzida a partir de dados de seq de RNA em massa. Eu gostaria de encontrar genes que eram não expressos em um grupo de amostras e foram expressos em outro grupo.

O problema, claro, é que nem todos efetivamente genes não expressos terão 0 contagens devido a erros de sequenciamento ou porque foram expressos em um pequeno subconjunto de células.

Estou interessado em soluções usando R.


Fundo

A tecnologia de sequenciamento de alto rendimento está se tornando rapidamente o método padrão para medir os níveis de expressão de RNA (também conhecido como RNA-seq) [1]. O advento de tecnologias de sequenciamento rápido, juntamente com custos reduzidos, permitiu o perfil detalhado dos níveis de expressão gênica, impactando quase todos os campos das ciências biológicas e agora está sendo adotado para uso clínico [2]. A tecnologia de RNA-seq permite a identificação detalhada de isoformas de genes, eventos de translocação, variações de nucleotídeos e modificações de bases pós-transcricionais [3]. Um dos principais objetivos desses experimentos é identificar os genes diferencialmente expressos em duas ou mais condições. Esses genes são selecionados com base em uma combinação de limiar de mudança de expressão e corte de pontuação, que geralmente são baseados em P valores gerados por modelagem estatística.

O nível de expressão de cada unidade de RNA é medido pelo número de fragmentos sequenciados que mapeiam para o transcrito, que se espera correlacionar diretamente com seu nível de abundância. Essa medida é fundamentalmente diferente dos métodos baseados em sondas genéticas, como microarrays. No RNA-seq, o sinal de expressão de um transcrito é limitado pela profundidade do sequenciamento e depende dos níveis de expressão de outros transcritos, enquanto nos métodos baseados em array as intensidades das sondas são independentes umas das outras. Isso, assim como outras diferenças técnicas, tem motivado o desenvolvimento de um número crescente de algoritmos estatísticos que implementam uma variedade de abordagens para normalização e detecção de expressão diferencial (DE). Abordagens típicas usam Poisson ou distribuições binomiais negativas para modelar os dados de contagem de genes e uma variedade de procedimentos de normalização (ver [4] para uma revisão).

Neste estudo de comparação, avaliamos alguns dos pacotes de software de expressão diferencial mais comumente usados ​​e disponíveis gratuitamente: Cuffdiff [5], edgeR [6], DESeq [7], PoissonSeq [8], baySeq [9] e limma [ 10] adaptado para uso de RNA-seq. Usamos dois conjuntos de dados de referência: o primeiro é o conjunto de dados Sequencing Quality Control (SEQC), que inclui amostras replicadas do RNA de referência de corpo inteiro humano e RNA de referência do cérebro humano, juntamente com controles de pico de RNA. Essas amostras são parte do estudo MAQC para benchmarking de tecnologia de microarray [11, 12], bem como o esforço SEQC para caracterizar a tecnologia RNA-seq e incluem cerca de 1.000 genes que foram validados por TaqMan qPCR. O segundo conjunto de dados são dados de RNA-seq de réplicas biológicas de três linhas de células que foram caracterizadas como parte do projeto ENCODE [13]. Nossa análise se concentrou em uma série de medidas que são mais relevantes para a detecção de expressão diferencial de genes a partir de dados de RNA-seq: i) normalização de dados de contagem ii) sensibilidade e especificidade de detecção de DE iii) desempenho no subconjunto de genes que são expressos em uma condição, mas não têm expressão detectável na outra condição e, finalmente, iv) os efeitos da redução da profundidade de sequenciamento e do número de repetições na detecção da expressão diferencial. É importante ressaltar que esta avaliação não aborda o problema relacionado e importante de detecção de expressão diferencial de isoformas e identificação de novos transcritos. Em vez disso, a avaliação é restrita ao caso específico de detecção de DE com base em modelos de genes unificados.

Nossos resultados demonstram diferenças substanciais entre os métodos, tanto em termos de especificidade quanto de sensibilidade para a detecção de genes expressos diferencialmente. Na maioria dos benchmarks, Cuffdiff teve um desempenho menos favorável com um maior número de falsos positivos sem qualquer aumento na sensibilidade. Nossos resultados demonstram conclusivamente que a adição de amostras replicadas fornece um poder de detecção de DE substancialmente maior do que o aumento da profundidade da sequência. Portanto, a inclusão de mais amostras replicadas em experimentos de RNA-seq é sempre preferível em vez de aumentar o número de leituras sequenciadas.

Bases teóricas

Um ponto de partida conveniente para comparar diferentes métodos de análise de RNA-seq é uma matriz de contagem simples N do n × m Onde N eu j é o número de leituras atribuídas ao gene eu no experimento de sequenciamento j (isto é, contagens de leitura). Tais matrizes podem ser produzidas a partir de dados de alinhamento usando ferramentas como HTSeq [15], Picard [16], BEDTools [17], featureCounts [18] ou Cufflinks [19]. O estudo apresentado aqui não aborda as sutilezas importantes ao calcular contagens de genes, em particular qual modelo de gene usar, como contar leituras sobrepondo regiões intrônicas e o uso de leituras mapeadas de forma ambígua. Em vez disso, o foco está na comparação entre métodos, dada uma matriz de contagem de expressão fixa. Para Cuffdiff, que usa um método de quantificação diferente e não compatível com os outros, usamos o método conjunto Cufflinks e para todos os outros métodos usamos HTSeq. É importante reconhecer que o número de leituras que se sobrepõem a um gene eu não é uma medida direta da expressão do gene. Em vez disso, a medida de contagem onde µ eu j e eu eu são a expressão esperada e o comprimento do gene, respectivamente. Portanto, há um viés de comprimento claro ao medir a expressão gênica por RNA-seq [20]. Um efeito desse viés é reduzir a capacidade de detectar a expressão diferencial entre genes mais curtos simplesmente pela falta de cobertura, uma vez que o poder dos testes estatísticos envolvendo dados de contagem diminui com um número menor de contagens [21, 22].

A análise da expressão gênica diferencial de dados de RNA-seq geralmente consiste em três componentes: normalização de contagens, estimativa de parâmetro do modelo estatístico e testes de expressão diferencial. Nesta seção, fornecemos um breve histórico das abordagens implementadas pelos vários algoritmos que executam essas três etapas. Limitamos nossa discussão ao caso mais comum de medir a expressão diferencial entre duas condições celulares ou fenótipos, embora alguns dos pacotes possam testar diferenças multiclasse ou experimentos multifatoriais onde várias condições biológicas e diferentes protocolos de sequenciamento estão incluídos.

Normalização

A primeira dificuldade a abordar ao trabalhar com dados de sequenciamento são as grandes diferenças no número de leituras produzidas entre diferentes execuções de sequenciamento, bem como vieses técnicos introduzidos por protocolos de preparação de biblioteca, plataformas de sequenciamento e composições de nucleotídeos [23]. Os procedimentos de normalização tentam explicar essas diferenças para facilitar comparações precisas entre os grupos de amostra. Uma normalização intuitiva é dividir a contagem de genes simplesmente pelo número total de leituras em cada biblioteca, ou leituras mapeadas, conforme introduzido pela primeira vez por Mortazavi et al. [1], um procedimento de normalização denominado leituras por quilobase por milhão de leituras (RPKM). Uma deficiência dessa abordagem é que a representação proporcional de cada gene depende dos níveis de expressão de todos os outros genes. Freqüentemente, uma pequena fração de genes é responsável por grandes proporções das leituras sequenciadas e pequenas mudanças de expressão nesses genes altamente expressos distorcem a contagem de genes expressos humildemente sob esse esquema. Isso pode resultar em expressão diferencial errônea [24, 25]. Uma variação de RPKM, denominada fragmentos por quilobase de exon por milhão de leituras mapeadas (FPKM), foi introduzida por Trapnell et al. para acomodar leituras emparelhadas [19], no entanto, isso tem a mesma limitação de alterações de acoplamento nos níveis de expressão entre todos os genes. DESeq calcula um fator de escala para uma determinada amostra, calculando a mediana da razão, para cada gene, de sua contagem de leitura sobre sua média geométrica em todas as amostras. Em seguida, ele usa a suposição de que a maioria dos genes não é DE e usa essa mediana de razões para obter o fator de escala associado a essa amostra. Cuffdiff estende isso realizando primeiro o escalonamento da biblioteca intra-condição e, em seguida, um segundo escalonamento entre as condições. Cuffdiff também tenta explicar as mudanças nos níveis de isoforma explicitamente por normalização específica de transcrição adicional que estima a abundância de cada isoforma.

Outros procedimentos de normalização tentam usar um subconjunto de genes expressos de forma estável ou normalizar dentro de amostras replicadas para ajustar globalmente os tamanhos da biblioteca. As médias aparadas dos valores M (TMM) de Robinson e Oshlack [25], que é implementado em edgeR, calcula um fator de escala entre dois experimentos usando a média ponderada do subconjunto de genes após a exclusão de genes que exibem contagens de leitura média alta e genes que têm grandes diferenças na expressão. Outra abordagem é somar contagens de genes até o quantil superior de 25% para normalizar os tamanhos das bibliotecas, conforme proposto por Bullard et al. [24] e é a normalização padrão no pacote baySeq. O pacote PoissonSeq usa uma estimativa de adequação para definir um conjunto de genes menos diferenciado entre duas condições, que é então usado para calcular os fatores de normalização da biblioteca. A normalização de quantis garante que as contagens em todas as amostras tenham a mesma distribuição empírica, classificando as contagens de cada amostra e definindo os valores para serem iguais à média dos quantis de todas as amostras [26]. Essa normalização é amplamente usada em matrizes de expressão e é implementada no pacote limma. Recentemente, uma nova função de normalização denominada voom projetada especificamente para dados de RNA-seq foi adicionada ao pacote limma. Ele executa uma regressão LOWESS para estimar a relação média-variância e transforma as contagens de leitura para a forma de log apropriada para modelagem linear [27].

Modelagem estatística da expressão gênica

Se os experimentos de sequenciamento forem considerados como amostras aleatórias de leituras de um conjunto fixo de genes, então uma representação natural das contagens de leitura de genes é a distribuição de Poisson da forma em que n é o número de contagens de leitura e λ é um número real igual ao número esperado de leituras de fragmentos da transcrição. Uma propriedade importante da distribuição de Poisson é que a variância é igual à média, que é igual λ. No entanto, na realidade, a variância da expressão gênica em múltiplas réplicas biológicas é maior do que seus valores médios de expressão [28-30]. Para resolver este problema de superdispersão, métodos como edgeR e DESeq usam a distribuição binomial negativa relacionada (NB), onde a relação entre a variância ν e significa µ é definido como ν = µ + αμ 2 onde α é o fator de dispersão.

A estimativa desse fator é uma das diferenças fundamentais entre os pacotes edgeR e DESeq. estimativas edgeR α como uma combinação ponderada de dois componentes: um efeito de dispersão específico do gene e um efeito de dispersão comum calculado a partir de todos os genes. DESeq, por outro lado, divide a estimativa da variância em uma combinação da estimativa de Poisson (ou seja, a expressão média do gene) e um segundo termo que modela a variabilidade da expressão biológica. Cuffdiff calcula um modelo de variância separado para genes de isoforma única e genes de isoforma múltipla. A variância da expressão de isoforma única é calculada de forma semelhante ao DESeq e a variância de isoformas múltiplas é modelada por um modelo de mistura de binômios negativos usando os parâmetros de distribuição beta como pesos de mistura. baySeq implementa um modelo bayesiano completo de distribuições binomiais negativas em que os parâmetros de probabilidade anteriores são estimados por amostragem numérica dos dados. PoissonSeq modela as contagens de genes N eu, j como uma variável de Poisson em que a média µ eu, j da distribuição é representado pelo log de relacionamento log-linear µ eu j = log d j + log β eu + γ eu y j Onde d j representa o tamanho normalizado da biblioteca, β eu é o nível de expressão do gene eu e γ eu é a correlação do gene eu com condição y j (note que em [8] os subscritos eu e j são amostras e genes, respectivamente). Se a expressão do gene eu não está correlacionado com a amostra j classe (ou seja, não há diferença significativa no gene eu expressão entre duas condições), então γ eu é zero.

Teste de expressão diferencial

A estimativa dos parâmetros para o respectivo modelo estatístico é seguida do teste de expressão diferencial, o cálculo da significância da mudança na expressão do gene. eu entre duas condições. Tanto o edgeR quanto o DESeq usam uma variação do teste exato de Fisher adotado para a distribuição NB, portanto, eles retornam o exato P valores calculados a partir das probabilidades derivadas. Cuffdiff usa as estatísticas de teste T = E[registro(y)] / Var [log (y)], Onde y é a razão das contagens normalizadas entre duas condições, e essa razão segue aproximadamente uma distribuição normal, portanto, um teste t é usado para calcular o P valor para DE. limma usa uma estatística t moderada para calcular P valores nos quais tanto o erro padrão quanto os graus de liberdade são modificados [10]. O erro padrão é moderado entre os genes com um fator de encolhimento, que efetivamente pega emprestado informações de todos os genes para melhorar a inferência em qualquer gene único. Os graus de liberdade também são ajustados por um termo que representa o a priori número de graus de liberdade para o modelo. A abordagem baySeq estima dois modelos para cada gene, um assumindo nenhuma expressão diferencial e um segundo assumindo expressão diferencial usando os dois grupos de amostra. A probabilidade posterior do modelo de DE, dados os dados observados, é usada para identificar genes expressos diferencialmente. No método PoissonSeq, o teste de expressão diferencial é simplesmente um teste para a significância do γ eu termo (isto é, correlação de gene eu expressão com as duas condições), que é avaliada por estatísticas de pontuação. Por experimentos de simulação, foi mostrado que essas estatísticas de pontuação seguem uma distribuição qui-quadrada, que é usada para derivar P valores para DE. Todos os métodos usam abordagens padrão para correção de múltiplas hipóteses (por exemplo, Benjamini-Hochberg) com exceção de PoissonSeq, que implementou uma nova estimativa de taxa de descoberta falsa (FDR) para dados de contagem que é baseada em permutação.


Perfis transcriptômicos baseados em sequenciamento de RNA de desenvolvimento de lentes embrionárias para descoberta de genes de catarata

Cataratas congênitas isoladas ou sindrômicas são defeitos de desenvolvimento heterogêneos, tornando a identificação dos genes associados um desafio. No passado, microarrays de expressão de lente de camundongo foram aplicados com sucesso em ferramentas de bioinformática (por exemplo, iSyTE) para facilitar a descoberta de genes associados à catarata humana. Para desenvolver um novo recurso para geneticistas, relatamos perfis de sequenciamento de RNA de alto rendimento (RNA-seq) de lentes de camundongos em estágios embrionários principais (E) 10.5 (fosso da lente), E12.5 (diferenciação celular de fibra primária), E14.5 e E16.5 (diferenciação de células de fibra secundária). Esses estágios capturam eventos importantes à medida que a lente se desenvolve de um placódio invaginante para um tecido transparente. Anteriormente, a expressão "enriquecida em lente" baseada em subtração de corpo de embrião inteiro in silico (WB) foi eficaz na priorização de genes ligados à catarata. Para aplicar uma abordagem análoga, geramos novos conjuntos de dados WB RNA-seq de camundongo e mostramos que a subtração WB in silico de conjuntos de dados de RNA-seq de lente identifica com sucesso genes-chave com base na expressão enriquecida com lente. Em expressão de ≥2 contagens por milhão, ≥1,5 log2 limite de enriquecimento de vezes (p & lt 0,05), lente E10.5 exibe 1401 genes enriquecidos (17% de genes expressos em lente), lente E12.5 exibe 1937 genes enriquecidos (22% de genes expressos em lente), lente E14.5 exibe 2514 genes enriquecidos (31% dos genes expressos na lente) e a lente E16.5 exibe 2.745 genes enriquecidos (34% dos genes expressos na lente). A análise da via biológica identificou genes associados ao desenvolvimento de lentes, regulação da transcrição e vias de sinalização, entre outros grupos funcionais. Além disso, esses novos dados de RNA-seq confirmaram a alta expressão de genes ligados à catarata e identificaram novos reguladores potenciais no cristalino. Finalmente, desenvolvemos novas trilhas de anotação UCSC Genome Brower específicas do estágio de lente e as tornamos publicamente acessíveis através do iSyTE (https://research.bioinformatics.udel.edu/iSyTE/) para visualização amigável da expressão / enriquecimento do gene da lente para priorizar genes de dados de alto rendimento de casos de catarata.

Bonecos

Fig. 1 .. Desenho experimental de RNA-seq para camundongo ...

Fig. 1 .. Desenho experimental de RNA-seq para análise de transcriptoma de lente embrionária de camundongo.

Um fluxograma experimental ...

Fig. 2 .. Análise de qualidade de dados de RNA-seq ...

Fig. 2 .. Análise de qualidade de dados de RNA-seq e eficácia de em sílico Subtração WB.

Fig. 3 .. Em sílico Subtração WB e expressão ...

Fig. 3 .. Em sílico A subtração de WB e a análise de expressão identificam efetivamente os genes ligados à catarata e ...

Fig. 4 .. Dinâmica de expressão de genes enriquecidos em lentes em ...

Fig. 4. Dinâmica de expressão de genes enriquecidos em lentes no desenvolvimento embrionário de camundongos.

Genes enriquecidos em lentes em diferentes embrionários ...

Fig. 5 .. Identificação da Ontologia Genética (GO),…

Fig. 5 .. Identificação da Ontologia Genética (GO), via KEGG e termos da Interpro para genes enriquecidos em lentes.

Fig. 6 .. Clusters de expressão de desenvolvimento revelam dinâmica ...

Fig. 6. Os agrupamentos de expressão de desenvolvimento revelam a dinâmica de genes enriquecidos em lentes.

Clustering baseado em algoritmo de árvore auto-organizável (SOTA) ...

Fig. 7 ..Trilhas personalizadas baseadas em dados RNA-seq ...

Fig. 7 .. Trilhas personalizadas baseadas em dados de RNA-seq para o navegador UCSC para visualizar a expressão da lente e ...


Discussão

Neste relatório, apresentamos um RNA-Seq Atlas (Seq-Atlas) para Glycine max usando o sequenciamento Illumina de última geração do transcriptoma de soja. Uma das questões em aberto sobre o método RNA-Seq é o que fazer com sequências de leitura curtas que mapeiam para vários locais em um genoma. Esta questão é particularmente relevante no genoma paleopoliplóide de G. max, que passou por duas rodadas de eventos de duplicação em grande escala no último

60 My que resultou em até quatro regiões de sintaxe dentro da maior parte do genoma [6]. Estudos anteriores indicaram o potencial de sub-representação do número total de contagens de um gene, especialmente em famílias de genes intimamente relacionados [15]. Descobrimos que, desde que estivéssemos cientes das armadilhas potenciais da sub-representação das contagens de genes, informações valiosas sobre a expressão gênica e a relação funcional dos genes poderiam ser obtidas apenas com leituras mapeáveis ​​de maneira única.

Dada a nossa compreensão limitada da complexidade total do genoma da soja, é gratificante que apenas uma pequena porcentagem (3,5%) das leituras mapeadas exclusivamente estavam localizadas fora dos modelos de genes previstos. Isso sugere que a anotação inicial da sequência do genoma da soja capturou a maior parte da atividade transcricional. Usando as informações adicionais sobre as regiões transcricionalmente ativas, o refinamento dos modelos de genes existentes e a capacidade de identificar novos modelos de genes serão melhorados.

Em uma análise da expressão específica do gene em vários tecidos, um dos desafios é superar a grande faixa dinâmica de contagens de expressão gerada pela tecnologia de sequenciamento de próxima geração para identificar genes com perfis de expressão gerais semelhantes. Os dados apresentados aqui têm uma faixa dinâmica de expressão gênica superior a seis ordens de magnitude. Embora um log2-transformação pode reduzir significativamente a faixa dinâmica, um agrupamento hierárquico no log2-dados transformados [11, 16, 17] têm o potencial de perder genes com perfis de expressão gênica altamente semelhantes, mas com expressão gênica significativamente mais baixa ou mais alta em cada tecido. Para identificar todos os genes com perfis de expressão gênica semelhantes, um teste exato de Fisher com uma correção de FDR de 0,05 para um determinado gene foi realizado nas contagens de expressão bruta entre cada tecido e todos os outros tecidos, resultando em uma descrição completa da mudança na expressão gênica. Como o teste Exato de Fisher normaliza para contagens totais no cálculo e a comparação foi entre contagens do mesmo gene e, portanto, têm o mesmo comprimento de gene, as contagens brutas (normalização pré-RPKM) foram usadas. Um agrupamento hierárquico de expressão gênica com base na direção da mudança na expressão e se falha ou não na hipótese nula de que os níveis de expressão são os mesmos entre dois tecidos identifica todos os genes com perfis de expressão semelhantes, independentemente dos níveis de expressão em cada tecido.

Na análise da expressão gênica específica de tecido (Figura 1), determinamos que o padrão geral de expressão gênica se enquadrou em três grupos (Figura 1): tecidos subterrâneos, sementes e aéreos. A semelhança entre este agrupamento usando RNA-Seq e o agrupamento de tecidos transcricionalmente semelhantes em Medicago [11] usando a tecnologia Affymetrix GeneChip valida ainda mais esse resultado. Os tecidos da soja são agrupados por estruturas de plantas intimamente relacionadas: os nódulos são células corticais de raízes modificadas, cada estágio da semente faz parte do desenvolvimento da semente e as vagens, as cascas e as flores são folhas modificadas [41, 42]. Além disso, os estágios de desenvolvimento da semente são mais semelhantes aos tecidos aéreos do que aos tecidos subterrâneos, pois as sementes são mais semelhantes às vagens do que às raízes.

Embora a semelhança do perfil de expressão não implique necessariamente uma função semelhante, pode fornecer uma visão sobre redes co-reguladas de genes. Agrupamentos de genes que são expressos de forma semelhante em tecidos específicos ou estágios de desenvolvimento podem fornecer uma dica quanto ao papel funcional dos genes sem função molecular conhecida. Em um esforço para dividir os dados em partes gerenciáveis, primeiro identificamos genes que foram expressos significativamente em sementes sobre os outros dois grupos de tecido: subterrâneo e aéreo. Em seguida, realizamos uma análise de agrupamento hierárquico para identificar subclados interessantes de genes com perfis de expressão semelhantes no desenvolvimento de sementes. Muitos dos desafios em exibir e interpretar um dendrograma (arquivo adicional 15) foram superados combinando o dendrograma com o log2Os boxplots baseados em cada tecido (Figura 6) resultaram em três clados. O clado 1, o clado 2-1 e o clado 2-2 contêm genes com aumento significativo na transcrição, principalmente nos estágios iniciais, tardios e intermediários de desenvolvimento da semente. Um teste exato de Fisher com uma correção de Bonferroni foi realizado nas categorias de GOslim para genes dos três clados para determinar quais categorias de GOslim estavam sobre-representadas em comparação com as categorias de GOslim para todos os genes no genoma. O clado de desenvolvimento inicial da semente foi super-representado na atividade da beta-glucuronidase, atividade da galactosiltransferase, constituintes estruturais dos ribossomos e atividade da glutamato desidrogenase. O clado de desenvolvimento de semente intermediário (2-2) foi super-representado na atividade da leucocianidina oxigenase, enquanto o clado de desenvolvimento tardio de semente foi super-representado na atividade de reservatório de nutrientes.

Uma vez que a proteína da semente está negativamente correlacionada com o conteúdo e rendimento do óleo da semente [43], os genes com uma função GOslim da atividade do reservatório de nutrientes podem fornecer informações sobre o processo de enchimento da semente. Para entender melhor a extensão do agrupamento de genes com atividade de reservatório de nutrientes no clado de desenvolvimento tardio da semente e para determinar sua relação com o enchimento da semente, identificamos todos os genes (143) em G. max com uma função molecular GOslim correspondente à atividade do reservatório de nutrientes (arquivo adicional 18). Destes genes, 83 são transcricionalmente ativos em nosso conjunto de dados, com uma contagem de transcrição total maior que dois em todos os tecidos. Destes genes transcricionalmente ativos, 19 são encontrados em quatro subclados do clado de desenvolvimento tardio da semente (Figura 6a: números em quadrados). Doze dos genes com atividade de reservatório de nutrientes são encontrados no subclado 2-1: G (Figura 6b). Esses genes são altamente expressos com uma contagem total de transcrição normalizada RPKM em todos os tecidos variando de 39 a 62.401 contagens. Além disso, os genes identificados no clado 2-1 com uma função molecular Goslim de atividade de reservatório de nutrientes são parte do processo de enchimento de sementes, pois a maioria desses genes tem funções baseadas na sequência de consenso provisória de Dana Farber [29, 44] que inclui glicinina , beta-conglicinina e proteína de ligação à sacarose (arquivo adicional 19). Uma vez que os outros genes no clado de desenvolvimento tardio da semente identificados acima têm perfis de expressão semelhantes a esses 19 genes, é provável que existam outros genes no clado de desenvolvimento tardio da semente e, em particular, genes no subclado 2-1: G que possuem ou papéis complementares no enchimento de sementes. É necessária uma análise de dados adicional para elucidar como os outros genes no clado de desenvolvimento tardio da semente se relacionam com os genes do reservatório de nutrientes identificados pelo GOslim e como a compreensão do processo de enchimento da semente melhorará a qualidade, o conteúdo e o rendimento da proteína da semente. Este atlas RNA-Seq fornece um ponto de partida para tal análise.

Como um exemplo final para demonstrar o poder de combinar um atlas de RNA-Seq com a sequência genômica, considere os genes da lipoxigenase de soja (LOXs) [45]. As enzimas lipoxigenases atuam sobre os ácidos graxos poliinsaturados para formar hidroperóxidos de ácidos graxos poliinsaturados que podem ser convertidos em aldeídos e álcoois, o que resulta em uma qualidade de sabor inferior na soja [46, 47]. Genótipos nulos foram identificados em experimentos de irradiação gama que nocautearam os três genes da lipoxigenase: LOX1, LOX2 e LOX3, que se expressaram durante o desenvolvimento da semente [48, 49]. LOX1 e LOX2 estão ligados e encontrados no cromossomo 13, enquanto LOX3 está localizado no cromossomo 15 [45]. o G. max Seq-atlas confirma que para os 72 genes da lipoxigenase (arquivo adicional 20) identificados no genoma da soja e designados com uma função molecular GOslim da atividade da lipoxigenase (GO: 0016165), apenas 3 genes são altamente e significativamente expressos durante o desenvolvimento da semente com base em um Teste exato de Fisher com correção de FDR de 0,05 durante o desenvolvimento da semente. Os genes são: Glyma13g42310, Glyma13g42320 e Glyma15g03030 (Figura 6: números em círculos). Os dados do Seq-Atlas e o último lançamento do genoma apóiam a ligação estreita entre LOX1 e LOX2 no cromossomo 13 - apenas aproximadamente 7.000 pares de bases separam os dois genes. Embora as identidades desses genes da lipoxigenase tenham sido determinadas antes do conhecimento da sequência genômica e do acesso ao sequenciamento de próxima geração [50], não é difícil imaginar como o atlas RNA-Seq poderia ser usado para aumentar a eficiência da descoberta científica.


Métodos

Cultura de células

As células T47D A1–2 foram cultivadas conforme descrito anteriormente 34, autenticadas por perfis de STR e testadas negativas para micoplasma. Para os tratamentos de dexametasona, as células foram inicialmente cultivadas por 24 h em soro reduzido, meio sem hormônio (MEM livre de vermelho de fenol (Gibco 51200) com FBS tratado com 5% de carvão / Dextran (Atlanta S11650), 1X penicilina / estreptomicina (Sigma P0781 ), HEPES a 1% (Sigma H0887), Glutamax 1X (Gibco 35050) e G418 250 ug / ml (Gibco 10131)). Posteriormente, dexametasona (Sigma D4902) foi adicionada ao meio a 100 nM para todos os pontos de tempo de tratamento e 9,5 × 10 −4% de etanol foi usado como controle de veículo. Para ambos os experimentos de RNAseq, as células de controle foram tratadas com etanol por 18 h.

Isolamento de RNA, síntese de cDNA, QPCR e RNAseq

Os kits Qiagen RNeasy com tratamento com DNase na coluna foram usados ​​de acordo com as instruções do fabricante para isolar o RNA total de células A1–2 em triplicado biológico. O Thermofisher SuperScript III foi usado seguindo as instruções do fabricante com Oligo (dT) para sintetizar DNA complementar (cDNA) a partir de 1 ng de RNA total. BioRad ssoAdvanced Universal SYBR Green Supermix foi usado para QPCR. Antes do RNAseq em massa, a qualidade do RNA foi confirmada no Agilent Bioanalyzer 2100 usando o kit RNA 6000 RNA Pico. Bibliotecas de RNAseq totais foram geradas no NIEHS Epigenomics Core usando Ribo-Zero Gold e sequenciadas na Illumina NovaSeq 6000 para mais de 140 milhões de leituras por replicação biológica. Os dados de RNAseq em massa foram processados ​​conforme descrito anteriormente 16.

RNAseq de célula única

Para criar suspensões de células únicas para scRNAseq, as células foram dissociadas usando 0,25% de Tripsina-EDTA (Gibco 25200). O isolador de célula única BioRad ddSeq foi usado para criar emulsões de célula única e o kit Illumina SureCell WTA 3 ′ Library Prep foi usado para processar amostras e gerar bibliotecas de scRNAseq. Para detectar o maior número possível de genes, as bibliotecas foram sequenciadas em alta profundidade na Illumina Novaseq 6000 para obter & gt200 milhões de leituras brutas para cada biblioteca. As leituras foram analisadas quanto ao código de barras da célula e informações de UMI, filtrando as leituras de acordo com as seguintes restrições: não. O segmento de código de barras de 6 nucleotídeos é mais do que uma distância de Hamming de 1 de um bloco de código de barras válido, todos os segmentos de ligante estão dentro de 1 comprimento de nucleotídeo do valor esperado, o UMI tem um comprimento de 8 nucleotídeos e apenas 1 incompatibilidade é permitida dentro dos segmentos ACG-GAC que flanqueiam o UMI. As leituras com as informações do código de barras da célula passando por esses filtros foram então alinhadas ao genoma de referência hg19 usando STAR v2.5.2b 35 e anotadas com o Gencode abrangente v28lift37. Leituras exclusivamente mapeadas foram subsequentemente atribuídas a genes usando featureCounts v1.5.3 36. Os duplicados de PCR foram removidos usando umi_tools 37 no modo por gene. O algoritmo de chamada de joelho de umi_tools foi usado para identificar um ponto de corte de contagem total de transcritos apropriado para células em cada amostra, resultando em um total de 4001 células. As células com porcentagem do gene mitocondrial & gt5% foram removidas e as contagens de células foram equilibradas entre os pontos de tempo por sub-configuração aleatória de cada ponto de tempo para 400 células. O ciclone do pacote scran 38,39 foi usado para determinar as pontuações do ciclo celular para todas as células, e o Seurat v3 40 foi então usado para normalizar e dimensionar os dados de modo que os impactos das contagens de transcritos, porcentagem mitocondrial e pontuações do ciclo celular fossem regredido. Seurat e ggplot2 41 foram usados ​​para visualização de dados.

Estatística e reprodutibilidade

O RNAseq em massa foi realizado usando triplicados biológicos independentes de células A1–2 em cada ponto de tempo do tratamento com Dex. Para scRNAseq, a análise de dados foi realizada usando um conjunto de amostragem reduzida aleatoriamente de 400 células A1–2 para cada ponto de tempo. A análise estatística foi realizada em R usando Seurat para scRNAseq e usando Limma-Voom para RNAseq em massa. Três réplicas biológicas foram usadas para RT-QPCRs e ChIP-QPCRs representados na Fig. 2f, h, e p-valores foram calculados no Excel usando heterocedástico unilateral t-testes. As correlações de Pearson da expressão gênica representadas na Fig. 5 foram calculadas em R.

Resumo do relatório

Mais informações sobre o desenho da pesquisa estão disponíveis no Nature Research Reporting Summary vinculado a este artigo.


Um guia para iniciantes e # x2019s para análise de dados de sequenciamento de RNA

Desde as primeiras publicações que cunharam o termo RNA-seq (Sequenciamento de RNA) apareceu em 2008, o número de publicações contendo dados de RNA-seq cresceu exponencialmente, atingindo um recorde histórico de 2.808 publicações em 2016 (PubMed). Com essa riqueza de dados de RNA-seq sendo gerados, é um desafio extrair o significado máximo desses conjuntos de dados e, sem as habilidades e conhecimentos apropriados, há o risco de má interpretação desses dados. No entanto, uma compreensão geral dos princípios subjacentes a cada etapa da análise de dados de RNA-seq permite que os investigadores sem experiência em programação e bioinformática analisem criticamente seus próprios conjuntos de dados, bem como os dados publicados. Nossos objetivos na presente revisão são quebrar as etapas de uma análise típica de RNA-seq e destacar as armadilhas e pontos de verificação ao longo do caminho que são vitais para cientistas de bancada e pesquisadores biomédicos realizando experimentos que usam RNA-seq.

O sequenciamento de RNA (RNA-seq) foi introduzido pela primeira vez em 2008 (1–4) e na última década tornou-se mais amplamente usado devido aos custos decrescentes e à popularização de núcleos de sequenciamento de recursos compartilhados em muitas instituições de pesquisa. A crescente popularidade do RNA-seq levou a uma necessidade crescente de conhecimento em bioinformática e recursos computacionais. Para que os cientistas de bancada analisem e processem corretamente grandes conjuntos de dados, eles precisam entender os princípios e limitações da bioinformática que vêm com o complexo processo de análise de RNA-seq. Embora a análise de RNA-seq possa ser incrivelmente poderosa e possa revelar muitas novas descobertas interessantes, ela difere das análises usuais a que os cientistas de bancada estão acostumados, pois vem como um conjunto de dados muito grande que não pode ser interpretado sem uma análise extensa.

O protocolo de RNA-seq começa com a conversão de RNA, seja total, enriquecido por mRNA, ou esgotado de rRNA, em cDNA. Após a fragmentação, a ligação do adaptador e a ligação do índice, cada fragmento de cDNA é subsequentemente sequenciado ou "lido" usando uma plataforma de alto rendimento. Os dados lidos brutos são então demultiplexados, alinhados e mapeados para genes para gerar uma tabela de contagens brutas, ponto em que os dados geralmente são entregues ao pesquisador de bancada para iniciar sua própria análise. Ainda não existe um verdadeiro consenso sobre o pipeline mais apropriado para o processamento de dados RNA-seq, no entanto, existem inúmeras ferramentas semiautomáticas online disponíveis, como BaseSpace (Illumina), MetaCore (Thomson Reuters) ou Bluebee (Lexogen). Embora essas ferramentas gerem gráficos de análise de componente principal (PCA), exibam mapas de calor e executem análises de expressão gênica diferencial sem a ajuda de um bioinformático, elas não permitem que os usuários avaliem totalmente a qualidade de seus dados, determinem a precisão de suas próprias análises e adaptar a análise à sua questão biológica, o que pode levar a interpretações errôneas do conjunto de dados. É importante para os investigadores compreender como abordar seu conjunto de dados, apreciar as características de seu conjunto de dados e observar as deficiências nos dados que podem limitar a capacidade de tirar conclusões. Além disso, é imperativo que cada conjunto de dados seja analisado de novo, no sentido de que limites e métodos devem ser adaptados novamente, o que não pode ser alcançado com o uso de aplicativos ou ferramentas online genéricas.

Para os fins deste artigo de métodos, usamos um conjunto de dados de exemplo de um experimento dentro de nosso grupo de pesquisa em que macrófagos alveolares murinos ingênuos foram comparados com aqueles isolados de pulmões transplantados 2 e 24 horas após a fusão. Apresentamos nossa análise usando este conjunto de dados para descrever uma abordagem amigável para análise de RNA-seq para um cientista de bancada.

Masculino Cx3cr1 gfp / + camundongos em um fundo C57BL / 6 e camundongos BALB / c de tipo selvagem com idades entre 12-14 semanas foram usados. Todos os camundongos foram alojados em instalações livres de patógenos específicos. Todos os reagentes foram certificados como livres de endotoxinas pelo fabricante. Todos os estudos foram conduzidos de acordo com as diretrizes do Comitê de Uso e Cuidado de Animais da Universidade Northwestern.

Os transplantes foram realizados entre pares alogênicos incompatíveis doador-receptor, conforme descrito anteriormente (5). Especificamente, pulmões de doadores de Cx3cr1 gfp / + camundongos foram usados ​​como aloenxertos e implantados em recipientes BALB / c de tipo selvagem. Em resumo, os camundongos doadores foram heparinizados e lavados anterógrados através da artéria pulmonar, a traqueia foi ligada após os pulmões serem recrutados e, em seguida, o bloco coração-pulmão foi colhido e mantido a 4 ° C por um período de 2 horas de isquemia fria. As anastomoses para o transplante único de pulmão esquerdo foram realizadas com a técnica com cuff por meio de toracotomia esquerda, o pulmão foi reperfundido e recruta- do e a seguir a toracotomia foi fechada por planos. Os camundongos foram desmamados do ventilador e extubados durante a recuperação, uma vez que estavam ambulatoriais. Em pontos de tempo especificados após a reperfusão, os camundongos receptores foram mortos e o aloenxerto de pulmão foi colhido.

Os pulmões foram processados ​​para suspensões de células únicas conforme descrito anteriormente (5). Resumidamente, o ventrículo direito foi lavado com 10 ml de solução salina balanceada de Hanks gelada, em seguida, os pulmões foram infiltrados com uma mistura de digestão de tecidos contendo colagenase D (Roche) e DNase I (Roche). Foi realizada uma combinação de dissociação mecânica usando o GentleMACS (Miltenyi Biotec) e digestão enzimática a 37 ° C por 30 minutos. As amostras foram então enriquecidas usando microesferas CD45 (Miltenyi Biotec) e sistema AutoMACS (Miltenyi Biotec) antes da coloração de anticorpo.

Ver Tabela E1 no suplemento de dados para anticorpos e diluições usados ​​para coloração de suspensão de célula única e Figura E1 para a estratégia de passagem para classificação de macrófagos alveolares. As células foram classificadas em tampão de classificação de células magneticamente ativado a 4 ° C usando um citômetro de fluxo de quatro laser BD FACSAria II SORP (BD Biosciences).

As células recém-selecionadas foram sedimentadas imediatamente, ressuspensas em 100 μl de tampão de extração PicoPure (Thermo Fisher Scientific) e, em seguida, armazenadas a -80 ° C. O isolamento de RNA foi realizado usando o kit de isolamento de RNA PicoPure (Thermo Fisher Scientific), e amostras com RNA de alta qualidade (número de integridade de RNA, & ​​gt7.0) medido usando o 4200 TapeStation (Agilent Technologies) foram usadas para a preparação da biblioteca. O mRNA foi obtido a partir de RNA total usando kits de isolamento magnético de mRNA NEBNext Poly (A) (New England BioLabs) e as bibliotecas de cDNA foram subsequentemente preparadas usando o kit NEBNext Ultra DNA Library Prep para Illumina (New England BioLabs). As bibliotecas foram sequenciadas em uma plataforma NextSeq 500 usando um kit de sequenciamento de alto rendimento de ponta única de 75 ciclos (Illumina). O sequenciamento rendeu bibliotecas com um tamanho médio de 8 milhões de leituras após o alinhamento. A análise de RNA-seq foi baseada em leituras alinhadas exclusivamente.

As leituras foram demultiplexadas (bcl2fastq), e os arquivos fastq foram alinhados ao genoma do camundongo mm10 (TopHat2 [6]) e mapeados para genes (HTSeq [7]) usando a anotação do gene Ensembl. Comparações de pares entre as várias condições foram executadas usando um modelo log-linear generalizado binomial negativo por meio da função de ajuste glmLRT em edgeR (8, 9).

Os dados de RNA-seq relatados neste artigo foram depositados no Gene Expression Omnibus (GEO) do NCBI e estão acessíveis através do número de acesso da série GEO GSE116583.

Um dos principais objetivos da análise de RNA-seq é identificar genes diferencialmente expressos e co-regulados e inferir o significado biológico para estudos futuros. O material de origem pode ser células cultivadas em vitro, homogenatos de tecido inteiro ou células classificadas. A capacidade de interpretar os resultados depende do projeto experimental apropriado, da implementação de controles e da análise correta. Todo esforço deve ser feito para minimizar o efeito do lote, porque mudanças pequenas e descontroladas em um ambiente podem resultar na identificação de genes diferencialmente expressos (DEGs) não relacionados ao experimento projetado. Fontes de efeito em lote podem ocorrer durante o experimento, durante a preparação da biblioteca de RNA ou durante a execução de sequenciamento e incluem, mas não estão limitados aos listados na Tabela 1. Uma vez que um experimento bem projetado e controlado é realizado, uma abordagem estruturada para o O conjunto de dados permite o controle de qualidade seguido por uma análise imparcial dos dados. Na presente análise, usamos uma abordagem que inclui configuração de filtragem de contagem baixa, estabelecimento de um limite de ruído, verificação de valores discrepantes em potencial, execução de testes estatísticos apropriados para identificar DEGs, agrupamento de genes por padrão de expressão e teste de enriquecimento de ontologia genética (GO) . Para cada um desses componentes de análise, pretendemos destacar importantes pontos de verificação e controles de qualidade que irão agilizar e fortalecer a análise de dados, evitar vieses e permitir que os investigadores usem ao máximo seus conjuntos de dados.

Tabela 1. Fontes de efeito em lote e estratégias propostas para mitigá-los

Para este tutorial, usamos um conjunto de dados compreendendo três grupos de macrófagos alveolares que foram estudados em um modelo murino de transplante de pulmão durante as primeiras 24 horas de reperfusão. Esta abordagem (da qual não fazemos nenhuma reivindicação de originalidade e remetemos o leitor a uma excelente revisão de Conesa e colegas [10] delineando as principais etapas da análise de dados de RNA-seq) permite ao investigador sondar os dados de uma maneira imparcial de uma forma esforço para identificar assinaturas transcricionais e permitir análises posteriores a jusante.

Ao avaliar a variabilidade dentro do conjunto de dados, é preferível que a variabilidade intergrupo, representando diferenças entre as condições experimentais em comparação com as condições de controle, seja maior do que a variabilidade intragrupo, representando variabilidade técnica ou biológica. Uma visão geral dos dados permite a caracterização da variação entre as réplicas e se os grupos experimentais definidos pelo investigador mostram diferenças reais entre os grupos (um grupo é um conjunto de réplicas da mesma condição ou do mesmo tipo de célula). Uma forma de visualizar a variação em um conjunto de dados é por meio do PCA (11). O PCA pega um grande conjunto de dados como entrada e reduz o número de "dimensões" do gene a um conjunto mínimo de dimensões transformadas linearmente, refletindo a variação total do conjunto de dados. Os resultados são comumente apresentados como um gráfico bidimensional em que os dados são visualizados ao longo de eixos que descrevem a variação dentro do conjunto de dados, conhecido como o componentes principais (PCs). PC1 descreve a maior variação dentro dos dados, PC2 a segunda mais e assim por diante. A variação representada por cada PC pode ser calculada como uma porcentagem da variância total e visualizada por um scree plot. Se os dois primeiros PCs não capturarem a maior parte da variação, pode ser útil gerar gráficos PCA bidimensionais adicionais exibindo outros PCs. Desta forma, um gráfico de PCA pode ajudar a visualizar o agrupamento entre as repetições e auxiliar na identificação de outliers técnicos ou biológicos.

Outra abordagem para determinar a variabilidade inter e intragrupo é calcular a distância representada pela correlação entre as amostras. Duas medidas de correlação comumente usadas são o coeficiente de Pearson e o coeficiente de correlação de Spearman (12-14), que descrevem a direcionalidade e a força da relação entre duas variáveis. A correlação de Pearson reflete a relação linear entre duas variáveis ​​responsáveis ​​por diferenças em sua média e DP, enquanto a correlação de classificação de Spearman é uma medida não paramétrica usando os valores de classificação das duas variáveis. Quanto mais semelhantes forem os perfis de expressão para todas as transcrições entre duas amostras, maior será o coeficiente de correlação. Esses coeficientes de correlação são calculados entre todas as amostras e podem ser visualizados como uma tabela ou um mapa de calor, permitindo ao investigador avaliar se as réplicas (técnicas ou biológicas) se agrupam. Além de permitir uma avaliação da variabilidade, tanto o PCA quanto a análise de correlação de amostra podem ajudar a identificar outliers que não foram excluídos durante as etapas anteriores, como o alinhamento. Por exemplo, uma amostra que se alinhou bem e demonstrou uma boa profundidade de leitura pode chegar a esta etapa do pipeline, no entanto, um PCA ou análise de correlação pode identificar esta biblioteca como uma amostra errônea ou contaminada, agrupando o outlier dentro de outro grupo. Também é possível que uma amostra corretamente rotulada caia como um outlier biológico, como se tivesse sido coletada de um animal que se acreditava ter recebido um desafio, mas não apresentou sintomas. Em resumo, essas análises fornecem uma visão geral global de todas as amostras, permitem a determinação de outliers e apresentam dados em um formato fácil de digerir para o investigador e o leitor.

Usando nosso conjunto de dados de macrófagos alveolares, mostramos um gráfico de PCA e um mapa de calor da correlação de Pearson entre as amostras de macrófagos alveolares: ingênuo, transplante 2 horas pós-perfusão e grupos de amostra 24 horas pós-fusão pós-transplante (Figura 1A). Tanto o gráfico PCA quanto o mapa de calor de correlação de Pearson foram gerados usando leituras normalizadas por quilobases de transcrição por 1 milhão de leituras mapeadas (RPKM) contagens (Vejo C onta N ormalizada). O PCA demonstrou agrupamento esperado entre as réplicas dentro das amostras e grupos de amostra espalhados pelos dois PCs. PC1 responde por 68,1% da variação e PC2 responde por 20,3% adicionais. O gráfico de scree (Figura E2) confirmou que a maioria da variância dentro do conjunto de dados foi descrita pelos primeiros dois PCs. Embora o gráfico de PCA enfatize a variabilidade intergrupo, a análise de correlação de Pearson (Figura 1B) fornece uma visão geral de toda a variação entre as amostras, mostrando um valor de correlação de r & gt 0,9 (Tabela 2), consistente com cada grupo pertencente ao mesmo tipo de célula.

Figura 1. Avaliando a variabilidade inter e intragrupo. (UMA) Gráfico de análise de componente principal (PC) exibindo todas as 12 amostras ao longo de PC1 e PC2, que descrevem 68,1% e 20,3% da variabilidade, respectivamente, dentro do conjunto de dados de expressão. A análise de PC foi aplicada a dados de contagem normalizados (leituras por quilobases de transcrição por 1 milhão de leituras mapeadas) e de contagem transformada em log. (B) Gráfico de correlação de Pearson visualizando a correlação (r) valores entre as amostras. A barra de escala representa o intervalo dos coeficientes de correlação (r) exibido.


Etapa 6. Calcule uma pontuação do ciclo celular para cada célula

Isso pode ser usado para determinar se a heterogeneidade na fase do ciclo celular está conduzindo o layout tSNE / UMAP e / ou agrupamento. Isso pode ou não estar obscurecendo o sinal com o qual você se preocupa, dependendo de seus objetivos de análise e da natureza dos dados. (Se necessário, ele pode ser removido em uma etapa posterior.) Também é útil para determinar se certas populações de células são mais proliferativas do que outras. A lista de genes do ciclo celular e o método de pontuação foram retirados de Tirosh I, et al. (2016).


Fundo

A análise genômica de uma única célula tornou possível entender a heterogeneidade celular [1]. Os avanços na pesquisa da genômica de uma única célula também forneceram oportunidades sem precedentes na pesquisa biomédica, onde é importante identificar diferentes tipos de células pertinentes ao envelhecimento e à malignidade celular. Atualmente, eliminar completamente o câncer usando terapias direcionadas molecularmente ainda é uma meta distante para muitos tipos de malignidade. Assim, a investigação de células-tronco cancerosas raras que são resistentes à terapia e o estudo da heterogeneidade intratumoral com respostas diferenciais a drogas em subpopulações celulares distintas fornece uma base para a abordagem desse objetivo [2]. Nos últimos 5 anos, estudos de célula única que objetivaram a escala e a precisão do perfil de todo o genoma de DNA [3], RNA [4], proteína [5], epigenética [6], acessibilidade à cromatina [7] e outros eventos moleculares [8] atingiram dezenas de milhares de células para sequenciamento de RNA de célula única maciçamente paralelo [9] e milhões de células para medições de proteínas de assinatura de citometria de massa [10]. Métodos mais novos e melhores para a realização de análises de células únicas podem capturar a heterogeneidade da população de células, incluindo a natureza heterogênea do câncer, e facilitar a descoberta dos mecanismos moleculares subjacentes.

Embora a análise de dados de sequenciamento de RNA de célula única (scRNA-seq) nos forneça uma oportunidade de estudar a heterogeneidade das células e os genes que são expressos diferencialmente em condições biológicas, é um processo desafiador realizar a análise. Com o rápido aumento nos dados de scRNA-seq, os métodos computacionais precisam superar desafios que vão desde o tratamento de ruído técnico até a construção e caracterização de identidades de células e a análise de linhagem celular por meio da computação de matrizes esparsas de alta dimensão. Portanto, métodos de análise computacional inovadores, eficientes, robustos e escaláveis ​​são essenciais para esta nova fronteira.

Atualmente, o principal obstáculo na análise de dados de scRNA-seq, decorre da baixa eficiência de captura e expressão gênica estocástica, o que aumenta os eventos de dropout de genes em dados de scRNA-seq de todo o genoma. Designamos esses eventos de dropout como eventos de dados perdidos de dados de uma única célula. Estudos anteriores indicam que as taxas gerais de falta são consistentemente altas em alguns dados de célula única. Por exemplo, em uma célula de embrião de camundongo, a taxa de perda pode chegar a quase 30%, mesmo após a redução de ruído [11]. Com uma alta fração de dados ausentes, a exclusão direta dos dados ausentes pode resultar na perda de informações valiosas [12] . Para produzir uma melhor separação de diferentes tipos de células e revelar novas subpopulações biologicamente significativas, várias publicações relataram os dados ausentes como dados censurados e erro falso negativo [13,14,15]. Todas essas metodologias assumem a distribuição dos dados perdidos, entretanto, derivar distribuições de probabilidade adequadas é um problema difícil [12]. Em 2016, Regev et al. observaram que dados ausentes (falsos negativos), falsos positivos e esparsidade de dados podem afetar fortemente as estimativas de heterogeneidade celular, portanto, novos métodos, bem como a adaptação efetiva de algoritmos existentes, são necessários [1]. Além disso, a imputação tradicional de dados perdidos, como filtragem conjunta baseada em usuário e baseada em item, freqüentemente assume que as posições ausentes já são conhecidas na matriz [16]. No entanto, ainda há questões-chave sobre matrizes de expressão de scRNA-seq que precisam ser abordadas. Sem as informações de posição ausentes, os métodos de imputação de dados mencionados acima não podem ser utilizados.

Para resolver os principais problemas na imputação de valor ausente, propusemos um novo modelo com um método de aprendizado de máquina orientado por dados, ou seja, imputação ausente em seq de RNA de célula única (MISC). O MISC foi projetado para resolver três problemas: onde estão os dados ausentes? quantos dados estão faltando? e quais são seus valores ?. Seu início envolve a modelagem do problema para transformar a imputação de dados ausentes em dois problemas de aprendizado de máquina para detecção e imputação dos eventos de dados ausentes. Em seguida, foi proposto um modelo baseado em métodos de classificação e regressão para a solução dos problemas citados. Finalmente, avaliamos o método de imputação ausente em dois conjuntos de dados reais para estudos de diferenciação de células e detecção de tipo de célula.


Resultados

Sequenciamento do transcriptoma e de novo conjunto

Biblioteca de cDNA normalizada preparada a partir do RNA total de R. coreanus Frutos de miquel (20 DAF) foram submetidos à leitura de pares com a plataforma Illumina. Após a remoção das sequências do adaptador, leituras ambíguas e leituras de baixa qualidade (Q20 & # x0003c20), 54,08 milhões de leituras compreendem um total de 4.867.515.000 nucleotídeos (& # x0003e4.86 Gb) foram obtidas para a montagem (Tabela 1). O programa Trinity [32] foi usado para reunir todas as leituras de alta qualidade em um total de 44.619 contigs com um comprimento médio de 755 bp, com metade do comprimento total da montagem em contigs & # x0003e1.1 kb (N50 & # x0200a = & # x0200a1,155 bp). 23.393 contigs (46,57% dos contigs) tinham pelo menos 500 bp de comprimento e 22,97% dos contigs eram maiores que 1.000 bp (Figura S2A). Os contigs foram então reunidos em 43.723 unigenes, com um comprimento médio de 754 bp e um N50 de 1.153 bp, usando a ferramenta TIGR Gene Indices Clustering (Tabela 1). A distribuição dos unigenes seguiu de perto a distribuição dos contigs (Figura S2B). A conservação de sequência e análise de similaridade são bastante úteis para transferir conhecimento de plantas modelo para plantas não modelo [38]. O conjunto de transcrição de KB foi analisado quanto à semelhança com os conjuntos de dados de proteoma para 12 plantas, incluindo Arabidopsis, arroz, pepino e uva usando a pesquisa BLASTX (limite de corte de valor E de & # x022641E-05). Conforme mostrado na Figura 1, o transcriptoma KB mostrou 68,9% de similaridade com morango (Fragaria versa) seguido por 65% com videira (Vitis vinifera) Como previsto, menor número de transcritos KB exibiu similaridade com proteomas de monocotiledôneas (53,7-56%) em comparação com dicotiledôneas (59,4 & # x0201368,9%). A análise da árvore filogenética do conhecido Rosaceae membros da família com base em dados de sequência de DNA de regiões genômicas nucleares e de cloroplasto sugeriram que Fragaria e Rubus são membros da subfamília Rosoideae, que têm o mesmo número cromossômico basal de x & # x0200a = & # x0200a7 [1], [39]. Assim, a homologia significativa entre as proteínas do morango e o transcriptoma KB sugere uma relação evolutiva entre eles, e eles também podem ter funções conservadas.

Os transcritos KB (cobertura & # x0226570%) foram pesquisados ​​contra sequências de proteoma usando o programa BLASTX para analisar similaridade de sequência. A porcentagem de transcrições que mostram similaridade significativa (valor E & # x0003c10 & # x022125) na pesquisa BLASTx são mostradas representadas por um diagrama de barras.

Tabela 1

Framboesa preta coreana
Número total de leituras54,083,500
Nucleotídeos totais (nt)4,867,515,000
Porcentagem GC47.61
Porcentagem Q2096.11
Montagem passo a passo
Trindade
& # x02003Número total da trindade44,619
& # x02003 Comprimento de toda a trindade (nt)33,708,690
& # x02003 Tamanho médio da sequência da trindade (nt)755
& # x02003Trinity N50 (nt)1,155
Unigene
& # x02003Número total de unigenes43,723
& # x02003 Comprimento de todos os unigenes (nt)32,956,464
& # x02003 Tamanho médio da sequência de unigenes (nt)754
& # x02003Unigenes N50 (nt)1,153

Anotação funcional do transcriptoma da framboesa negra coreana

Para anotação de unigenes montados, sequências de genes distintas foram pesquisadas usando BLASTX contra banco de dados NCBI não redundante (NR) com um valor de corte E de 10 & # x022125. De todos os 43.723 unigenes KB, 29.955 (68,51% de todos os unigenes) tiveram resultados BLAST para proteínas conhecidas no banco de dados NR (Tabelas 2 e S2). Além disso, todos os unigenes foram alinhados ao banco de dados público de proteínas, incluindo Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups (COG) e Gene ontology (GO) por BLASTX (cut-off E-value de 10 & # x022125). Conforme resumido na Tabela 2, um total de 30.178 unigenes (69,02%) foram anotados desta maneira, enquanto o restante não pôde ser pareado com genes conhecidos devido à contaminação de sequência e limitação de genoma e informações de EST em Rubus. A contaminação de sequências humanas, bacterianas e virais foi investigada usando a versão baseada na web do DeconSeq (http://deconseq.sourceforge.net/) [40], com uma cobertura de consulta e limiar de identidade de sequência de 90%. Entre as leituras de sequenciamento KB, encontramos baixa contaminação (cerca de 1,11% do total de leituras), indicando que a maioria dos unigenes, que não corresponderam a genes conhecidos, podem ser genes específicos para Rubus gênero e Rosaceae família.

Mesa 2

Banco de dados público de proteínasNo. de hit unigenePercentagem
NR29,95568.51
KEGG15,97736.54
COG10,48023.97
Swiss-Prot21,96450.23
IR11,25325.74
Total30,17869.02

No caso da classificação funcional do gene, um total de 11.253 (25,74%) unigenes foram atribuídos com pelo menos um termo GO. Curiosamente, 18.989 unigenes foram atribuídos com pelo menos um termo GO na categoria biológica, 21.786 unigenes na categoria de componente celular e 11.193 unigenes na categoria de função molecular (Figura S3). Para avaliar ainda mais a integridade da biblioteca de transcriptomas KB, os unigenes montados foram pesquisados ​​contra COG. De 29.955 ocorrências de NR, 10.480 das sequências montadas foram atribuídas à classificação COG (Tabela 2).Entre as 25 categorias funcionais, as cinco categorias principais com cobertura máxima de unigenes são as seguintes & # x0201c Predição de função geral apenas & # x0201d (3.437 unigenes) associadas a funções fisiológicas e metabólicas básicas, & # x0201cTranscrição & # x0201d (2.013 unigenes), & # x0201cReplicação, recombinação e reparo & # x0201d (1.780 unigenes), & # x0201c Modificação Posttranslational, turnover de proteína, chaperones & # x0201d (1.525 unigenes) e & # x0201c Mecanismos de transdução de sinal & # x0201d (1.381 unigenes), enquanto apenas algumas estruturas & #Exacelltrular & # x0201 foram atribuídos x0201d e & # x0201c Estrutura nuclear & # x0201d (Figura S4). Além disso, 534 unigenes foram classificados no grupo de & # x0201cBiossíntese de metabólitos secundários, transporte e catabolismo & # x0201d. Juntos, esses resultados sugerem que o de novo unigenes montados de KB têm ampla cobertura de transcriptoma e fornecem um recurso valioso para facilitar a descoberta de novos genes envolvidos em processos fisiológicos e de desenvolvimento específicos.

Perfil de metabólito de frutas coreanas de framboesa preta em diferentes estágios de amadurecimento

Perfil de metabólito baseado em espectrometria de massa de KB em três estágios de amadurecimento, 15 dias pós-antese (DPA) (estágio 1, comumente usado na medicina tradicional à base de ervas), 20 DPA (estágio 2) e 25 DPA (estágio 3, usado como fruta fresca e para processamento), foi realizada para investigar as mudanças na composição metabólica (Tabela S3 e S4). Para apoiar o reconhecimento de padrões das diferenças metabólicas em amostras de KB de acordo com os diferentes estágios de amadurecimento, os conjuntos de dados foram analisados ​​estatisticamente por análise de componente principal (PCA) e análise discriminante de mínimos quadrados parciais (PLS-DA). Conforme mostrado na Figura 2, os gráficos de pontuação PLS-DA de ambos os conjuntos de dados apontam para o perfil de metabólito diferencial entre as amostras colhidas em diferentes estágios. Padrões semelhantes foram mostrados no PCA (Figura S5). Na análise de GC-IT-MS, os frutos do estágio 1 (dimensão PLS1 positiva) foram claramente separados do estágio 3 (dimensão PLS1 negativa) ao longo do PLS1 (21,0%), enquanto os frutos do estágio 2 foram separados dos estágios 1 e 3 ao longo do PLS2 (12,9% ), (Figura 2A). Um total de quinze metabólitos foram selecionados como variáveis ​​diferenciais usando o valor VIP (importância variável na projeção) (VIP & # x0003e0.7) e p valor (p& # x0003c0.05) do conjunto de dados PLS-DA (Tabela 3). Metabólitos primários, incluindo ácidos orgânicos, açúcares, ácidos graxos e aminoácidos, foram encontrados para influenciar a formação de agrupamentos de perfis de GC-IT-MS. O mapa de calor representa a distribuição diferencial de metabólitos durante o amadurecimento da fruta (Figura 3). Em comparação com a expansão dos frutos, o progresso do amadurecimento dos frutos tem sido caracterizado por uma diminuição na sacarose, ácidos orgânicos e GABA [41]. Da mesma forma, nosso achado demonstrou a diminuição da sacarose, a maioria dos aminoácidos, ácidos orgânicos e ácidos graxos, enquanto a frutose e a glicose aumentaram em concentração relativa durante o amadurecimento dos frutos, indicando que o KB sofreu grandes mudanças no metabolismo do carbono-nitrogênio durante o amadurecimento dos frutos.

Os conjuntos de dados obtidos por GC-IT-MS (A) e UPLC-Q-TOF-MS (B) foram analisados ​​por PLS-DA.

23 metabólitos que mostram mudanças estatisticamente significativas durante o processo de amadurecimento (VIP & # x0003e0.7 e p& # x0003c0.05) foram representados por Mapa de calor, com intensidades relativas indicadas pela escala de calor.

Tabela 3

Não.RT (min) 1 Valor VIPÍon identificado (m / z) 2 Metabólitos putativos 3 Derivatizado 4 p-valorID 5
14.411.21147ácido oxálico(TMS)20.005STD
26.081.15299ácido fosfórico(TMS)30.003STD
38.911.09147ácido málico(TMS)30.021STD
49.320.95147, 232Ácido L-aspártico(TMS)30.023STD
59.421.33156ácido piroglutâmico(TMS)20.005STD
69.461.17174& # x003b3-ácido aminobutírico (GABA)(TMS)30.002STD
712.661.46273Ácido Cítrico(TMS)40.005STD
813.161.33217D-frutoseMeOX (TMS)50.034STD
913.421.70148, 217D-glicoseMeOX (TMS)50.000STD
1014.061.24281ácido gálico(TMS)40.000STD
1115.021.07313Ácido palmíticoTMS0.013STD
1216.571.49339Ácido oleicoTMS0.001STD
1316.781.63117,341ácido esteáricoTMS0.000STD
1419.721.16361sacarose(TMS)80.004STD
1523.651.44237& # x003b1-tocoferolTMS0.000STD

O gráfico de pontuação resultante do PLS-DA da análise UPLC-Q-TOF-MS mostrou que os frutos de diferentes estágios de maturação foram claramente distinguidos pela combinação de PLS1 (27,4%) e PLS2 (10,3%), (Figura 2B). Oito metabólitos secundários derivados de flavonóides, incluindo flavonóis, flavonóis, proantocianidinas e antocianinas, foram metabólitos significativos para a análise relativa entre os três estágios de maturação das frutas (Tabela 4). A quantidade de flavonóis (catequina e epicatequina) e dímeros de proantocianidina tipo B diminuiu, enquanto os flavonóis (quercetina 3-O-rutinósido e quercetina-glucuronido) e antocianinas (cianidina 3-O-xilosilrutinósido e cianidina 3-O-rutinósido) aumentaram durante o progresso do amadurecimento (Figura 3). Em particular, verificou-se que a principal antocianina no amadurecimento KB são os derivados de cianidina, que foram altamente acumulados em frutas vermelhas escuras (estágio 3).

Tabela 4

UPLC-Q-TOF-MSLC-IT-MS / MS
Não.RT 1 Valor VIPMassa experimental [M-H] & # x02212Fórmula& # x00394ppm[M-H] & # x02212 [M + H] +Íons de fragmento de MS n (m / z)UV (nm)Tentativamente metabólitos 2 p-valorID 3
13.061.71725.1926C32H37O19& # x022120.4725727727 e # x0003e287 e # x0003e213276, 519cianidina 3-O-xilosilrutinósido0.000DB
23.131.72593.1504C27H29O15& # x022120.3593595595 & # x0003e287 & # x0003e213278, 518cianidina 3-O-rutinósido0.000STD
33.161.67577.1336C30H25O12& # x022121.7577579577 e # x0003e289 e # x0003e188248, 314dímero de proantocianidina B10.000DB
43.281.60289.0712C15H13O60.0289291291 e # x0003e123274, 306catequina0.000STD
53.381.62577.1348C30H25O120.3577579577 e # x0003e289 e # x0003e188242, 301dímero de proantocianidina B20.000DB
63.571.32289.0711C15H14O6& # x022120.3289291291 e # x0003e123272, 300epicatequina0.000STD
73.901.65609.1461C27H29O160.8609611609 e # x0003e301 e # x0003e151275, 351quercetina 3-O-rutinosídeo0.000STD
84.180.95477.0679C21H17O132.1477479477 e # x0003e301 e # x0003e179269,358quercetina-glucuronídeo0.002DB

Identificação de enzimas candidatas envolvidas na biossíntese de antocianinas

O acúmulo de derivados de cianidina durante o processo de amadurecimento nos levou a investigar se esse fenômeno é mediado pela alteração da expressão de genes da via metabólica dos flavonóides. A via de biossíntese de antocianina detalhada foi ilustrada na Figura 4A [42]. A via biossintética da antocianina é geralmente dividida em duas seções, a primeira e a posterior. De acordo com nossa previsão, com base na atribuição da via KEGG, encontramos 28 unigenes que codificam 23 enzimas putativas envolvidas na biossíntese de antocianinas da biblioteca de transcriptomas KB (Tabela 5). Embora a maioria deles tivesse sequências parciais de nucleotídeos, pudemos analisar seu padrão de expressão durante o processo de amadurecimento usando qRT-PCR. Conforme ilustrado na Figura 4B, o nível de expressão da maioria dos unigenes, como chalcona sintases (CHSs), chalcona flavanona isomerases (CHIs) e flavanona 3-hidroxilases (F3Hs) regulando a biossíntese de diidrokaempferol aumentou significativamente quando as frutas ficaram vermelho escuro (estágio 3) . Da mesma forma, o nível de transcrição aprimorado de DFRs (dihidroflavonol 4-redutases), LDOXs (dioxigenases de leucoantocianidina) e 3GTs (antocianidina 3-O-glucosiltransferases) unigenes também foi observada, o que leva à formação de antocianinas. Além disso, a expressão de F3 e # x02032H1 (flavonóide 3 & # x02032-hidroxilase1) mostrou um aumento dramático durante o processo de amadurecimento. Por outro lado, F3 & # x020325 & # x02032H (flavonóide 3 & # x020325 & # x02032-hidroxilase) o unigene foi regulado negativamente de 15 DPA para 25 DPA. Isso é consistente com o acúmulo de derivados de cianidina, mas não de derivados de delfinidina (Figura 4C). Uma especulação é que o aumento do nível de F3 e # x02032H1 pode ter acumulado derivados de cianidina como principais compostos de antocianina devido à competição enzimática pelo substrato diidrokaempferol.

(A) Via biossintética das antocianinas. Os nomes das enzimas foram abreviados da seguinte forma: chalcona sintase (CHS), chalcona isomerase (CHI), flavanona 3-hidroxilase (F3H), flavonóide 3 & # x02032-hidroxilase (F3 & # x02032H), flavonóide 3 & # x020325 & # x02032-hidroxilase (F3H) x020325 & # x02032H), dihidroflavonol 4-redutase (DFR), leucoantocianidina dioxigenase (LDOX) e antocianidina 3-O-glucosiltransferase (3GT). (B) Padrão de expressão de genes envolvidos na via de biossíntese de antocianinas. Os níveis de expressão dos genes do estágio 2 e 3 de amadurecimento foram comparados ao estágio 1. As médias e erros padrão foram calculados a partir de três medições independentes. Teste t de Student em comparação com o Estágio 1, * P & # x0003c0.05, ** P & # x0003c0.01 e *** P & # x0003c0.001. (C) Os conteúdos de antocianina determinados por UPLC-Q-TOF-MS foram representados por gráficos de box-whisker. Os metabolitos foram relativamente quantificados por suas áreas de pico m / z com o instrumento e processamento de dados. As barras de erro são os desvios padrão de três medições independentes.

Tabela 5

Rubus coreanus
GeneNúmero CENome do gene 1 NUMNPUNGSN Fragaria vesca 2
Chalcone sintaseEC: 2.3.1.74CHS22-FV7G02590
FV7G02600
Chalcona isomeraseEC: 5.5.1.6CHI44-FV7G31110
FV7G25290
FV2G25040
FV3G19830
Flavanona 3-hidroxilaseEC: 1.14.11.9F3H22 <"type": "entrez-nucleotide", "attrs": <"text": "EU078685", "term_id": "158578344" >> EU078685FV1G13680
<"type": "entrez-nucleotide", "attrs": <"text": "EU255776", "term_id": "166208404" >> EU255776FV6G51040
Diidroflavonol-4-redutaseEC: 1.1.1.219DFR54-FV2G34030
FV7G11200
FV3G23500
FV2G50610
Flavonóide 3 & # x02032-hidroxilaseEC: 1.14.13.21F3 e # x02032H53-FV5G12250
FV1G10430
FV0G15960
Flavonóide 3 & # x02032,5 & # x02032-hidroxilaseEC: 1.14.13.88F3 & # x020325 & # x02032H11-FV5G00680
Leucoantocianidina dioxigenaseEC: 1.14.11.19LDOX33-FV5G01390
FV5G18660
FV7G29560
Antocianidina 3-O-glucosiltransferaseEC: 2.4.1.1153GT64-FV7G33970
FV7G06650
FV5G37930
FV3G11420
Total 28232

Para gerar um mapa de fluxo da biossíntese de antocianinas durante o processo de amadurecimento, a análise do fluxo metabólico foi realizada usando a ferramenta YANA [43]. A mudança na biossíntese de antocianina resultou das expressões gênicas diferenciais durante o processo de amadurecimento foi mostrado na Figura S6, o que sugere que DFR4 e LDOX1 são enzimas principais para a biossíntese de derivados de cianidina.

Caracterização funcional da chalcona isomerase em R. coreanus Transcriptoma de Miquel

A distribuição de CHI em plantas superiores resulta em mudanças fenotípicas de cor [44] e # x02013 [46]. Além disso, foi sugerido que uma das principais limitações na via biossintética dos flavonóides no tomate é a falta de CHI expressão [47], [48], indicando a importância desta enzima para o fluxo da via dos flavonóides, incluindo a biossíntese de antocianinas. Da biblioteca do transcriptoma KB, encontramos um unigene CHI completo (RcMCHI2, Unigene18325), juntamente com três sequências de fragmentos que codificam diferentes proteínas CHI. Os CHIs são classificados principalmente em dois tipos, CHI tipo I e CHI tipo II [49]. Os CHIs do tipo I isomerizam apenas naringenina chalcona para produzir naringenina, enquanto os CHIs do tipo II podem catalisar naringenina chalcona e isoliquiritigenina em naringenina e liquiritigenina, respectivamente [49]. Uma árvore filogenética gerada pelo método de união de vizinhos com base nas sequências de aminoácidos exibiu que RcMCHI2 está evolutivamente relacionado ao grupo CHI tipo I (Figura S7). o Arabidopsis transparente testa 5-1 (tt5-1) mutante, sem atividade CHI, foi selecionado para caracterizar funcionalmente o RcMCHI2. RcMCHI2 transportando Myc-tag foi colocado sob o controle do promotor 35S e transformado no Arabidopsis tt5-1 mutante. A expressão do transgene foi confirmada por western blotting com anticorpo anti-marcador myc, e quatro linhas independentes foram selecionadas para investigações de fenótipo (Figura 5A). As sementes coletadas de linhagens transgênicas selecionadas apresentaram pigmentação característica do tipo selvagem. Arabidopsis, enquanto as sementes do Arabidopsis tt5-1 mutantes eram de cor amarela (Figura 5B). Sementes secas de Arabidopsis tt5-1 mutante em excitação com UV-A exibiu uma fluorescência azul esverdeada que foi resgatada com a expressão de RcMCHI2. Além disso, essas plantas transgênicas mostraram pigmentos de antocianina nos cotilédones e hipocótilos quando cultivadas sob condição de deficiência de nitrogênio, ao contrário de Arabidopsis tt5-1 mudas (Figura 5B). Além disso, a expressão ectópica de RcMCHI2 no Arabidopsis tt5-1 linhas acumulam quantidade semelhante de delfinidina 3-O-rutinósido e cianidina 3-O-rutinósido, em comparação com o tipo selvagem Arabidopsis mudas (Figura 5C Tabela S5). Essas descobertas sugerem fortemente que RcMCHI2 codifica a enzima CHI funcional, que apóia fortemente que nossa biblioteca de transcriptomas KB fornece uma base sólida para futuras pesquisas funcionais e genômicas.

(A) A expressão de MYC-RcMCHI2 foi confirmada por western blotting com anticorpo anti-marcador myc em plantas T1 com duas semanas de idade. (B) Restauração da pigmentação da antocianina e da cor do tegumento da semente no transgênico Arabidopsis tt5-1 mutantes expressando RcMCHI2. Mudas foram cultivadas em meio com baixo teor de nitrogênio para induzir o acúmulo de antocianina. As sementes secas maduras foram expostas a UV. (C) A análise de antocianinas em extratos de transgênicos RcMCHI2 Arabidopsis tt5-1 linhas. A área relativa de antocianinas em mudas T2 com 1 semana de idade foi calculada pelo software MassLynx. As barras de erro são os desvios padrão de três medições independentes. Teste t de Student comparado com tt5 mutante, *** P & # x0003c0.001.


DISCUSSÃO

Desenvolvemos aqui uma abordagem simples e intuitiva, AssociVar, para (a) detectar genuíno mutações do sequenciamento da população MinION, e (b) inferir o conjunto de haplótipos (cepas) presentes em uma população. Nossa abordagem é baseada na noção de que os erros de sequenciamento serão dispersos aleatoriamente ao longo das leituras, enquanto as mutações reais tendem a se associar a origens genéticas específicas. No caso em que as taxas de erros técnicos são altas (como ocorre com o MinION), isso permite focar na real diversidade genética que está oculta na vasta gama de erros técnicos gerados por este método. Notavelmente, nossa abordagem é geral o suficiente para que possa ser usada para qualquer tipo de sequenciamento de leitura longa.

Aplicamos AssociVar aos dados de sequenciamento de uma população evoluída de fagos onde o sequenciamento Illumina estava disponível, o que nos permitiu corroborar se as mutações que encontramos com base na análise dos dados do MinION por si só eram de fato reais. Surpreendentemente, todas, exceto uma das mutações de alta frequência observadas nos dados de p15A e p15B (& gt10%) foram obtidas usando AssociVar, apesar do fato de que o percentil 99 para erros técnicos foi tão alto quanto 43% (Tabela 1). Na verdade, apesar da taxa de deleção muito alta, AssociVar identificou com precisão a única mutação de deleção real presente em nossas populações, sugerindo uma sensibilidade muito alta do método. Nossa abordagem também mostra alta especificidade, com uma taxa de falsos positivos inferior a 0,1%. Finalmente, mostramos que o uso de uma abordagem ingênua baseada em um limite de frequência como um ponto de corte para separar mutações reais de erros resulta em taxas de falsos positivos extremamente altas, demonstrando o valor de nossa abordagem.

Originalmente, ao observar os dados na Figura 1, como uma primeira aproximação, parecia provável que as mutações com uma frequência semelhante seriam mutações compartilhadas nos mesmos genomas. Por conseguinte, tínhamos a hipótese de que pelo menos dois grupos de mutações na linha B (T1764- / G3114A / A1664G e A1611G / A1744G / T1440C / G1906A / A535G) estariam presentes nos mesmos genomas. Isso acabou sendo apenas parcialmente verdadeiro: as mutações com frequência semelhante às vezes estavam de fato nos mesmos genomas (por exemplo, T1764- / G3114A), mas às vezes não completamente (os dois anteriores e A1664G) (Tabela 2). Esses resultados ilustram a utilidade do MinION para resolver as relações entre as mutações e sua vantagem para diferenciar variantes com mutações exibindo frequências semelhantes.

Usamos ainda nossa abordagem para realizar a análise reversa: ao analisar o mRNA do gene enolase de levedura, nossa análise sugeriu que a população de mRNA sequenciada não era homogênea. Isso foi então precisamente verificado pelo sequenciamento Illumina da mesma população. Notavelmente, esta análise mostra que (a) AssociVar pode ser usado para analisar diferentes tipos de dados, variando de genomas de vírus a mRNA de qualquer organismo, e (b) AssociVar pode ser usado sem sequenciar uma sequência de controle. Notamos que isso requer mais cuidado, uma vez que nossa análise de MS2 mostrou que associações espúrias entre mutações podem ser criadas artificialmente pelo próprio processo de sequenciamento. O uso de AssociVar sem uma sequência de controle requer que o usuário especifique o limite da estatística qui-quadrado normalizada. Como acontece com todos os métodos, a especificidade do AssociVar tem o custo da sensibilidade e vice-versa (Figura 6). No entanto, parece que a melhor estratégia que podemos sugerir é usar um limite muito alto, o que é extremamente eficaz para variantes com frequência superior a 5 ou 10%.

É importante delinear as limitações de nossa abordagem. Notamos que não podemos distinguir haplótipos / cepas que diferem entre si em apenas uma posição, porque nosso método se baseia na associação entre duas posições contendo mutações reais. Da mesma forma, se duas cepas diferem em loci muito proximais, AssociVar também falhará, uma vez que filtramos associações entre mutações que estão separadas por & lt15 pares de bases. Postulamos que as associações presumivelmente artefatuais que observamos entre os loci proximais são induzidas pelo RNA (ou DNA) passando pelo poro do sequenciador. Finalmente, também observamos padrões específicos de mutações que foram reproduzidas entre nossa sequência de controle e as duas populações evoluídas de MS2. Isso sugere duas possibilidades: primeiro, talvez o contexto da sequência e / ou a estrutura secundária do RNA induzam erros específicos no MinION e, segundo, é possível que os genomas MS2 sofram modificações de RNA e essas sejam a causa desses erros específicos. O sequenciamento direto de RNA do MinION registra o sinal elétrico bruto produzido pelo RNA que atravessa os poros, e isso potencialmente oferece a oportunidade de identificar modificações de RNA usando uma ferramenta recentemente desenvolvida chamada Tombo (versão 1.5) por Oxford Nanopore (https: //nanoporetech.github .io / tombo /). Infelizmente, não foi possível determinar de forma conclusiva a presença ou o efeito das modificações do RNA e sua relação com as mutações associadas. Nossos resultados sugerem que Tombo ainda sofre de uma alta taxa de falsos positivos, enquanto a taxa de verdadeiros positivos do método ainda não foi determinada (38).O primeiro foi demonstrado aqui por um grande número de locais presumivelmente modificados encontrados no gene da levedura enolase, apesar do fato de que esse gene foi criado sinteticamente in vitro, onde as modificações provavelmente não ocorreriam. No entanto, analisamos nossas amostras de MS2 e encontramos um padrão semelhante de modificações presumíveis entre as três amostras de MS2, mas não houve correlação entre os locais com uma alta taxa de modificação e os locais com pontuações de chi normalizadas altas pela AssociVar (ver texto suplementar, Figuras suplementares S7 –S10). Embora não possamos descartar que as modificações de RNA são responsáveis ​​pelo padrão de erros no MinION, concluímos que mais pesquisas são necessárias para determinar quais fatores induzem esses erros.

Embora nosso método seja ideal para RNA direto ou sequenciamento direto de DNA, também usamos o método para cDNA que foi amplificado do RNA no caso da análise do vírus Zika (30) (Figura Suplementar S4). Quando tentamos reconstruir os haplótipos conhecidos presentes nesta amostra, nosso método não teve sucesso total para recapitular os haplótipos (dados não mostrados). Uma possível explicação para isso é que durante a etapa de amplificação, ou sequências quiméricas de ambas as cepas foram criadas, ou ocorreu recombinação por PCR, quebrando algumas das ligações entre os locais. Em tais casos, o uso de AssociVar é limitado à detecção de mutações apenas, e isso sugere ainda que o sequenciamento direto de RNA / DNA pode ser preferível.

Para resumir, antecipamos que, devido à sua facilidade de uso e vantagens listadas acima, o sequenciamento de leitura longa direta usando MinION será cada vez mais valioso no campo da genética de vírus e em diversos campos adicionais, como estudos de transcriptoma, genética do câncer e microbiologia. A abordagem da AssociVar que sugerimos aqui é simples e aplicável a qualquer organismo e, como tal, esperamos que seja uma adição útil à caixa de ferramentas de genômica em vários campos.


Assista o vídeo: W20: Single-Cell RNA-Seq Analysis with Python - Day 3 (Novembro 2021).