Una introducción completa a su genoma con SciPy Stack

La bioinformática es un campo interdisciplinario que desarrolla métodos y herramientas de software para analizar y comprender datos biológicos.

Dicho de manera más simple, simplemente puede pensar en ello como ciencia de datos para biología.

Entre los muchos tipos de datos biológicos, los datos genómicos son uno de los más analizados. Especialmente con el rápido avance de las tecnologías de secuenciación de ADN (NGS) de próxima generación, el volumen de datos genómicos ha crecido exponencialmente. Según Stephens, Zachary D et al., la adquisición de datos genómicos está en la escala de exabyte por año.

En esta publicación, demuestro un ejemplo de análisis de un archivo GFF3 para el genoma humano con SciPy Stack. Generic Feature Format Version 3 (GFF3) es el formato de archivo de texto estándar actual para almacenar características genómicas. En particular, en esta publicación aprenderá cómo usar la pila SciPy para responder las siguientes preguntas sobre el genoma humano:

1. ¿Cuánto del genoma está incompleto?
2. ¿Cuántos genes hay en el genoma?
3. ¿Cuánto dura un gen típico?
4. ¿Cómo es la distribución de genes entre los cromosomas?

El último archivo GFF3 para el genoma humano se puede descargar desde aquí. El archivo README que viene en este directorio proporciona una breve descripción de este formato de datos, y aquí se encuentra una especificación más detallada.

Usaremos Pandas, un componente principal de la pila SciPy que proporciona estructuras de datos rápidas, flexibles y expresivas, para manipular y comprender el archivo GFF3.

Configuración

Lo primero es lo primero, configuremos un entorno virtual con la pila SciPy instalada. Este proceso puede llevar mucho tiempo si se construye manualmente desde el código fuente, ya que la pila involucra muchos paquetes, algunos de los cuales dependen del código FORTRAN o C externo. Aquí, recomiendo usar Miniconda, que hace que la configuración sea muy fácil.

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

El indicador -b en la línea bash le dice que se ejecute en modo por lotes. Después de usar los comandos anteriores para instalar Miniconda correctamente, inicie un nuevo entorno virtual para genómica y luego instale la pila SciPy.

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

Tenga en cuenta que solo hemos especificado los 3 paquetes que vamos a usar en esta publicación. Si desea que todos los paquetes aparezcan en la pila de SciPy, simplemente agréguelos al final del comando de conda create .

Si no está seguro del nombre exacto de un paquete, pruebe conda search . Activemos el entorno virtual e iniciemos IPython.

 source activate venv/ ipython

IPython es un reemplazo significativamente más poderoso de la interfaz del intérprete de Python predeterminado, por lo que cualquier cosa que solía hacer en el intérprete de Python predeterminado también se puede hacer en IPython. Recomiendo encarecidamente a todos los programadores de Python, que aún no han usado IPython, que lo prueben.

Descargar el archivo de anotaciones

Con nuestra configuración ahora completa, descarguemos el archivo de anotación del genoma humano en formato GFF3.

Se trata de 37 MB, un archivo muy pequeño en comparación con el contenido de información de un genoma humano, que es de unos 3 GB en texto plano. Esto se debe a que el archivo GFF3 solo contiene la anotación de las secuencias, mientras que los datos de la secuencia generalmente se almacenan en otro formato de archivo llamado FASTA. Si está interesado, puede descargar FASTA aquí, pero no usaremos los datos de secuencia en este tutorial.

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

El prefijo ! le dice a IPython que este es un comando de shell en lugar de un comando de Python. Sin embargo, IPython también puede procesar algunos comandos de shell de uso frecuente como ls , pwd , rm , mkdir , rmdir , ¡incluso sin un prefijo ! .

Echando un vistazo al encabezado del archivo GFF, verá muchas líneas de metadatos/pragmas/directivas que comienzan con ## o #! .

De acuerdo con README, ## significa que los metadatos son estables, mientras que #! significa que es experimental.

Más adelante también verá ### , que es otra directiva con un significado aún más sutil basado en la especificación.

Se supone que los comentarios legibles por humanos deben estar después de un solo # . Para simplificar, trataremos todas las líneas que comiencen con # como comentarios y simplemente las ignoraremos durante nuestro análisis.

 ##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

La primera línea indica que la versión del formato GFF utilizada en este archivo es 3.

A continuación hay resúmenes de todas las regiones de secuencia. Como veremos más adelante, dicha información también se puede encontrar en la parte del cuerpo del archivo.

Las líneas que comienzan con #! muestra información sobre la construcción particular del genoma, GRCh38.p7, a la que se aplica este archivo de anotaciones.

Genome Reference Consortium (GCR) es un consorcio internacional que supervisa las actualizaciones y mejoras de varios ensamblajes de genoma de referencia, incluidos los de humanos, ratones, peces cebra y pollos.

Explorando este archivo, aquí están las primeras líneas de anotación.

 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 ...

Las columnas son seqid , fuente , tipo , inicio , fin , puntuación , cadena , fase , atributos . Algunos de ellos son muy fáciles de entender. Tome la primera línea como ejemplo:

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

Esta es la anotación del primer cromosoma con un seqid de 1, que comienza desde la primera base hasta la base 24.895.622.

En otras palabras, el primer cromosoma tiene una longitud de unos 25 millones de bases.

Nuestro análisis no necesitará información de las tres columnas con un valor de . (es decir, puntaje, hebra y fase), por lo que simplemente podemos ignorarlos por ahora.

La última columna de atributos dice que el cromosoma 1 también tiene tres nombres de alias, a saber, CM000663.2, chr1 y NC_000001.11. Básicamente, así es como se ve un archivo GFF3, pero no los inspeccionaremos línea por línea, por lo que es hora de cargar todo el archivo en Pandas.

Pandas es una buena opción para trabajar con el formato GFF3 porque es un archivo delimitado por tabuladores, y Pandas tiene muy buen soporte para leer archivos tipo CSV.

Tenga en cuenta que una excepción al formato delimitado por tabuladores es cuando GFF3 contiene ##FASTA .

De acuerdo con la especificación, ##FASTA indica el final de una porción de anotación, a la que seguirá una o más secuencias en formato FASTA (no delimitado por tabuladores). Pero este no es el caso del archivo GFF3 que vamos a analizar.

 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)

La última línea anterior carga todo el archivo pandas.read_csv .

Dado que no es un archivo CSV estándar, debemos personalizar un poco la llamada.

Primero, informamos a Pandas de la falta de disponibilidad de la información del encabezado en el GFF3 con header=None , y luego especificamos el nombre exacto para cada columna con names=col_names .

Si no se especifica el argumento de names , Pandas utilizará números incrementales como nombres para cada columna.

sep='\t' le dice a Pandas que las columnas están separadas por tabuladores en lugar de por comas. El valor de sep= en realidad puede ser una expresión regular (regex). Esto se vuelve útil si el archivo en cuestión usa diferentes separadores para cada columna (oh, sí, eso sucede). comment='#' significa que las líneas que comienzan con # se consideran comentarios y se ignorarán.

compression='gzip' le dice a Pandas que el archivo de entrada está comprimido con gzip.

Además, pandas.read_csv tiene un amplio conjunto de parámetros que permite leer diferentes tipos de formatos de archivo similares a CSV.

El tipo del valor devuelto es un DataFrame , que es la estructura de datos más importante en Pandas, utilizada para representar datos 2D.

Pandas también tiene una estructura de datos de Series y Panel para datos 1D y 3D, respectivamente. Consulte la documentación para obtener una introducción a las estructuras de datos de Pandas.

Echemos un vistazo a las primeras entradas con el 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

La salida está bien formateada en un formato tabular con cadenas más largas en la columna de atributos reemplazada parcialmente con ... .

Puede configurar Pandas para que no omita cadenas largas con pd.set_option('display.max_colwidth', -1) . Además, Pandas tiene muchas opciones que se pueden personalizar.

A continuación, obtengamos información básica sobre este marco de datos con el 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

Esto muestra que el GFF3 tiene 2 601 848 líneas anotadas y cada línea tiene nueve columnas.

Para cada columna, también muestra sus tipos de datos.

Ese inicio y final son de tipo int64, números enteros que representan posiciones en el genoma.

Las otras columnas son todas de tipo object , lo que probablemente significa que sus valores consisten en una mezcla de enteros, flotantes y cadenas.

El tamaño de toda la información es de aproximadamente 178,7+ MB almacenados en la memoria. Este resulta ser más compacto que el archivo sin comprimir, que tendrá unos 402 MB. A continuación se muestra una verificación rápida.

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

Desde una vista de alto nivel, hemos cargado todo el archivo GFF3 en un objeto DataFrame en Python, y todo nuestro análisis siguiente se basará en este único objeto.

Ahora, veamos de qué se trata la primera columna 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 es una forma de acceder a los datos de la columna desde un marco de datos. Otra forma es df['seqid'] , que es una sintaxis más general, porque si el nombre de la columna es una palabra clave reservada de Python (por ejemplo class ) o contiene un . o carácter de espacio, la primera forma ( df.seqid ) no funcionará.

El resultado muestra que hay 194 secuencias únicas, que incluyen cromosomas 1 a 22, X, Y y ADN mitocondrial (MT), así como otras 169 secuencias.

Los seqids que comienzan con KI y GL son secuencias de ADN, o andamios, en el genoma que no se han ensamblado con éxito en los cromosomas.

Para aquellos que no están familiarizados con la genómica, esto es importante.

Aunque el primer borrador del genoma humano salió hace más de 15 años, el genoma humano actual aún está incompleto. La dificultad para ensamblar estas secuencias se debe en gran medida a las complejas regiones repetitivas del genoma.

A continuación, echemos un vistazo a la columna de origen.

El LÉAME dice que la fuente es un calificador de texto libre destinado a describir el algoritmo o el procedimiento operativo que generó esta característica.

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

Este es un ejemplo del uso del método value_counts , que es extremadamente útil para un conteo rápido de variables categóricas.

Del resultado, podemos ver que hay siete valores posibles para esta columna, y la mayoría de las entradas en el archivo GFF3 provienen de havana, ensembl y ensembl_havana.

Puede obtener más información sobre el significado de estas fuentes y las relaciones entre ellas en esta publicación.

Para simplificar las cosas, nos centraremos en las entradas de las fuentes GRCh38, havana, ensembl y ensembl_havan.a.

¿Cuánto del genoma está incompleto?

La información sobre cada cromosoma completo se encuentra en las entradas de la fuente GRCh38, así que primero filtraremos el resto y asignaremos el resultado filtrado a una nueva variable 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...

Filtrar es fácil en Pandas.

Si inspecciona el valor evaluado de la expresión df.source == 'GRCh38' , es una serie de valores True y False para cada entrada con el mismo índice que df . Pasarlo a df[] solo devolverá aquellas entradas donde sus valores correspondientes sean True.

Hay 194 claves en df[] para las cuales df.source == 'GRCh38' .

Como vimos anteriormente, también hay 194 valores únicos en la columna seqid , lo que significa que cada entrada en gdf corresponde a un seqid en particular.

Luego, seleccionamos al azar 10 entradas con el método de sample para observarlas más de cerca.

Puedes ver que las secuencias no ensambladas son de tipo supercontig mientras que las otras son de cromosoma. Para calcular la fracción del genoma que está incompleta, primero necesitamos saber la longitud del genoma completo, que es la suma de las longitudes de todos los 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

En el fragmento anterior primero, hicimos una copia de gdf con .copy() . De lo contrario, el gdf original es solo una porción de df , y modificarlo directamente resultaría en SettingWithCopyWarning (consulte aquí para obtener más detalles).

Luego calculamos la longitud de cada entrada y la volvemos a agregar a gdf como una nueva columna llamada "longitud". La longitud total resulta ser de unos 3100 millones, y la fracción de secuencias sin ensamblar es de alrededor del 0,37 %.

Así es como funciona el corte en los dos últimos comandos.

Primero, creamos una lista de cadenas que cubre todos los seqids de secuencias bien ensambladas, que son todos cromosomas y mitocondrias. Luego usamos el método isin para filtrar todas las entradas cuyos seqid están en la lista chrs .

Se agrega un signo menos ( - ) al comienzo del índice para invertir la selección, porque en realidad queremos todo lo que no está en la lista (es decir, queremos los que no están ensamblados y comienzan con KI y GL)...

Nota: Dado que las secuencias ensambladas y no ensambladas se distinguen por la columna de tipo, la última línea se puede reescribir alternativamente de la siguiente manera para obtener los mismos resultados.

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

¿Cuántos genes hay?

Aquí nos enfocamos en las entradas de source ensembl, havana y ensembl_havana ya que es donde pertenecen la mayoría de las anotaciones.

 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

El método isin se usa nuevamente para filtrar.

Luego, un recuento rápido de valores muestra que la mayoría de las entradas son exón, secuencia de codificación (CDS) y región no traducida (UTR).

Estos son elementos subgénicos, pero lo que buscamos principalmente es el recuento de genes. Como se muestra, hay 42.470, pero queremos saber más.

Específicamente, ¿cuáles son sus nombres y qué hacen? Para responder a estas preguntas, debemos observar detenidamente la información en la columna 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)

Tienen el formato de una lista separada por punto y coma de pares de etiquetas y valores. La información que más nos interesa es el nombre del gen, la ID del gen y la descripción, y los extraeremos con expresión 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)

Primero, extraemos los nombres de los genes.

En la expresión regular Name=(?P<gene_name>.+?); , +? se usa en lugar de + porque queremos que no sea codicioso y que la búsqueda se detenga en el primer punto y coma; de lo contrario, el resultado coincidirá hasta el último punto y coma.

Además, la expresión regular se compila primero con re.compile en lugar de usarse directamente como en re.search para un mejor rendimiento porque la aplicaremos a miles de cadenas de atributos.

extract_gene_name sirve como una función auxiliar para usar en pd.apply , que es el método que se usa cuando se necesita aplicar una función en cada entrada de un marco de datos o serie.

En este caso particular, queremos extraer el nombre del gen para cada entrada en ndf.attributes y agregar los nombres nuevamente a ndf en una nueva columna llamada gene_name .

Los identificadores de genes y la descripción se extraen de manera similar.

 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)

La expresión regular para RE_GENE_ID es un poco más específica ya que sabemos que cada gene_id debe comenzar con ENSG , donde ENS significa conjunto y G significa gen.

Para las entradas que no tienen ninguna descripción, devolveremos una cadena vacía. Después de extraer todo, ya no usaremos la columna de atributos, así que suéltela para mantener las cosas ordenadas y limpias con el 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...

En la llamada anterior, los attributes indican la columna específica que queremos eliminar.

axis=1 significa que estamos soltando una columna en lugar de una fila ( axis=0 por defecto).

inplace=True significa que la eliminación se realiza en el marco de datos en sí mismo en lugar de devolver una nueva copia con la columna especificada eliminada.

Una mirada rápida muestra que la columna de atributos gene_name .head gene_id y desc .

Por curiosidad, veamos si todos los gene_id y gene_name son ú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,)

Sorprendentemente, el número de nombres de genes es menor que el de ID de genes, lo que indica que algún gene_name debe corresponder a múltiples ID de genes. Averigüemos cuáles son.

 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 las entradas por el valor de gene_name , luego contamos la cantidad de elementos en cada grupo con .count() .

Si inspecciona la salida de ndf.groupby('gene_name').count() , todas las columnas se cuentan para cada grupo, pero la mayoría de ellas tienen los mismos valores.

Tenga en cuenta que los valores de NA no se tendrán en cuenta al contar, así que solo cuente la primera columna, seqid (usamos .ix[:, 0] para asegurarnos de que no haya valores de NA).

Luego ordene los valores de conteo con .sort_values ​​e invierta el orden con .ix[::-1] .

Como resultado, el nombre de un gen se puede compartir con hasta siete identificaciones de genes.

 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]

Una mirada más cercana a todos los genes SCARNA20 muestra que, de hecho, todos son diferentes.

Si bien comparten el mismo nombre, están ubicados en diferentes posiciones del genoma.

Sus descripciones, sin embargo, no parecen muy útiles para distinguirlos.

El punto aquí es saber que los nombres de genes no son únicos para todas las identificaciones de genes, y aproximadamente el 0,15% de los nombres son compartidos por múltiples genes.

¿Cuánto dura un gen típico?

Similar a lo que hicimos cuando tratábamos de comprender lo incompleto del genoma, podemos agregar fácilmente una columna 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 algunas estadísticas simples basadas en los valores de longitud:

• La longitud media de un gen es de unas 36.000 bases.

• La longitud media de un gen es de unas 5200 bases.

• Las longitudes mínima y máxima de los genes son de aproximadamente ocho y 2,3 millones de bases, respectivamente.

Debido a que la media es mucho mayor que la mediana, implica que la distribución de tallas está sesgada hacia la derecha. Para tener una apariencia más concreta, tracemos la distribución.

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

Pandas proporciona una interfaz simple para matplotlib para que el trazado sea muy útil con DataFrames o series.

En este caso, dice que queremos un gráfico de histograma ( kind='hist' ) con 50 contenedores, y que el eje y esté en una escala logarítmica ( logy=True ).

Del histograma, podemos ver que la mayoría de los genes están dentro del primer contenedor. Sin embargo, algunas longitudes de genes pueden tener más de dos millones de bases. Descubramos cuáles son:

 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 puede ver, el gen más largo se llama CNTNAP2, que es la abreviatura de proteína asociada a la contactina 2. Según su página de wikipedia,

Este gen abarca casi el 1,6% del cromosoma 7 y es uno de los genes más grandes del genoma humano.

¡Por supuesto! Lo acabamos de comprobar nosotros mismos. Por el contrario, ¿qué pasa con los genes más pequeños? Resulta que pueden ser tan cortos como ocho 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

Las longitudes de los dos casos extremos están separadas por cinco órdenes de magnitud (2,3 millones frente a 8), lo que es enorme y puede ser un indicador del nivel de diversidad de la vida.

Un solo gen se puede traducir a muchas proteínas diferentes a través de un proceso llamado empalme alternativo, algo que no hemos explorado. Dicha información también se encuentra dentro del archivo GFF3, pero fuera del alcance de esta publicación.

Distribución de genes entre cromosomas

Lo último que me gustaría discutir es la distribución de genes entre los cromosomas, que también sirve como ejemplo para presentar el método .merge para combinar dos tramas de datos. Intuitivamente, los cromosomas más largos probablemente albergan más genes. Veamos si eso es cierto.

 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

Tomamos prestada la variable chrs de la sección anterior y la usamos para filtrar las secuencias no ensambladas. Según el resultado, el cromosoma 1 más grande tiene la mayor cantidad de genes. Si bien el cromosoma Y tiene la menor cantidad de genes, no es el cromosoma más pequeño.

Tenga en cuenta que parece que no hay genes en la mitocondria (MT), lo cual no es cierto.

Un poco más de filtrado en el primer DataFrame df devuelto por pd.read_csv muestra que todos los genes MT provienen de la fuente insdc (que se filtraron antes al generar edf donde solo consideramos fuentes de havana, ensembl o 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. Alrededor del 0,37% del genoma humano aún está incompleto a pesar de que el primer borrador salió hace más de 15 años.
2. Hay alrededor de 42,000 genes en el genoma humano según este archivo GFF3 en particular que usamos.
3. La longitud de un gen puede variar desde unas pocas docenas hasta más de dos millones de bases.
4. Los genes no están distribuidos uniformemente entre los cromosomas. En general, cuanto más grande es el cromosoma, más genes alberga, pero para un subconjunto de cromosomas, la correlación puede ser negativa.

El archivo GFF3 es muy rico en información de anotaciones, y solo hemos arañado la superficie. Si está interesado en una mayor exploración, aquí hay algunas preguntas con las que puede jugar:

1. ¿Cuántas transcripciones tiene típicamente un gen? ¿Qué porcentaje de genes tiene más de una transcripción?
2. ¿Cuántas isoformas tiene típicamente una transcripción?
3. ¿Cuántos exones, CDS y UTR tiene típicamente una transcripción? ¿Qué tamaños son?
4. ¿Es posible categorizar los genes según su función como se describe en la columna de descripción?