Cómo hacer una transcripción a un archivo de mapeo de genes

Necesito una transcripción del archivo de mapeo de genes para Salmon. Soy consciente de la anotación bioconductor paquetes que pueden hacer este trabajo. Sin embargo, estaba trabajando en una especie que no tiene la anotación en un formato de paquete (voy a usar Drosphila como ejemplo para esta publicación de blog). Tuve que ir y obtener el archivo gtf e hice un archivo de este tipo desde cero.

Por favor, lea el especificaciones de esos dos formatos de archivo.

Use el comando Unix para hacer una transcripción al archivo de mapeo de genes desde el archivo gtf

Vemos que los tipos de características son bastante diferentes, aunque ambos son archivos de anotaciones para la misma especie. los gtf El archivo está relativamente bien formateado, y podemos hacer una transcripción al archivo de mapeo de genes fácilmente usando la línea de comandos de Unix.

zless -S ~/Downloads/Drosophila_melanogaster.BDGP6.95.gtf.gz | grep -v "#" | awk '$3=="transcript"' | cut -f9 | tr -s ";" " " | awk '' | sort | uniq |  sed 's/"//g' | tee tx2gene_ensemble.tsv| head
#
#
#
#
#
#
#
#
#
#

hmm… ¿por qué la primera línea tiene ambos genes en las dos columnas?… control de cordura:

zless -S ~/Downloads/Drosophila_melanogaster.BDGP6.95.gtf.gz | grep "FBgn0013687" | less -S
## mitochondrion_genome FlyBase gene    14917   19524   .   +   .   gene_id "FBgn0013687"; gene_name "mt:ori"; gene_source "FlyBase"; gene_biotype "pseudogene";
## mitochondrion_genome FlyBase transcript  14917   19524   .   +   .   gene_id "FBgn0013687"; transcript_id "FBgn0013687"; gene_name "mt:ori"; gene_source "FlyBase"; gene_biotype "pseudogene"; transcript_source "FlyBase"; transcript_biotype "pseudogene";
## mitochondrion_genome FlyBase exon    14917   19524   .   +   .   gene_id "FBgn0013687"; transcript_id "FBgn0013687"; exon_number "1"; gene_name "mt:ori"; gene_source "FlyBase"; gene_biotype "pseudogene"; transcript_source "FlyBase"; transcript_biotype "pseudogene"; exon_id "FBgn0013687-E1";

De hecho, está en el archivo gtf original.

Usar gffutilspara hacer una transcripción al archivo de mapeo de genes desde el archivo gff

los gff El archivo no está tan bien definido. Todavía se pueden usar algunos trucos de Unix para obtener el archivo tx2gene.tsv de un archivo gff, pero puede ser bastante incómodo, especialmente para archivos gff de otras especies que no están bien anotadas. En cambio, usemos gffutilsun paquete de python para hacer lo mismo.

Instalar en pc gffutils en terminales:

source activate snakemake
conda install gffutils

Tenga en cuenta que estoy ejecutando python a través de Rsutdio/ Primero lea cómo configurar la ruta de python para reticulate a https://rstudio.github.io/reticulate/articles/versions.html leer más sobre https://cran.r-project.org/web/packages/reticulate/vignettes/versions.html

De alguna manera, tengo que crear un .Rprofile en la misma carpeta de .Rproj archivo con la siguiente línea para usar mi entorno conda de snakemake que es python3:

Sys.setenv(PATH = paste("/anaconda3/envs/snakemake/bin/", Sys.getenv("PATH"), sep=":"))

library(reticulate)


py_discover_config()
#
#
#
#
#
#
#
#
#
#
#
#


import sys
print(sys.version)
#
#
import gffutils
import itertools
import os
os.listdir()
db = gffutils.create_db("GCF_000001215.4_Release_6_plus_ISO1_MT_genomic.gff.gz", ":memory:", force = True,merge_strategy="merge", id_spec=)
list(db.featuretypes())

for mRNA in itertools.islice(db.features_of_type('mRNA'), 10):
        print(mRNA['transcript_id'][0], mRNA['gene'][0])
        
        

#
#
#
#
#
#
#
#
#
#
tx_and_gene=[]
with open("tx2gene_NCBI.tsv", "w") as f:
        for feature in db.all_features():
                transcript = feature.attributes.get('transcript_id', [None])[0]
                gene = feature.attributes.get('gene', [None])[0]
                if gene and transcript and ([transcript, gene] not in tx_and_gene):
                        tx_and_gene.append([transcript, gene])
                        f.write(transcript + "t" + gene + "n")

Estas líneas de códigos no son difíciles de escribir. Se necesita más tiempo para leer la documentación del paquete y comprender cómo usarlo. Un problema con la bioinFORMÁTICA es que hay tantos formatos de archivo diferentes. Para empeorar las cosas, incluso para el formato de archivo gff, muchos archivos no siguen la especificación exacta. Puedes probarlo en http://daler.github.io/gffutils/ejemplos.html.

Fuente del artículo

¿Que te ha parecido?

Deja un comentario