Mapeo de secuencias de ARNm de longitud completa: el blog de Dave Tang

He usado BLAT para alinear secuencias de ARNm de longitud completa hace mucho tiempo. Dado que BLAT ha estado fuera durante más de 20 años, me preguntaba qué herramienta de alineación moderna podría usar como reemplazo. minimapa2 me vino a la mente y en esta publicación lo uso para mapear algunas secuencias de transcripción conocidas en el genoma y comparar los resultados del mapeo con las coordenadas informadas.

Usaré Docker por el bien de la reproducibilidad; específicamente estoy usando una imagen de propósito general que he preparado previamente. Él Dockerfile proporciona todos los detalles sobre la imagen. Si desea seguir la publicación y no desea utilizar Docker, las distribuciones de Linux más actualizadas debería trabaja bien.

# pull image
docker pull davetang/build:1.2.1

# start a container
docker run --rm -it -v $(pwd):/work -u $(stat -c "%u:%g" $) davetang/build:1.2.1 /bin/bash

Si ves «¡No tengo nombre!» en el aviso después de iniciar un contenedor Docker, simplemente ignórelo. No hay entrada en /etc/contraseña para su ID de usuario y grupo. Queremos ejecutar el contenedor con el mismo ID de usuario y grupo que el sistema operativo host, de modo que tengamos permiso para leer y escribir cualquier archivo generado dentro del contenedor Docker.

Ahora, para comenzar, necesitamos el genoma hg38 para mapear nuestras transcripciones y podemos descargarlo de UCSC.

wget -c https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz

Usaremos GENCODE secuencias de transcripción.

wget -c https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_41/gencode.v41.transcripts.fa.gz

Descargaremos y extraeremos minimapa2 de GitHub.

wget https://github.com/lh3/minimap2/releases/download/v2.24/minimap2-2.24_x64-linux.tar.bz2
tar -xjf minimap2-2.24_x64-linux.tar.bz2

Usaremos herramientas de cama para convertir el archivo SAM generado por minimap2 en un archivo BED.

wget https://github.com/arq5x/bedtools2/releases/download/v2.30.0/bedtools.static.binary -O bedtools
chmod 755 bedtools

también necesitaremos SAMherramientas para realizar algún filtrado en el archivo SAM.

wget https://github.com/samtools/samtools/releases/download/1.15.1/samtools-1.15.1.tar.bz2
tar -xjf samtools-1.15.1.tar.bz2
mkdir tools
my_prefix=$(pwd)/tools
cd samtools-1.15.1
./configure --prefix=$
make && make install
cd ..

Ahora que tenemos todos los archivos requeridos, ejecutaremos minimap2 usando un ajuste preestablecido específico para alineaciones empalmadas largas. Él -a La opción especifica la salida SAM y el -t La opción especifica el número de subprocesos a utilizar.

minimap2-2.24_x64-linux/minimap2 -ax splice -t 8 hg38.fa.gz gencode.v41.transcripts.fa.gz > gencode.v41.transcripts.sam

Convertiremos la salida SAM a CAMA12 formato. Queremos este formato porque queremos comparar fácilmente las alineaciones del minimapa2 con las alineaciones informadas en el Explorador del genoma de la UCSC.

Minimapa2 devuelve alineaciones secundarias, por lo que usaremos samtools para mantener solo la/s alineación/es primaria/s.

tools/bin/samtools view -F 256 -b gencode.v41.transcripts.sam | ./bedtools bamtobed -splitD -bed12 -i - | gzip > gencode.v41.transcripts.minimap2.bed.gz

Para descargar el archivo BED de las transcripciones de GENCODE en el navegador del genoma de la UCSC, utilice el Explorador de tablas. Seleccione lo siguiente en el Explorador de tablas:

  • clado: mamífero
  • genoma: Humano
  • asamblea: Dic. 2013 (GRCh38/hg38)
  • grupo: Genes y predicciones de genes
  • pista: GENCODE V41
  • mesa: gen conocido
  • región: genoma
  • formato de salida: BED – datos extensibles del navegador
  • Nombre del archivo de salida: gencode.v41.transcripts.bed.gz
  • tipo de archivo devuelto: gzip comprimido

Luego haga clic obtener salida y obtener CAMA. (Puede descargar el archivo BED de mi servidor si no tiene ganas de usar el Navegador de tablas).

Para comparar los resultados, debemos hacer coincidir la entrada BED para la misma transcripción y comparar las siguientes columnas:

  • cromo
  • chromStart
  • cromoEnd
  • número de bloques
  • tamaños de bloque
  • bloqueEmpieza

Por ejemplo, si examinamos la alineación de minimap2 para la transcripción ENST00000413465.6, podemos ver que se asigna a chr17 comenzando en la posición 7661778 y terminando en 7676594. La alineación contiene 7 «bloques», que corresponden a los exones. Los tamaños de exón son 236, 110, 113, 184, 279, 22 y 74 y la posición inicial de los exones es 0, 12402, 13080, 13274, 14215, 14603 y 14742. Para obtener las coordenadas genómicas de los exones, añade el bloqueEmpieza hacia chromStart para comenzar y agregar el bloqueEmpieza, tamaños de bloquey chromStart para llegar al final.

zcat gencode.v41.transcripts.minimap2.bed.gz | grep ENST00000413465.6 | cut -f1-3,10-12
chr17   7661778 7676594 7       236,110,113,184,279,22,74       0,12402,13080,13274,14215,14603,14742

Aquí está la alineación informada del Navegador del genoma de UCSC, que es idéntica a la alineación informada por minimap2 (aparte de la coma final al final de bloqueEmpieza y tamaños de bloque).

zcat gencode.v41.transcripts.bed.gz | grep ENST00000413465.6 | cut -f1-3,10-12
chr17   7661778 7676594 7       236,110,113,184,279,22,74,      0,12402,13080,13274,14215,14603,14742,

Iba a escribir un script de Perl para comparar todas las alineaciones, pero como la mayoría de la gente ya no usa Perl, escribí un script de Python en su lugar. Soy un aficionado a Python, así que avíseme si cometí algún error o cómo puedo mejorar mi código. El siguiente código comparará todas las alineaciones que no sean multimapping entre minimap2 y las reportadas en el UCSC Genome Browser.

#!/usr/bin/env python3

import pandas as pd

# BED12 format
bed_col_names = [
   'chrom',
   'chrom_start',
   'chrom_end',
   'name',
   'score',
   'strand',
   'thick_start',
   'thick_end',
   'item_rgb',
   'block_count',
   'block_sizes',
   'block_starts'
]

# tab-delimited files are loaded with read_csv but with sep = "t"
minimap2 = pd.read_csv("gencode.v41.transcripts.minimap2.bed.gz", sep = "t", names = bed_col_names)
ucsc = pd.read_csv("gencode.v41.transcripts.bed.gz", sep = "t", names = bed_col_names)

# remove some columns that will not be used
# whether to drop labels from the index (axis = 0) or columns (axis = 1)
minimap2 = minimap2.drop(['score', 'thick_start', 'thick_end', 'item_rgb'], axis = 1)
ucsc = ucsc.drop(['score', 'thick_start', 'thick_end', 'item_rgb'], axis = 1)

# keep only the transcript ID as the transcript identifier in the minimap2 results
minimap2["name"] = minimap2["name"].str.split("|").str[0]

# set index column, like setting row names
minimap2 = minimap2.set_index(minimap2['name'])
ucsc = ucsc.set_index(ucsc['name'])

# remove the trailing commas in block sizes and starts
ucsc["block_sizes"] = ucsc["block_sizes"].str.replace(",$", "", regex = True)
ucsc["block_starts"] = ucsc["block_starts"].str.replace(",$", "", regex = True)

# transcripts that map to more than one loci
mm_multimappers = minimap2["name"].duplicated()
print("Number of transcripts that minimap2 mapped to one or more loci:")
print(mm_multimappers.value_counts(), "n")
ucsc_multimappers = ucsc["name"].duplicated()
print("Number of transcripts that UCSC mapped to one or more loci:")
print(ucsc_multimappers.value_counts(), "n")

# remove multimappers
# use keep = False to store duplicated transcripts including the first occurrence
mm_multimappers = minimap2["name"].duplicated(keep = False)
ucsc_multimappers = ucsc["name"].duplicated(keep = False)
# the tilde switches the booleans, like !
minimap2 = minimap2[~mm_multimappers]
ucsc = ucsc[~ucsc_multimappers]

# test case using TP53
tp53 = 'ENST00000413465.6'
tp53_minimap2 = minimap2[minimap2['name'] == tp53]
tp53_ucsc = ucsc[ucsc['name'] == tp53]
print("Comparing", tp53, "between minimap2 and UCSC")
print(tp53_minimap2 == tp53_ucsc, "n")

# reorder rows of UCSC to match minimap2
transcripts = minimap2['name']
ucsc = ucsc.reindex(transcripts)
print("Number of transcripts that were mapped identically between minimap2 and UCSC:")
mm_vs_ucsc = (minimap2 == ucsc)
# https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.all.html
all_cols_true = mm_vs_ucsc.all(axis = 'columns')
print(all_cols_true.value_counts(), "n")

# sanity check by manually checking some results
identical_cases = all_cols_true[all_cols_true].head()
nonidentical_cases = all_cols_true[all_cols_true == False].head()

print("Some transcripts that were mapped identically between minimap2 and UCSC:")
print(minimap2[minimap2['name'].isin(identical_cases.index)])
print(ucsc[ucsc['name'].isin(identical_cases.index)], "n")

print("Some transcripts that were mapped differently between minimap2 and UCSC:")
print(minimap2[minimap2['name'].isin(nonidentical_cases.index)])
print(ucsc[ucsc['name'].isin(nonidentical_cases.index)],"n")

quit()

La salida inicial indicará el número de multimappers.

Number of transcripts that minimap2 mapped to one or more loci:
False    249612
True       2429
Name: name, dtype: int64

Number of transcripts that UCSC mapped to one or more loci:
False    272178
True        174
Name: name, dtype: int64

El siguiente conjunto de resultados mostrará resultados que comparan ENST00000413465.6, que vimos anteriormente que era idéntico.

Comparing ENST00000413465.6 between minimap2 and UCSC
                   chrom  chrom_start  chrom_end  name  strand  block_count  block_sizes  block_starts
name
ENST00000413465.6   True         True       True  True    True         True         True          True

El recuento de todas las comparaciones se muestra a continuación e indica que alrededor del 91 % de las alineaciones eran exactamente iguales.

Number of transcripts that were mapped identically between minimap2 and UCSC:
True     225261
False     22179
dtype: int64

Finalmente, el último conjunto de resultados muestra algunos ejemplos de transcripciones mapeadas de manera idéntica y diferente.

Some transcripts that were mapped identically between minimap2 and UCSC:
                  chrom  chrom_start  chrom_end               name strand  block_count   block_sizes block_starts
name
ENST00000456328.2  chr1        11868      14409  ENST00000456328.2      +            3  359,109,1189   0,744,1352
ENST00000473358.1  chr1        29553      31097  ENST00000473358.1      +            3   486,104,122  0,1010,1422
ENST00000461467.1  chr1        35244      36073  ENST00000461467.1      -            2       237,353        0,476
ENST00000606857.1  chr1        52472      53312  ENST00000606857.1      +            1           840            0
ENST00000642116.1  chr1        57597      64116  ENST00000642116.1      +            3   56,157,1201  0,1102,5318
                  chrom  chrom_start  chrom_end               name strand  block_count   block_sizes block_starts
name
ENST00000456328.2  chr1      11868.0    14409.0  ENST00000456328.2      +          3.0  359,109,1189   0,744,1352
ENST00000473358.1  chr1      29553.0    31097.0  ENST00000473358.1      +          3.0   486,104,122  0,1010,1422
ENST00000461467.1  chr1      35244.0    36073.0  ENST00000461467.1      -          2.0       237,353        0,476
ENST00000606857.1  chr1      52472.0    53312.0  ENST00000606857.1      +          1.0           840            0
ENST00000642116.1  chr1      57597.0    64116.0  ENST00000642116.1      +          3.0   56,157,1201  0,1102,5318

Some transcripts that were mapped differently between minimap2 and UCSC:
                   chrom  chrom_start  chrom_end               name strand  block_count                               block_sizes                                       block_starts
name
ENST00000450305.2   chr1        12009      13670  ENST00000450305.2      +            8                   48,5,43,85,80,2,150,218                    0,50,175,603,965,1046,1215,1443
ENST00000488147.1   chr1        14403      29570  ENST00000488147.1      -           11  100,31,152,159,198,136,137,147,99,154,37  0,604,1392,2203,2454,2829,3202,3511,3864,10334...
ENST00000469289.1  chr19        71873      72718  ENST00000469289.1      +            2                                   401,134                                              0,711
ENST00000607096.1   chr9        30143      30281  ENST00000607096.1      +            1                                       138                                                  0
ENST00000417324.1  chr19        76162      77690  ENST00000417324.1      -            3                               621,205,361                                         0,723,1167
                  chrom  chrom_start  chrom_end               name strand  block_count                              block_sizes                                       block_starts
name
ENST00000450305.2  chr1      12009.0    13670.0  ENST00000450305.2      +          6.0                      48,49,85,78,154,218                            0,169,603,965,1211,1443
ENST00000488147.1  chr1      14403.0    29570.0  ENST00000488147.1      -         11.0  98,34,152,159,198,136,137,147,99,154,37  0,601,1392,2203,2454,2829,3202,3511,3864,10334...
ENST00000469289.1  chr1      30266.0    31109.0  ENST00000469289.1      +          2.0                                  401,134                                              0,709
ENST00000607096.1  chr1      30365.0    30503.0  ENST00000607096.1      +          1.0                                      138                                                  0
ENST00000417324.1  chr1      34553.0    36081.0  ENST00000417324.1      -          3.0                              621,205,361                                         0,723,1167

Resumen

Minimap2 (usando el preajuste de lectura larga empalmado) funciona bastante bien con el mapeo de secuencias de ARNm de longitud completa. El 91 % de las alineaciones comparables fue idéntica entre minimapa2 y las alineaciones informadas en el navegador del genoma de la UCSC.

El trabajo adicional debería desglosar cómo cada alineación es diferente. Minimap2 también podría ejecutarse con diferentes parámetros, ya que la siguiente alineación de minimap2 parece dudosa con dos bloques de 5 y 2 nucleótidos de longitud.

ENST00000450305.2   chr1        12009      13670  ENST00000450305.2      +            8                   48,5,43,85,80,2,150,218                    0,50,175,603,965,1046,1215,1443

También se podría usar un alineador adicional, ya que no se garantiza que la alineación informada en el navegador del genoma de UCSC sea la alineación definitiva.

Finalmente, pasaré más tiempo aprendiendo sobre minimap2 ya que creo que puede usarse como mi nueva herramienta de alineación multipropósito.

Imprimir amigable, PDF y correo electrónico

Fuente del artículo

Deja un comentario