Questão:
Derive um GTF contendo genes codificadores de proteínas de um arquivo GTF com Exons e CDS
Tom Kelly
2018-12-10 17:48:42 UTC
view on stackexchange narkive permalink

Por que preciso de um arquivo compatível

Estou tentando executar o velocyto com o pacote R para analisar a velocidade do RNA (trajetórias celulares) com dados RNASeq de uma única célula. Realizei uma análise de célula única a partir de dados 10x Genomics usando cellranger.

Alinhei com sucesso as leituras para obter arquivos de tear e importei-os para R. Posso obter a velocidade desses arquivos seguindo vinhetas. No entanto, não posso reproduzir a análise de velocidade do RNA com base na "estrutura do gene".

Estou trabalhando com um organismo diferente do exemplo (não humano ou camundongo), portanto os dados de anotação fornecidos no exemplo não trabalhos. Eu tenho um arquivo GTF para a última anotação deste organismo. No entanto, ele contém apenas “exon” e “CDS” como recursos. Essa parece ser a origem do problema. A função “find.ip.sites” requer um GTF com “features” = “gene” e um dos “atributos” como “protein_coding”. Esses requisitos são embutidos em código na função velocyto.R.

Eu tenho os seguintes arquivos GTF do conjunto de dados AtRTD2. Os rótulos dos cromossomos correspondem aos arquivos de tear em R. 

  Chr1 TAIR10 exon 3631 3913. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 exon 3996 4276. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 exon 4486 4605. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 exon 4706 5095. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 exon 5174 5326. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 exon 5439 5899. +. transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 3760 3913. + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1 TAIR10 CDS 3996 4276. + 2 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 4486 4605. + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 4706 5095. + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 5174 5326. + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 TAIR10 CDS 5439 5630. + 0 transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010"; Chr1 Araport11 exon 6788 7069. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Observe "ARV1"; Chr1 Araport11 exon 7157 7450. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 7564 7649. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 7762 7835. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 7942 7987. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 8236 8325. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 8417 8464. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 exon 8571 8737. -. transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 CDS 7315 7450. - 1 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Observe "ARV1";
Chr1 Araport11 CDS 7564 7649. - 0 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 CDS 7762 7835. - 2 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 CDS 7942 7987. - 0 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1"; Chr1 Araport11 CDS 8236 8325. - 0 transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Nota "ARV1";

As especificações para um GTF ou GFF2 são:

Campos

Os campos devem ser separados por tabulação. Além disso, todos, exceto o campo final em cada linha de recurso, devem conter um valor; colunas "vazias" devem ser denotadas com um '.'

  1. seqname - nome do cromossomo ou estrutura; os nomes dos cromossomos podem ser fornecidos com ou sem o prefixo 'chr'. Nota importante: o seqname deve ser usado dentro do Ensembl, ou seja, um nome de cromossomo padrão ou um identificador do Ensembl, como uma ID de cadafalso, sem qualquer conteúdo adicional, como espécie ou montagem. Veja o exemplo de saída do GFF abaixo.

  2. fonte - nome do programa que gerou este recurso ou a fonte de dados (banco de dados ou nome do projeto)

  3. recurso - nome do tipo de recurso, por exemplo Gene, Variation, Similarity

  4. início - Posição inicial do elemento, com numeração sequencial começando em 1.

  5. final - Posição final do elemento, com numeração sequencial começando em 1.

  6. pontuação - Um valor de ponto flutuante.

  7. vertente - definido como + (avançar) ou - (retroceder).

  8. frame - Um de '0', '1' ou '2'. '0' indica que a primeira base do recurso é a primeira base de um códon, '1' que a segunda base é a primeira base de um códon e assim por diante.

  9. attribute - Uma lista separada por ponto-e-vírgula de pares de valores de tags, fornecendo informações adicionais sobre cada recurso.

O que estou procurando para

Existe uma maneira de derivar um arquivo GTF compatível contendo genes codificadores de proteínas dos exons e CDS? Quero produzir um GTF que contenha genes nas características e codificação_proteica nos atributos. É possível fazer isso com ferramentas ou scripts existentes?

O que tentei até agora

Posso modificar o código-fonte da função “find.ip.sites” para executar meu arquivo GTF com esses recursos ausentes. No entanto, isso requer a execução de funções internas do pacote escrito em Rcpp e significa que meu fluxo de trabalho substitui futuras atualizações do pacote. Executar o resto da vinheta retorna erros, pois nenhum íntron longo o suficiente foi identificado (definir os limites mais baixos também é incompatível com os modelos lineares gerais chamados). Portanto, acho melhor gerar um arquivo GTF ou GFF3 compatível do que alertar o código-fonte das funções.

Embora destinadas a arquivos GTF, as funções do pacote podem importar arquivos GFF3, apesar de serem descritas para GTF na documentação. Eu tentei gerar um arquivo GFF3 usando gffread de Cufflinks e substituindo “feature” = “mRNA” por “gene” e adicionando “protein_coding”. Isso também retorna erros ao executar o algoritmo de velocidade. Ele não funciona para o arquivo GTF usado como entrada para cellranger ou para o arquivo gerado por ele. Existe uma maneira de anotar genes que codificam proteínas e limites de íntron / exon com base em um arquivo GTF contendo apenas anotações de exon e CDS?

Com botões de punho 2.2.1, um GFF3 foi gerado 

  gffread - E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gffsed -i '/ ARNm / s / gene_name = AT / gene_type = protein_coding; gene_name = AT / g' AtRTD2_19April2016.gff'sed -i 's / ARNm / gene / g 'AtRTD2_19April2016.gff

Estes são os arquivos GFF3 que experimentei: 

  # gffread - E AtRTD2_19April2016.gtf -o- ## gff-version3Chr1 TAIR10 gene 3631 5899. +. ID = AT1G01010.1; geneID = AT1G01010; gene_type = proteína_coding; gene_name = AT1G01010Chr1 TAIR10 exon 3631 3913. +. Pai = AT1G01010.1Chr1 TAIR10 exon 3996 4276. +. Pai = AT1G01010.1Chr1 TAIR10 exon 4486 4605. +. Pai = AT1G01010.1Chr1 TAIR10 exão 4706 5095. +. Pai = AT1G01010.1Chr1 TAIR10 exon 5174 5326. +. Pai = AT1G01010.1Chr1 TAIR10 exon 5439 5899. +. Pai = AT1G01010.1Chr1 TAIR10 CDS 3760 3913. + 0 Pai = AT1G01010.1Chr1 TAIR10 CDS 3996 4276. + 2 Pai = AT1G01010.1Chr1 TAIR10 CDS 4486 4605. + 0 Pai = AT1G01010.1Chr1 TAIR10 CDS 4706 5095. + 0 Pai = AT1G01010.1Chr1 TAIR10 CDS 5174 5326. + 0 Pai = AT1G01010.1Chr1 TAIR10 CDS 5439 5630. + 0 Pai = gene AT1G01010.1Chr1 Araport11 6788 8737. +. ID = AT1G01020_P2; geneID = AT1G01020; gene_type = proteína_coding; gene_name = AT1G01020 Chr1 Araport11 exon 6788 7069. -. Pai = AT1G01020_P2Chr1 Araport11 exon 7157 7450. -. Pai = AT1G01020_P2Chr1 Araport11 exon 7564 7649. -. Pai = AT1G01020_P2Chr1 Araport11 exon 7762 7835. -. Pai = AT1G01020_P2Chr1 Araport11 exon 7942 7987. -. Pai = AT1G01020_P2Chr1 Araport11 exon 8236 8325. -. Pai = AT1G01020_P2Chr1 Araport11 exão 8417 8464. -. Pai = AT1G01020_P2Chr1 Araport11 exon 8571 8737. -. Pai = AT1G01020_P2Chr1 Araport11 CDS 7315 7450. - 1 Pai = AT1G01020_P2Chr1 Araport11 CDS 7564 7649. - 0 Pai = AT1G01020_P2Chr1 Araport11 CDS 7762 7835. - 2 Pai = AT1G01020_P2
Chr1 Araport11 CDS 7942 7987. - 0 Pai = AT1G01020_P2Chr1 Araport11 CDS 8236 8325. - 0 Parent = AT1G01020_P2

As especificações para um GTF ou GFF2 são:

As especificações para um GF3 são:

Campos

Os campos devem ser separados por tabulação. Além disso, todos, exceto o campo final em cada linha de recurso, devem conter um valor; colunas "vazias" devem ser denotadas com um '.'

  1. seqid - nome do cromossomo ou estrutura; os nomes dos cromossomos podem ser fornecidos com ou sem o prefixo 'chr'. Nota importante: o ID de seq deve ser usado dentro do Ensembl, ou seja, um nome de cromossomo padrão ou um identificador do Ensembl, como um ID de cadafalso, sem qualquer conteúdo adicional, como espécie ou montagem. Veja o exemplo de saída GFF abaixo.

  2. fonte - nome do programa que gerou este recurso, ou a fonte de dados (banco de dados ou nome do projeto)

  3. tipo - tipo de característica. Deve ser um termo ou acesso da ontologia de sequência SOFA

  4. início - posição inicial do recurso, com numeração de sequência começando em 1.

  5. final - Posição final do elemento, com numeração sequencial começando em 1.

  6. pontuação - Um valor de ponto flutuante.

  7. vertente - definido como + (para frente) ou - (reverso).

  8. fase - Um de '0', '1' ou '2'. '0' indica que a primeira base do recurso é a primeira base de um códon, '1' que a segunda base é a primeira base de um códon e assim por diante.

  9. atributos - uma lista separada por ponto-e-vírgula de pares de valor de tag, fornecendo informações adicionais sobre cada recurso. Algumas dessas tags são predefinidas, por exemplo, ID, nome, alias, pai - consulte a documentação GFF para obter mais detalhes.

Provavelmente ajudaria se você pudesse [editar] sua pergunta e nos mostrar algumas linhas de cada um dos dois arquivos GTF para que possamos ver o que você quer dizer com mais clareza.
Você provavelmente pode usar AWK para isso, é o recurso de "gene" semelhante a "CDS" ou é do "início do primeiro exon" ao "final do último exon".
Acho que a pergunta e o "O que estou procurando" são um pouco confusos. Especificamente, parece que você deseja gerar um GTF em "O que estou procurando", mas deseja obter a sequência de aminoácidos de proteínas de um GTF no título. Existe alguma maneira de esclarecer ou tornar o título e "O que estou procurando" consistente?
Temos um [GTF completo para Arabidopsis em Plantas Ensembl] (ftp://ftp.ensemblgenomes.org/pub/plants/current/gtf/arabidopsis_thaliana/).
Obrigado, @Emily_Ensembl,, consegui obter o velocyto para calcular a velocidade "global" do RNA com base nesta sequência de referência. No caso de alguém ter esse problema, o arquivo loom precisa ser recalculado com o mesmo arquivo GTF com o utilitário shell para ser compatível com os comandos R.
Um responda:
Ian Sudbery
2018-12-11 16:23:38 UTC
view on stackexchange narkive permalink

No seu caso, eu definitivamente sugeriria seguir o conselho de @Ensembl e usar o Arabidopsis GTF da Ensembl. Mas para referência futura, se um Ensembl GTF não estivesse disponível, você poderia construir algo assim usando a classe gtf do cgat

O módulo cgat e as dependências podem ser instaladas seguindo as instruções aqui. Especificamente, você pode baixar e executar o script de instalação. Isso configurará um ambiente conda dedicado para executar as versões de dependências necessárias. Documentação adicional pode ser encontrada aqui. 

  # download do script de instalação: curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/ master / install.sh # instale a versão de desenvolvimento (recomendado, nenhuma versão de produção ainda): bash install.sh --devel

Isso requer python3 e anaconda para instalar. Você pode precisar executar este script no ambiente conda: 

  # show available environmentconda info --envs # activate by pathource activate / home / user / local / bin / conda-install / envs / cgat-a

Em seguida, execute o seguinte como um script python3. Copie o conteúdo em um arquivo como convert_gtf.py . Em seguida, execute-o no terminal do seu arquivo gtf: 

  python convert_gtf.py genes.gtf > genes_new.gtf  
  #python 3.6.4import sysfrom cgat import GTFfrom cgatcore import iotoolsfor gene in GTF.flat_gene_iterator (GTF.iterator (iotools.open_file (sys.argv [1])), strict = False): gene_start = min ( e.start para e no gene) gene_end = max (e.end para e no gene) se houver (e.feature == "CDS" para e no gene): gene_type = "protein_coding" else: gene_type = "non_coding" para exon no gene: exon_line.setAttribute ("gene_biotype", gene_type) gene_line = GTF.Entry (). fromGTF (gene [0]) gene_line.feature = "gene" gene_line.start = gene_start gene_line.end = gene_end
# Não está claro qual é o transcript_id para uma linha genética. # Tecnicamente falando, os GTFs do Ensembl não seguem o padrão GTF #, pois as entradas do GTF devem incluir transcript_ids, mas # os do conjunto não. if hasattr (gene_line, "gene_name"): gene_line.attributes ["gene_name"] = gene [0] .gene_name else: gene_line.attributes ["gene_name"] = gene [0] .gene_id gene_line.attributes ["gene_biotype"] = gene_type print (str (gene_line) + "\ n") para exon no gene: print (str (exon_line) + "\ n")
Isso é ótimo saber. Visto que isso se refere ao uso de ferramentas em R, seria útil ter instruções para instalar a biblioteca e executar para aqueles não familiarizados com Python.
Eu testei este script e fiz pequenas alterações. Isso é o que funcionou para mim (para ocultar o GTF: não tenho certeza se ele é compatível com o velocyto ainda).
Obrigado pelas correções. As coisas mudaram recentemente e eu perco a noção de em que versão as coisas estão. Acredito que a equipe CGAT está trabalhando para simplificar a instalação. Alguns pontos: não tenho certeza do que você está tentando alcançar com o bit `if hasattr (gene_line," gene_name) `: presumivelmente, se já tiver um atributo` gene_name`, você não precisa definir um ? Observe também que as linhas resultantes não terão um `transcript_id`. Embora isso também seja verdade para os GTFs do Ensembl, significa que esses GTFs não significam a especificação do GTF.
Tive problemas com nomes de genes ausentes e, muitas vezes, TAIR IDs são usados ​​em vez de símbolos de genes. Portanto, isso pode ser específico para Arabidopsis. O mais importante são as correções para o loop para exons, pois essas linhas não foram definidas.
Não é que eu não entenda que você possa querer adicionar um `gene_name` aos atributos, é apenas que atualmente você testa a presença de um` gene_name` e, em seguida, define-o para seu valor já definido se ele já tiver 1. Isso parece redundante.


ⓘ
Estas perguntas e respostas foram traduzidas automaticamente do idioma inglês.O conteúdo original está disponível em stackexchange, que agradecemos pela licença cc by-sa 4.0 sob a qual é distribuído.
Loading...