Uma introdução abrangente ao seu genoma com a pilha SciPy

A bioinformática é um campo interdisciplinar que desenvolve métodos e ferramentas de software para analisar e entender dados biológicos.

De forma mais simples, você pode simplesmente pensar nisso como ciência de dados para biologia.

Entre os muitos tipos de dados biológicos, os dados genômicos são um dos mais amplamente analisados. Especialmente com o rápido avanço das tecnologias de sequenciamento de DNA (NGS) de próxima geração, o volume de dados genômicos vem crescendo exponencialmente. De acordo com Stephens, Zachary D et al., a aquisição de dados genômicos está na escala de exabyte por ano.

Neste post, demonstro um exemplo de análise de um arquivo GFF3 para o genoma humano com o SciPy Stack. Formato de Recurso Genérico Versão 3 (GFF3) é o formato de arquivo de texto padrão atual para armazenar recursos genômicos. Em particular, neste post você aprenderá como usar a pilha SciPy para responder as seguintes perguntas sobre o genoma humano:

1. Quanto do genoma está incompleto?
2. Quantos genes existem no genoma?
3. Quanto tempo dura um gene típico?
4. Como é a distribuição de genes entre os cromossomos?

O arquivo GFF3 mais recente para o genoma humano pode ser baixado aqui. O arquivo README que vem neste diretório fornece uma breve descrição desse formato de dados, e uma especificação mais completa é encontrada aqui.

Usaremos o Pandas, um componente importante da pilha SciPy que fornece estruturas de dados rápidas, flexíveis e expressivas, para manipular e entender o arquivo GFF3.

Configuração

Primeiramente, vamos configurar um ambiente virtual com a pilha SciPy instalada. Este processo pode ser demorado se construído a partir da fonte manualmente, pois a pilha envolve muitos pacotes – alguns dos quais dependem de código externo FORTRAN ou C. Aqui, recomendo usar o Miniconda, que facilita bastante a configuração.

 wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh -b

O sinalizador -b na linha bash informa para executar no modo de lote. Depois que os comandos acima forem usados ​​para instalar o Miniconda com sucesso, inicie um novo ambiente virtual para genômica e instale a pilha SciPy.

 mkdir -p genomics cd genomics conda create -p venv ipython matplotlib pandas

Observe que especificamos apenas os 3 pacotes que usaremos neste post. Se você quiser todos os pacotes listados na pilha SciPy, basta anexá-los ao final do comando conda create .

Se você não tiver certeza do nome exato de um pacote, tente conda search . Vamos ativar o ambiente virtual e iniciar o IPython.

 source activate venv/ ipython

O IPython é um substituto significativamente mais poderoso para a interface do interpretador Python padrão, portanto, tudo o que você costumava fazer no interpretador python padrão também pode ser feito no IPython. Eu recomendo a todos os programadores Python, que ainda não usaram o IPython, que experimentem.

Baixe o arquivo de anotação

Com nossa configuração agora concluída, vamos baixar o arquivo de anotação do genoma humano no formato GFF3.

É cerca de 37 MB, um arquivo muito pequeno comparado ao conteúdo de informação de um genoma humano, que é cerca de 3 GB em texto simples. Isso porque o arquivo GFF3 contém apenas a anotação das sequências, enquanto os dados da sequência geralmente são armazenados em outro formato de arquivo chamado FASTA. Se você estiver interessado, você pode baixar o FASTA aqui, mas não usaremos os dados de sequência neste tutorial.

 !wget ftp://ftp.ensembl.org/pub/release-85/gff3/homo_sapiens/Homo_sapiens.GRCh38.85.gff3.gz

O prefixado ! informa ao IPython que este é um comando shell em vez de um comando Python. No entanto, o IPython também pode processar alguns comandos shell usados ​​com frequência, como ls , pwd , rm , mkdir , rmdir mesmo sem prefixo ! .

Dando uma olhada no cabeçalho do arquivo GFF, você verá muitas linhas de metadados/pragmas/diretivas começando com ## ou #! .

De acordo com o README, ## significa que os metadados são estáveis, enquanto #! significa que é experimental.

Mais tarde você também verá ### , que é outra diretiva com um significado ainda mais sutil baseado na especificação.

Comentários legíveis por humanos devem estar após um único # . Para simplificar, trataremos todas as linhas que começam com # como comentários e simplesmente as ignoraremos durante nossa análise.

 ##gff-version 3 ##sequence-region 1 1 248956422 ##sequence-region 10 1 133797422 ##sequence-region 11 1 135086622 ##sequence-region 12 1 133275309 ... ##sequence-region MT 1 16569 ##sequence-region X 1 156040895 ##sequence-region Y 2781480 56887902 #!genome-build GRCh38.p7 #!genome-version GRCh38 #!genome-date 2013-12 #!genome-build-accession NCBI:GCA_000001405.22 #!genebuild-last-updated 2016-06

A primeira linha indica que a versão do formato GFF usada neste arquivo é 3.

A seguir estão os resumos de todas as regiões da sequência. Como veremos mais adiante, tais informações também podem ser encontradas na parte do corpo do arquivo.

As linhas que começam com #! mostra informações sobre a construção específica do genoma, GRCh38.p7, à qual este arquivo de anotação se aplica.

O Genome Reference Consortium (GCR) é um consórcio internacional, que supervisiona atualizações e melhorias para vários conjuntos de genoma de referência, incluindo aqueles para humanos, camundongos, peixes-zebra e frangos.

Examinando este arquivo, aqui estão as primeiras linhas de anotação.

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11 ### 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 1 . biological_region 10678 10687 0.999 + . logic_name=eponine 1 . biological_region 10681 10688 0.999 - . logic_name=eponine ...

As colunas são seqid , source , type , start , end , score , strand , phase , attribute . Alguns deles são muito fáceis de entender. Tome a primeira linha como exemplo:

 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_000001.11

Esta é a anotação do primeiro cromossomo com uma seqid de 1, que começa da primeira base até a 24.895.622ª base.

Em outras palavras, o primeiro cromossomo tem cerca de 25 milhões de bases de comprimento.

Nossa análise não precisará de informações das três colunas com valor . (ou seja, pontuação, cadeia e fase), então podemos simplesmente ignorá-los por enquanto.

A última coluna de atributos diz que o cromossomo 1 também tem três nomes de alias, ou seja, CM000663.2, chr1 e NC_000001.11. É basicamente assim que um arquivo GFF3 se parece, mas não vamos inspecioná-los linha por linha, então é hora de carregar o arquivo inteiro no Pandas.

O Pandas é adequado para lidar com o formato GFF3 porque é um arquivo delimitado por tabulação, e o Pandas tem um suporte muito bom para a leitura de arquivos do tipo CSV.

Observe que uma exceção ao formato delimitado por tabulação é quando o GFF3 contém ##FASTA .

De acordo com a especificação, ##FASTA indica o fim de uma porção de anotação, que será seguida por uma ou mais sequências no formato FASTA (não delimitado por tabulação). Mas este não é o caso do arquivo GFF3 que vamos analisar.

 In [1]: import pandas as pd In [2]: pd.__version__ Out[2]: '0.18.1' In [3]: col_names = ['seqid', 'source', 'type', 'start', 'end', 'score', 'strand', 'phase', 'attributes'] Out[3]: df = pd.read_csv('Homo_sapiens.GRCh38.85.gff3.gz', compression='gzip', sep='\t', comment='#', low_memory=False, header=None, names=col_names)

A última linha acima carrega todo o arquivo pandas.read_csv .

Como não é um arquivo CSV padrão, precisamos personalizar um pouco a chamada.

Primeiro, informamos ao Pandas a indisponibilidade de informações de cabeçalho no GFF3 com header=None e, em seguida, especificamos o nome exato para cada coluna com names=col_names .

Se o argumento de names não for especificado, o Pandas usará números incrementais como nomes para cada coluna.

sep='\t' diz ao Pandas que as colunas são separadas por tabulação em vez de separadas por vírgula. O valor para sep= pode realmente ser uma expressão regular (regex). Isso se torna útil se o arquivo em mãos usar separadores diferentes para cada coluna (oh sim, isso acontece). comment='#' significa que as linhas que começam com # são consideradas comentários e serão ignoradas.

compression='gzip' diz ao Pandas que o arquivo de entrada é compactado com gzip.

Além disso, pandas.read_csv possui um rico conjunto de parâmetros que permite que diferentes tipos de formatos de arquivo semelhantes a CSV sejam lidos.

O tipo do valor retornado é um DataFrame , que é a estrutura de dados mais importante no Pandas, usada para representar dados 2D.

O Pandas também possui uma estrutura de dados Series e Panel para dados 1D e 3D, respectivamente. Consulte a documentação para uma introdução às estruturas de dados do Pandas.

Vamos dar uma olhada nas primeiras entradas com o método .head .

 In [18]: df.head() Out[18]: seqid source type start end score strand phase attributes 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 1 1 . biological_region 10469 11240 1.3e+03 . . external_name=oe %3D 0.79;logic_name=cpg 2 1 . biological_region 10650 10657 0.999 + . logic_name=eponine 3 1 . biological_region 10655 10657 0.999 - . logic_name=eponine 4 1 . biological_region 10678 10687 0.999 + . logic_name=eponine

A saída é bem formatada em um formato tabular com strings mais longas na coluna de atributos parcialmente substituídas por ... .

Você pode definir Pandas para não omitir strings longas com pd.set_option('display.max_colwidth', -1) . Além disso, o Pandas possui muitas opções que podem ser personalizadas.

A seguir, vamos obter algumas informações básicas sobre esse dataframe com o método .info .

 In [20]: df.info() <class 'pandas.core.frame.DataFrame'> RangeIndex: 2601849 entries, 0 to 2601848 Data columns (total 9 columns): seqid object source object type object start int64 end int64 score object strand object phase object attributes object dtypes: int64(2), object(7) memory usage: 178.7+ MB

Isso mostra que o GFF3 tem 2.601.848 linhas anotadas e cada linha tem nove colunas.

Para cada coluna, também mostra seus tipos de dados.

Esse início e fim são do tipo int64, inteiros representando posições no genoma.

As outras colunas são todas do tipo object , o que provavelmente significa que seus valores consistem em uma mistura de inteiros, floats e strings.

O tamanho de todas as informações é de cerca de 178,7+ MB armazenados na memória. Isso acaba sendo mais compacto que o arquivo descompactado, que terá cerca de 402 MB. Uma verificação rápida é mostrada abaixo.

 gunzip -c Homo_sapiens.GRCh38.85.gff3.gz > /tmp/tmp.gff3 && du -s /tmp/tmp.gff3 402M /tmp/tmp.gff3

De uma visão de alto nível, carregamos todo o arquivo GFF3 em um objeto DataFrame em Python e todas as nossas análises a seguir serão baseadas nesse único objeto.

Agora, vamos ver do que se trata a primeira coluna seqid.

 In [29]: df.seqid.unique() Out[29]: array(['1', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '2', '20', '21', '22', '3', '4', '5', '6', '7', '8', '9', 'GL000008.2', 'GL000009.2', 'GL000194.1', 'GL000195.1', ... 'KI270757.1', 'MT', 'X', 'Y'], dtype=object) In [30]: df.seqid.unique().shape Out[30]: (194,)

df.seqid é uma maneira de acessar dados de coluna de um dataframe. Outra maneira é df['seqid'] , que é uma sintaxe mais geral, porque se o nome da coluna for uma palavra-chave reservada do Python (por exemplo class ) ou contiver um arquivo . ou caractere de espaço, a primeira maneira ( df.seqid ) não funcionará.

A saída mostra que existem 194 seqids únicos, que incluem cromossomos 1 a 22, X, Y e DNA mitocondrial (MT), bem como 169 outros seqids.

Os seqids começando com KI e GL são sequências de DNA – ou andaimes – no genoma que não foram montados com sucesso nos cromossomos.

Para aqueles que não estão familiarizados com genômica, isso é importante.

Embora o primeiro esboço do genoma humano tenha saído há mais de 15 anos, o genoma humano atual ainda está incompleto. A dificuldade em montar essas sequências é em grande parte devido a regiões repetitivas complexas no genoma.

Em seguida, vamos dar uma olhada na coluna de origem.

O README diz que a fonte é um qualificador de texto livre destinado a descrever o algoritmo ou procedimento operacional que gerou esse recurso.

 In [66]: df.source.value_counts() Out[66]: havana 1441093 ensembl_havana 745065 ensembl 228212 . 182510 mirbase 4701 GRCh38 194 insdc 74

Este é um exemplo do uso do método value_counts , que é extremamente útil para uma contagem rápida de variáveis ​​categóricas.

A partir do resultado, podemos ver que existem sete valores possíveis para esta coluna, e a maioria das entradas no arquivo GFF3 vem de havana, ensembl e ensembl_havana.

Você pode aprender mais sobre o que essas fontes significam e as relações entre elas neste post.

Para simplificar, vamos nos concentrar nas entradas das fontes GRCh38, havana, ensembl e ensembl_havan.a.

Quanto do genoma está incompleto?

As informações sobre cada cromossomo inteiro estão nas entradas da fonte GRCh38, então vamos primeiro filtrar o resto e atribuir o resultado filtrado a uma nova variável gdf .

 In [70]: gdf = df[df.source == 'GRCh38'] In [87]: gdf.shape Out[87]: (194, 9) In [84]: gdf.sample(10) Out[84]: seqid source type start end score strand phase attributes 2511585 KI270708.1 GRCh38 supercontig 1 127682 . . . ID=supercontig:KI270708.1;Alias=chr1_KI270708v... 2510840 GL000208.1 GRCh38 supercontig 1 92689 . . . ID=supercontig:GL000208.1;Alias=chr5_GL000208v... 990810 17 GRCh38 chromosome 1 83257441 . . . ID=chromosome:17;Alias=CM000679.2,chr17,NC_000... 2511481 KI270373.1 GRCh38 supercontig 1 1451 . . . ID=supercontig:KI270373.1;Alias=chrUn_KI270373... 2511490 KI270384.1 GRCh38 supercontig 1 1658 . . . ID=supercontig:KI270384.1;Alias=chrUn_KI270384... 2080148 6 GRCh38 chromosome 1 170805979 . . . ID=chromosome:6;Alias=CM000668.2,chr6,NC_00000... 2511504 KI270412.1 GRCh38 supercontig 1 1179 . . . ID=supercontig:KI270412.1;Alias=chrUn_KI270412... 1201561 19 GRCh38 chromosome 1 58617616 . . . ID=chromosome:19;Alias=CM000681.2,chr19,NC_000... 2511474 KI270340.1 GRCh38 supercontig 1 1428 . . . ID=supercontig:KI270340.1;Alias=chrUn_KI270340... 2594560 Y GRCh38 chromosome 2781480 56887902 . . . ID=chromosome:Y;Alias=CM000686.2,chrY,NC_00002...

A filtragem é fácil no Pandas.

Se você inspecionar o valor avaliado da expressão df.source == 'GRCh38' , é uma série de valores True e False para cada entrada com o mesmo índice que df . Passar para df[] retornará apenas as entradas em que seus valores correspondentes são True.

Existem 194 chaves em df[] para as quais df.source == 'GRCh38' .

Como vimos anteriormente, também existem 194 valores únicos na coluna seqid , o que significa que cada entrada em gdf corresponde a um seqid específico.

Em seguida, selecionamos aleatoriamente 10 entradas com o método de sample para examinar mais de perto.

Você pode ver que as sequências não montadas são do tipo supercontig enquanto as demais são de cromossomo. Para calcular a fração incompleta do genoma, primeiro precisamos saber o comprimento de todo o genoma, que é a soma dos comprimentos de todos os seqids.

 In [90]: gdf = gdf.copy() In [91]: gdf['length'] = gdf.end - gdf.start + 1 In [93]: gdf.head() Out[93]: seqid source type start end score strand phase attributes length 0 1 GRCh38 chromosome 1 248956422 . . . ID=chromosome:1;Alias=CM000663.2,chr1,NC_00000... 248956421 235068 10 GRCh38 chromosome 1 133797422 . . . ID=chromosome:10;Alias=CM000672.2,chr10,NC_000... 133797421 328938 11 GRCh38 chromosome 1 135086622 . . . ID=chromosome:11;Alias=CM000673.2,chr11,NC_000... 135086621 483370 12 GRCh38 chromosome 1 133275309 . . . ID=chromosome:12;Alias=CM000674.2,chr12,NC_000... 133275308 634486 13 GRCh38 chromosome 1 114364328 . . . ID=chromosome:13;Alias=CM000675.2,chr13,NC_000... 114364327 In [97]: gdf.length.sum() Out[97]: 3096629532 In [99]: chrs = [str(_) for _ in range(1, 23)] + ['X', 'Y', 'MT'] In [101]: gdf[-gdf.seqid.isin(chrs)].length.sum() / gdf.length.sum() Out[101]: 0.0037021917421198327

No trecho acima, primeiro, fizemos uma cópia de gdf com .copy() . Caso contrário, o gdf original é apenas uma fatia de df , e modificá-lo diretamente resultaria em SettingWithCopyWarning (veja aqui para mais detalhes).

Em seguida, calculamos o comprimento de cada entrada e adicionamos de volta ao gdf como uma nova coluna chamada “comprimento”. O comprimento total é de cerca de 3,1 bilhões e a fração de sequências não montadas é de cerca de 0,37%.

Aqui está como o fatiamento funciona nos dois últimos comandos.

Primeiro, criamos uma lista de strings que abrange todos os seqids de sequências bem montadas, que são todos cromossomos e mitocôndrias. Em seguida, usamos o método isin para filtrar todas as entradas cujos seqid estão na lista chrs .

Um sinal de menos ( - ) é adicionado ao início do índice para reverter a seleção, porque na verdade queremos tudo o que não está na lista (ou seja, queremos os desmontados começando com KI e GL)…

Nota: Como as sequências montadas e desmontadas são diferenciadas pela coluna de tipo, a última linha pode ser reescrita como a seguir para obter os mesmos resultados.

 gdf[(gdf['type'] == 'supercontig')].length.sum() / gdf.length.sum()

Quantos genes existem?

Aqui nos concentramos nas entradas da fonte ensembl, havana e ensembl_havana, pois são onde a maioria das entradas de anotação pertence.

 In [109]: edf = df[df.source.isin(['ensembl', 'havana', 'ensembl_havana'])] In [111]: edf.sample(10) Out[111]: seqid source type start end score strand phase attributes 915996 16 havana CDS 27463541 27463592 . - 2 ID=CDS:ENSP00000457449;Parent=transcript:ENST0... 2531429 X havana exon 41196251 41196359 . + . Parent=transcript:ENST00000462850;Name=ENSE000... 1221944 19 ensembl_havana CDS 5641740 5641946 . + 0 ID=CDS:ENSP00000467423;Parent=transcript:ENST0... 243070 10 havana exon 13116267 13116340 . + . Parent=transcript:ENST00000378764;Name=ENSE000... 2413583 8 ensembl_havana exon 144359184 144359423 . + . Parent=transcript:ENST00000530047;Name=ENSE000... 2160496 6 havana exon 111322569 111322678 . - . Parent=transcript:ENST00000434009;Name=ENSE000... 839952 15 havana exon 76227713 76227897 . - . Parent=transcript:ENST00000565910;Name=ENSE000... 957782 16 ensembl_havana exon 67541653 67541782 . + . Parent=transcript:ENST00000379312;Name=ENSE000... 1632979 21 ensembl_havana exon 37840658 37840709 . - . Parent=transcript:ENST00000609713;Name=ENSE000... 1953399 4 havana exon 165464390 165464586 . + . Parent=transcript:ENST00000511992;Name=ENSE000... In [123]: edf.type.value_counts() Out[123]: exon 1180596 CDS 704604 five_prime_UTR 142387 three_prime_UTR 133938 transcript 96375 gene 42470 processed_transcript 28228 ... Name: type, dtype: int64

O método isin é usado novamente para filtragem.

Em seguida, uma contagem rápida de valores mostra que a maioria das entradas são exon, sequência de codificação (CDS) e região não traduzida (UTR).

Estes são elementos de subgenes, mas estamos procurando principalmente a contagem de genes. Como mostrado, são 42.470, mas queremos saber mais.

Especificamente, quais são seus nomes e o que eles fazem? Para responder a essas perguntas, precisamos examinar atentamente as informações na coluna de atributos.

 In [127]: ndf = edf[edf.type == 'gene'] In [173]: ndf = ndf.copy() In [133]: ndf.sample(10).attributes.values Out[133]: array(['ID=gene:ENSG00000228611;Name=HNF4GP1;biotype=processed_pseudogene;description=hepatocyte nuclear factor 4 gamma pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:35417];gene_id=ENSG00000228611;havana_gene=OTTHUMG00000016986;havana_version=2;logic_name=havana;version=2', 'ID=gene:ENSG00000177189;Name=RPS6KA3;biotype=protein_coding;description=ribosomal protein S6 kinase A3 [Source:HGNC Symbol%3BAcc:HGNC:10432];gene_id=ENSG00000177189;havana_gene=OTTHUMG00000021231;havana_version=5;logic_name=ensembl_havana_gene;version=12', 'ID=gene:ENSG00000231748;Name=RP11-227H15.5;biotype=antisense;gene_id=ENSG00000231748;havana_gene=OTTHUMG00000018373;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000227426;Name=VN1R33P;biotype=unitary_pseudogene;description=vomeronasal 1 receptor 33 pseudogene [Source:HGNC Symbol%3BAcc:HGNC:37353];gene_id=ENSG00000227426;havana_gene=OTTHUMG00000154474;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000087250;Name=MT3;biotype=protein_coding;description=metallothionein 3 [Source:HGNC Symbol%3BAcc:HGNC:7408];gene_id=ENSG00000087250;havana_gene=OTTHUMG00000133282;havana_version=3;logic_name=ensembl_havana_gene;version=8', 'ID=gene:ENSG00000177108;Name=ZDHHC22;biotype=protein_coding;description=zinc finger DHHC-type containing 22 [Source:HGNC Symbol%3BAcc:HGNC:20106];gene_id=ENSG00000177108;havana_gene=OTTHUMG00000171575;havana_version=3;logic_name=ensembl_havana_gene;version=5', 'ID=gene:ENSG00000249784;Name=SCARNA22;biotype=scaRNA;description=small Cajal body-specific RNA 22 [Source:HGNC Symbol%3BAcc:HGNC:32580];gene_id=ENSG00000249784;logic_name=ncrna;version=1', 'ID=gene:ENSG00000079101;Name=CLUL1;biotype=protein_coding;description=clusterin like 1 [Source:HGNC Symbol%3BAcc:HGNC:2096];gene_id=ENSG00000079101;havana_gene=OTTHUMG00000178252;havana_version=7;logic_name=ensembl_havana_gene;version=16', 'ID=gene:ENSG00000229224;Name=AC105398.3;biotype=antisense;gene_id=ENSG00000229224;havana_gene=OTTHUMG00000152025;havana_version=1;logic_name=havana;version=1', 'ID=gene:ENSG00000255552;Name=LY6G6E;biotype=protein_coding;description=lymphocyte antigen 6 complex%2C locus G6E (pseudogene) [Source:HGNC Symbol%3BAcc:HGNC:13934];gene_id=ENSG00000255552;havana_gene=OTTHUMG00000166419;havana_version=1;logic_name=ensembl_havana_gene;version=7'], dtype=object)

Eles são formatados como uma lista separada por ponto e vírgula de pares de valor de tag. A informação que mais nos interessa é o nome do gene, ID do gene e descrição, e vamos extraí-los com expressão regular (regex).

 import re RE_GENE_NAME = re.compile(r'Name=(?P<gene_name>.+?);') def extract_gene_name(attributes_str): res = RE_GENE_NAME.search(attributes_str) return res.group('gene_name') ndf['gene_name'] = ndf.attributes.apply(extract_gene_name)

Primeiro, extraímos os nomes dos genes.

Na regex Name=(?P<gene_name>.+?); , +? é usado em vez de + porque queremos que não seja ganancioso e deixe a pesquisa parar no primeiro ponto e vírgula; caso contrário, o resultado corresponderá ao último ponto e vírgula.

Além disso, o regex é compilado primeiro com re.compile em vez de ser usado diretamente como em re.search para melhor desempenho, porque o aplicaremos a milhares de strings de atributos.

extract_gene_name serve como uma função auxiliar a ser usada em pd.apply , que é o método a ser usado quando uma função precisa ser aplicada em cada entrada de um dataframe ou série.

Nesse caso específico, queremos extrair o nome do gene para cada entrada em ndf.attributes e adicionar os nomes de volta ao ndf em uma nova coluna chamada gene_name .

Os IDs e a descrição dos genes são extraídos de maneira semelhante.

 RE_GENE_ID = re.compile(r'gene_id=(?P<gene_id>ENSG.+?);') def extract_gene_id(attributes_str): res = RE_GENE_ID.search(attributes_str) return res.group('gene_id') ndf['gene_id'] = ndf.attributes.apply(extract_gene_id) RE_DESC = re.compile('description=(?P<desc>.+?);') def extract_description(attributes_str): res = RE_DESC.search(attributes_str) if res is None: return '' else: return res.group('desc') ndf['desc'] = ndf.attributes.apply(extract_description)

A regex para RE_GENE_ID é um pouco mais específica, pois sabemos que todo gene_id deve começar com ENSG , onde ENS significa ensembl e G significa gene.

Para entradas que não possuem descrição, retornaremos uma string vazia. Depois que tudo for extraído, não usaremos mais a coluna de atributos, então vamos soltá-la para manter as coisas limpas com o método .drop :

 In [224]: ndf.drop('attributes', axis=1, inplace=True) In [225]: ndf.head() Out[225]: seqid source type start end score strand phase gene_id gene_name desc 16 1 havana gene 11869 14409 . + . ENSG00000223972 DDX11L1 DEAD/H-box helicase 11 like 1 [Source:HGNC Sym... 28 1 havana gene 14404 29570 . - . ENSG00000227232 WASH7P WAS protein family homolog 7 pseudogene [Sourc... 71 1 havana gene 52473 53312 . + . ENSG00000268020 OR4G4P olfactory receptor family 4 subfamily G member... 74 1 havana gene 62948 63887 . + . ENSG00000240361 OR4G11P olfactory receptor family 4 subfamily G member... 77 1 ensembl_havana gene 69091 70008 . + . ENSG00000186092 OR4F5 olfactory receptor family 4 subfamily F member...

Na chamada acima, attributes indicam a coluna específica que queremos eliminar.

axis=1 significa que estamos descartando uma coluna em vez de uma linha ( axis=0 por padrão).

inplace=True significa que o descarte é operado no próprio DataFrame em vez de retornar uma nova cópia com a coluna especificada descartada.

Uma rápida olhada em .head mostra que a coluna de atributos realmente desapareceu, e três novas colunas: gene_name , gene_id e desc foram adicionadas.

Por curiosidade, vamos ver se todos gene_id e gene_name são únicos:

 In [232]: ndf.shape Out[232]: (42470, 11) In [233]: ndf.gene_id.unique().shape Out[233]: (42470,) In [234]: ndf.gene_name.unique().shape Out[234]: (42387,)

Surpreendentemente, o número de nomes de genes é menor que o de IDs de genes, indicando que alguns gene_name devem corresponder a vários IDs de genes. Vamos descobrir quais são.

 In [243]: count_df = ndf.groupby('gene_name').count().ix[:, 0].sort_values().ix[::-1] In [244]: count_df.head(10) Out[244]: gene_name SCARNA20 7 SCARNA16 6 SCARNA17 5 SCARNA15 4 SCARNA21 4 SCARNA11 4 Clostridiales-1 3 SCARNA4 3 C1QTNF9B-AS1 2 C11orf71 2 Name: seqid, dtype: int64 In [262]: count_df[count_df > 1].shape Out[262]: (63,) In [263]: count_df.shape Out[263]: (42387,) In [264]: count_df[count_df > 1].shape[0] / count_df.shape[0] Out[264]: 0.0014863047632528842

Agrupamos todas as entradas pelo valor de gene_name e contamos o número de itens em cada grupo com .count() .

Se você inspecionar a saída de ndf.groupby('gene_name').count() , todas as colunas serão contadas para cada grupo, mas a maioria delas terá os mesmos valores.

Observe que os valores de NA não serão considerados na contagem, portanto, tome apenas a contagem da primeira coluna, seqid (usamos .ix[:, 0] para garantir que não haja valores de NA).

Em seguida, classifique os valores de contagem com .sort_values ​​e inverta a ordem com .ix[::-1] .

Como resultado, um nome de gene pode ser compartilhado com até sete IDs de gene.

 In [255]: ndf[ndf.gene_name == 'SCARNA20'] Out[255]: seqid source type start end score strand phase gene_id gene_name desc 179399 1 ensembl gene 171768070 171768175 . + . ENSG00000253060 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 201037 1 ensembl gene 204727991 204728106 . + . ENSG00000251861 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 349203 11 ensembl gene 8555016 8555146 . + . ENSG00000252778 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 718520 14 ensembl gene 63479272 63479413 . + . ENSG00000252800 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 837233 15 ensembl gene 75121536 75121666 . - . ENSG00000252722 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1039874 17 ensembl gene 28018770 28018907 . + . ENSG00000251818 SCARNA20 Small Cajal body specific RNA 20 [Source:RFAM%3BAcc:RF00601] 1108215 17 ensembl gene 60231516 60231646 . - . ENSG00000252577 SCARNA20 small Cajal body-specific RNA 20 [Source:HGNC Symbol%3BAcc:HGNC:32578]

Um olhar mais atento a todos os genes SCARNA20 mostra que eles são realmente todos diferentes.

Enquanto eles compartilham o mesmo nome, eles estão localizados em diferentes posições do genoma.

Suas descrições, no entanto, não parecem muito úteis para distingui-los.

O ponto aqui é saber que os nomes dos genes não são exclusivos para todos os IDs dos genes e cerca de 0,15% dos nomes são compartilhados por vários genes.

Quanto tempo dura um gene típico?

Semelhante ao que fizemos quando estávamos tentando entender a incompletude do genoma, podemos facilmente adicionar uma coluna de length a ndf :

 In [277]: ndf['length'] = ndf.end - ndf.start + 1 In [278]: ndf.length.describe() Out[278]: count 4.247000e+04 mean 3.583348e+04 std 9.683485e+04 min 8.000000e+00 25% 8.840000e+02 50% 5.170500e+03 75% 3.055200e+04 max 2.304997e+06 Name: length, dtype: float64

.describe() calcula algumas estatísticas simples com base nos valores de comprimento:

• O comprimento médio de um gene é de cerca de 36.000 bases

• O comprimento médio de um gene é cerca de 5.200 bases de comprimento

• Os comprimentos mínimo e máximo dos genes são cerca de oito e 2,3 milhões de bases, respectivamente.

Como a média é muito maior que a mediana, isso implica que a distribuição do comprimento é assimétrica para a direita. Para ter uma visão mais concreta, vamos plotar a distribuição.

 import matplotlib as plt ndf.length.plot(kind='hist', bins=50, logy=True) plt.show()

O Pandas fornece uma interface simples para o matplotlib para tornar a plotagem muito útil com DataFrames ou séries.

Neste caso, ele diz que queremos um gráfico de histograma ( kind='hist' ) com 50 bins, e que o eixo y esteja em uma escala logarítmica ( logy=True ).

A partir do histograma, podemos ver que a maioria dos genes está dentro do primeiro compartimento. No entanto, alguns comprimentos de genes podem ser superiores a dois milhões de bases. Vamos descobrir quais são:

 In [39]: ndf[ndf.length > 2e6].sort_values('length').ix[::-1] Out[39]: seqid source type start end score strand phase gene_name gene_id desc length 2309345 7 ensembl_havana gene 146116002 148420998 . + . CNTNAP2 ENSG00000174469 contactin associated protein-like 2 [Source:HG... 2304997 2422510 9 ensembl_havana gene 8314246 10612723 . - . PTPRD ENSG00000153707 protein tyrosine phosphatase%2C receptor type ... 2298478 2527169 X ensembl_havana gene 31097677 33339441 . - . DMD ENSG00000198947 dystrophin [Source:HGNC Symbol%3BAcc:HGNC:2928] 2241765 440886 11 ensembl_havana gene 83455012 85627922 . - . DLG2 ENSG00000150672 discs large MAGUK scaffold protein 2 [Source:H... 2172911 2323457 8 ensembl_havana gene 2935353 4994972 . - . CSMD1 ENSG00000183117 CUB and Sushi multiple domains 1 [Source:HGNC ... 2059620 1569914 20 ensembl_havana gene 13995369 16053197 . + . MACROD2 ENSG00000172264 MACRO domain containing 2 [Source:HGNC Symbol%... 2057829

Como você pode ver, o gene mais longo é chamado CNTNAP2, que é a abreviação de contactina associada à proteína 2. De acordo com sua página da wikipedia,

Este gene abrange quase 1,6% do cromossomo 7 e é um dos maiores genes do genoma humano.

De fato! Acabamos de verificar isso por nós mesmos. Em contraste, e quanto aos menores genes? Acontece que eles podem ser tão curtos quanto oito bases.

 In [40]: ndf.sort_values('length').head() Out[40]: seqid source type start end score strand phase gene_name gene_id desc length 682278 14 havana gene 22438547 22438554 . + . TRDD1 ENSG00000223997 T cell receptor delta diversity 1 [Source:HGNC... 8 682282 14 havana gene 22439007 22439015 . + . TRDD2 ENSG00000237235 T cell receptor delta diversity 2 [Source:HGNC... 9 2306836 7 havana gene 142786213 142786224 . + . TRBD1 ENSG00000282431 T cell receptor beta diversity 1 [Source:HGNC ... 12 682286 14 havana gene 22449113 22449125 . + . TRDD3 ENSG00000228985 T cell receptor delta diversity 3 [Source:HGNC... 13 1879625 4 havana gene 10238213 10238235 . - . AC006499.9 ENSG00000271544 23

Os comprimentos dos dois casos extremos estão separados por cinco ordens de magnitude (2,3 milhões vs. 8), o que é enorme e pode ser uma indicação do nível de diversidade da vida.

Um único gene pode ser traduzido para muitas proteínas diferentes por meio de um processo chamado splicing alternativo, algo que não exploramos. Tais informações também estão dentro do arquivo GFF3, mas fora do escopo deste post.

Distribuição de genes entre cromossomos

A última coisa que gostaria de discutir é a distribuição de genes entre cromossomos, que também serve como exemplo para a introdução do método .merge para combinar dois DataFrames. Intuitivamente, cromossomos mais longos provavelmente hospedam mais genes. Vamos ver se isso é verdade.

 In [53]: ndf = ndf[ndf.seqid.isin(chrs)] In [54]: chr_gene_counts = ndf.groupby('seqid').count().ix[:, 0].sort_values().ix[::-1] Out[54]: chr_gene_counts seqid 1 3902 2 2806 11 2561 19 2412 17 2280 3 2204 6 2154 12 2140 7 2106 5 2002 16 1881 X 1852 4 1751 9 1659 8 1628 10 1600 15 1476 14 1449 22 996 20 965 13 872 18 766 21 541 Y 436 Name: source, dtype: int64

Pegamos emprestada a variável chrs da seção anterior e a usamos para filtrar as sequências não montadas. Com base na saída, o maior cromossomo 1, de fato, tem o maior número de genes. Embora o cromossomo Y tenha o menor número de genes, não é o menor cromossomo.

Observe que parece não haver genes na mitocôndria (MT), o que não é verdade.

Um pouco mais de filtragem no primeiro DataFrame df retornado por pd.read_csv mostra que todos os genes MT são de source insdc (que foram filtrados antes ao gerar edf onde consideramos apenas fontes de havana, ensembl ou ensembl_havana).

 In [134]: df[(df.type == 'gene') & (df.seqid == 'MT')] Out[134]: seqid source type start end score strand phase attributes 2514003 MT insdc gene 648 1601 . + . ID=gene:ENSG00000211459;Name=MT-RNR1;biotype=M... 2514009 MT insdc gene 1671 3229 . + . ID=gene:ENSG00000210082;Name=MT-RNR2;biotype=M... 2514016 MT insdc gene 3307 4262 . + . ID=gene:ENSG00000198888;Name=MT-ND1;biotype=pr... 2514029 MT insdc gene 4470 5511 . + . ID=gene:ENSG00000198763;Name=MT-ND2;biotype=pr... 2514048 MT insdc gene 5904 7445 . + . ID=gene:ENSG00000198804;Name=MT-CO1;biotype=pr... 2514058 MT insdc gene 7586 8269 . + . ID=gene:ENSG00000198712;Name=MT-CO2;biotype=pr... 2514065 MT insdc gene 8366 8572 . + . ID=gene:ENSG00000228253;Name=MT-ATP8;biotype=p... 2514069 MT insdc gene 8527 9207 . + . ID=gene:ENSG00000198899;Name=MT-ATP6;biotype=p... 2514073 MT insdc gene 9207 9990 . + . ID=gene:ENSG00000198938;Name=MT-CO3;biotype=pr... 2514080 MT insdc gene 10059 10404 . + . ID=gene:ENSG00000198840;Name=MT-ND3;biotype=pr... 2514087 MT insdc gene 10470 10766 . + . ID=gene:ENSG00000212907;Name=MT-ND4L;biotype=p... 2514091 MT insdc gene 10760 12137 . + . ID=gene:ENSG00000198886;Name=MT-ND4;biotype=pr... 2514104 MT insdc gene 12337 14148 . + . ID=gene:ENSG00000198786;Name=MT-ND5;biotype=pr... 2514108 MT insdc gene 14149 14673 . - . ID=gene:ENSG00000198695;Name=MT-ND6;biotype=pr... 2514115 MT insdc gene 14747 15887 . + . ID=gene:ENSG00000198727;Name=MT-CYB;biotype=pr...

This example also shows how to combine two conditions during filtering with & ; the logical operator for “or” would be | .

Note that the parentheses around each condition are required, and this part of the syntax in Pandas is different from Python, which would have been literal and and or .

Next, let's borrow the gdf DataFrame from the previous section as a source for the length of each chromosome:

 In [61]: gdf = gdf[gdf.seqid.isin(chrs)] In [62]: gdf.drop(['start', 'end', 'score', 'strand', 'phase' ,'attributes'], axis=1, inplace=True) In [63]: gdf.sort_values('length').ix[::-1] Out[63]: seqid source type length 0 1 GRCh38 chromosome 248956422 1364641 2 GRCh38 chromosome 242193529 1705855 3 GRCh38 chromosome 198295559 1864567 4 GRCh38 chromosome 190214555 1964921 5 GRCh38 chromosome 181538259 2080148 6 GRCh38 chromosome 170805979 2196981 7 GRCh38 chromosome 159345973 2514125 X GRCh38 chromosome 156040895 2321361 8 GRCh38 chromosome 145138636 2416560 9 GRCh38 chromosome 138394717 328938 11 GRCh38 chromosome 135086622 235068 10 GRCh38 chromosome 133797422 483370 12 GRCh38 chromosome 133275309 634486 13 GRCh38 chromosome 114364328 674767 14 GRCh38 chromosome 107043718 767312 15 GRCh38 chromosome 101991189 865053 16 GRCh38 chromosome 90338345 990810 17 GRCh38 chromosome 83257441 1155977 18 GRCh38 chromosome 80373285 1559144 20 GRCh38 chromosome 64444167 1201561 19 GRCh38 chromosome 58617616 2594560 Y GRCh38 chromosome 54106423 1647482 22 GRCh38 chromosome 50818468 1616710 21 GRCh38 chromosome 46709983 2513999 MT GRCh38 chromosome 16569

The columns that are not relevant to the analysis are dropped for clarity.

Yes, .drop can also take a list of columns and drop them altogether in one operation.

Note that the row with a seqid of MT is still there; we will get back to it later. The next operation we will perform is merge the two datasets based on the values of seqid.

 In [73]: cdf = chr_gene_counts.to_frame(name='gene_count').reset_index() In [75]: cdf.head(2) Out[75]: seqid gene_count 0 1 3902 1 2 2806 In [78]: merged = gdf.merge(cdf, on='seqid') In [79]: merged Out[79]: seqid source type length gene_count 0 1 GRCh38 chromosome 248956422 3902 1 10 GRCh38 chromosome 133797422 1600 2 11 GRCh38 chromosome 135086622 2561 3 12 GRCh38 chromosome 133275309 2140 4 13 GRCh38 chromosome 114364328 872 5 14 GRCh38 chromosome 107043718 1449 6 15 GRCh38 chromosome 101991189 1476 7 16 GRCh38 chromosome 90338345 1881 8 17 GRCh38 chromosome 83257441 2280 9 18 GRCh38 chromosome 80373285 766 10 19 GRCh38 chromosome 58617616 2412 11 2 GRCh38 chromosome 242193529 2806 12 20 GRCh38 chromosome 64444167 965 13 21 GRCh38 chromosome 46709983 541 14 22 GRCh38 chromosome 50818468 996 15 3 GRCh38 chromosome 198295559 2204 16 4 GRCh38 chromosome 190214555 1751 17 5 GRCh38 chromosome 181538259 2002 18 6 GRCh38 chromosome 170805979 2154 19 7 GRCh38 chromosome 159345973 2106 20 8 GRCh38 chromosome 145138636 1628 21 9 GRCh38 chromosome 138394717 1659 22 X GRCh38 chromosome 156040895 1852 23 Y GRCh38 chromosome 54106423 436

Since chr_gene_counts is still a Series object, which doesn't support a merge operation, we need to convert it to a DataFrame object first with .to_frame .

.reset_index() converts the original index (ie seqid ) into a new column and resets current index as 0-based incremental numbers.

The output from cdf.head(2) shows what it looks like. Next, we used the .merge method to combine the two DataFrame on the seqid column ( on='seqid' ).

After merging gdf and cdf , the MT entry is still missing. This is because, by default, .merge operates an inner join, while left join, right join, or outer join are available by tuning the how parameter.

Please refer to the documentation for more details.

Later, you may find that there is also a related .join method. .merge and .join are similar yet have different APIs.

According to the official documentation says

The related DataFrame.join method, uses merge internally for the index-on-index and index-on-column(s) joins, but joins on indexes by default rather than trying to join on common columns (the default behavior for merge). If you are joining on index, you may wish to use DataFrame.join to save yourself some typing.

Basically, .merge is more general-purpose and used by .join .

Finally, we are ready to calculate the correlation between chromosome length and gene_count .

 In [81]: merged[['length', 'gene_count']].corr() Out[81]: length gene_count length 1.000000 0.728221 gene_count 0.728221 1.000000

By default .corr calculates the Pearson correlation between all pairs of columns in a dataframe.

But we have only a single pair of columns in this case, and the correlation turns out to be positive – 0.73.

In other words, the larger the chromosome, the more likely it is to have more genes. Let's also plot the two columns after sorting the value pairs by length :

 ax = merged[['length', 'gene_count']].sort_values('length').plot(x='length', y='gene_count',) # add some margin to both ends of x axis xlim = ax.get_xlim() margin = xlim[0] * 0.1 ax.set_xlim([xlim[0] - margin, xlim[1] + margin]) # Label each point on the graph for (s, x, y) in merged[['seqid', 'length', 'gene_count']].sort_values('length').values: ax.text(x, y - 100, str(s)) 

As seen in image above, even though it is a positive correlation overall, it does not hold for all chromosomes. In particular, for Chromosome 17, 16, 15, 14, 13, the correlation is actually negative, meaning the number of genes on the chromosome decreases as the chromosome size increases.

Findings and Future Research

That ends our tutorial on the manipulation of an annotation file for human genome in GFF3 format with the SciPy stack. The tools we've mainly used include IPython, Pandas, and matplotlib. During the tutorial, not only have we learned some of the most common and useful operations in Pandas, we also answered some very interesting questions about our genome. In summary:

1. Cerca de 0,37% do genoma humano ainda está incompleto, embora o primeiro rascunho tenha saído há mais de 15 anos.
2. Existem cerca de 42.000 genes no genoma humano com base nesse arquivo GFF3 específico que usamos.
3. O comprimento de um gene pode variar de algumas dezenas a mais de dois milhões de bases.
4. Os genes não são distribuídos uniformemente entre os cromossomos. No geral, quanto maior o cromossomo, mais genes ele hospeda, mas para um subconjunto dos cromossomos, a correlação pode ser negativa.

O arquivo GFF3 é muito rico em informações de anotação e acabamos de arranhar a superfície. Se você estiver interessado em explorar mais, aqui estão algumas perguntas com as quais você pode brincar:

1. Quantas transcrições um gene normalmente tem? Que porcentagem de genes tem mais de 1 transcrição?
2. Quantas isoformas uma transcrição normalmente tem?
3. Quantos éxons, CDS e UTRs uma transcrição normalmente possui? Quais são os tamanhos?
4. É possível categorizar os genes com base em sua função, conforme descrito na coluna de descrição?