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/bashSi 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.gzUsaremos GENCODE secuencias de transcripción.
wget -c https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_41/gencode.v41.transcripts.fa.gzDescargaremos 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.bz2Usaremos 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 bedtoolstambié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.samConvertiremos 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.gzPara 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,14742Aquí 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: int64El 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 TrueEl 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: int64Finalmente, 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,1167Resumen
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,1443Tambié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.
