Enciclopedia de uso de Samtools

1. Perspectiva

Las funciones principales del comando de vista son: intercambiar archivos sam con archivos bam y luego realizar diversas operaciones en los archivos bam, como ordenar y extraer datos ( estas operaciones están en un archivo bam, por lo que esta operación no se puede realizar cuando la entrada es un archivo sam; finalmente, los datos ordenados o extraídos se envían en formato bam o sam (predeterminado).

Ventajas de los archivos bam: los archivos bam son archivos binarios y ocupan menos espacio en disco que los archivos de texto sam; las operaciones que utilizan archivos binarios bam son más rápidas.

En el comando de vista, la entrada (-t o -T) y la salida (-h) del encabezado del archivo sam (ID de secuencia) se controlan mediante parámetros separados.

Uso: samtools view[options]<in.bam>|<in.sam>[region1 [...]]

De forma predeterminada, si no, al agregar regiones se generarán todas las regiones.

Opciones:

-b salida BAM

? De forma predeterminada, la salida es un archivo en formato SAM. Este parámetro establece el formato BAM de salida.

¿El título impreso de la salida SAM

? De forma predeterminada, el archivo de salida en formato sam no tiene encabezado. Este parámetro establece la información del encabezado cuando se genera el archivo sam.

-H ¿Imprimir solo encabezado (sin alineación)

? ¿Solo genera el encabezado del archivo?

-¿La entrada de S es SAM

? De forma predeterminada, la entrada es un archivo BAM. Si la entrada es un archivo SAM, es mejor agregar este parámetro; de lo contrario, a veces se producen errores.

-u salida BAM sin comprimir (force -b)

? El uso de este parámetro requiere un parámetro -b, que ahorra tiempo, pero requiere más espacio en disco.

-c no imprime las alineaciones, solo las cuenta y las imprime

? total. ¿Todas las opciones de filtro como '-f ', '-F ' y '-q '? son todos tomados en cuenta.

? ¿La función estadística del trabajo de filtrado

-c solo imprime el recuento de registros coincidentes

-Archivo L? La alineación de salida se superpone con el archivo BED de entrada [vacío]

-¿archivo t? Nombre de referencia y lista de longitud (obligatorio) [vacío]

? Utilice un archivo de lista como entrada para el título.

-¿Archivo T? ¿Archivo de secuencia de referencia (force-S) [vacío]

? Utilice archivos fasta de secuencia como entrada para el encabezado.

-o archivo? Nombre del archivo de salida [stdout]

-F INT? Indicador de filtro, ¿0 significa no establecido [0]

? ¿Saltar alineación con bits en INT[0]

? El número 4 indica que la secuencia no se alinea con la secuencia de referencia.

? El número 8 indica que la secuencia emparejada de esta secuencia no se alinea con la secuencia de referencia.

? Función de filtro. Por ejemplo, F12 solo filtra mapas de doble extremo.

-qINT? Calidad mínima de mapeo [0]

En general, el valor de calidad mínimo para la comparación es 20, lo que indica una comparación única. Puede combinar esto con el parámetro -bF mencionado anteriormente para extraer resultados de comparación específicos.

Ejemplo:

Convertir archivo sam a archivo bam

vista samtools -bS ABC.Sam>abc.bam

Bang, Sam

samtools view -h -o out.sam out.bam

Extraer resultados de comparación con la secuencia de referencia

samtools view - bf4 ABC . .

F.bam

Para extraer los resultados de la comparación donde ambos pares de lecturas emparejadas están alineadas con la secuencia de referencia, solo es necesario utilizar como filtro los dos valores de 4+8 y 12. parámetros.

samtools view-bF 12 ABC . F12.bam

Extrae los resultados de la comparación que no coinciden con la secuencia de referencia.

samtools view-bf4 ABC.bam>abc.f.bam

Extraiga los resultados de la comparación del archivo bam a caffold1 y guárdelo en el formato de archivo sam.

samtools ver ABC . bam scaffold 1 > scaffold 1. sam

Extrae resultados de comparación en scaffold 1 que pueden comparar 30k con 100k.

samtools Ver ABC .bam scaffold 1:30000-100000 $ gt; scaffold 1_30k-100k.sam

Agregue encabezados a archivos sam o bam según el archivo fasta.

Sam tools view-T genome . fasta-h scaffold 1 . Sam & gt; scaffold 1.h.sam

2 Clasificación

Ordenar archivos bam. están ordenados. Algunos programas requieren archivos bam o sam ordenados, como stringtie, por lo que se debe utilizar sort. Cuando busque profundidad, también debe ordenar;;

Uso: Sam tools sort[-n][-m <maxMem>]<in.bam><out.prefix>?

De forma predeterminada, el parámetro de memoria -m es 500, 000, 000 o 500 M (abreviaturas como k, myg no son compatibles). Al procesar datos grandes, establecer un valor mayor puede ahorrar tiempo si hay suficiente memoria.

-n establece el método de clasificación para ordenar por ID de lecturas cortas. De forma predeterminada, se ordena de izquierda a derecha según la secuencia (es decir, el encabezado) y la posición de la secuencia en el archivo fasta.

Ejemplo:

Samtools sort Accept.bamaccept.sort finalmente genera Accept.sort.bam

Fusionar

Combina los dos. o más archivos bam ordenados fusionados en un archivo bam. Los archivos combinados han sido ordenados.

Uso:? Fusionar herramientas Sam[-NR][-h inh . Sam]<out.bam><in 1.bam><in2.bam>[...]

Opciones:- n? Ordenar por nombre de lectura

-r? Agregar etiqueta RG (inferida del nombre del archivo)

- ¿Y tú? Salida BAM sin comprimir

-f? Si está presente, sobrescriba la salida BAM

-1 ?Nivel de compresión 1

-R Street? ¿Fusionar archivos en el área especificada STR [todos]

-h archivo? Copie los encabezados del archivo a & ltout.bam & gt[in1.bam]

Ejemplo:

4. Índice

De forma predeterminada, primero debe archivos Bam. se clasifican antes de poder indexarse. De lo contrario, informará un error.

Después de la indexación, los archivos con sufijo. Se generará Bai para un procesamiento aleatorio rápido. Los archivos Bai son necesarios en muchos casos, especialmente en el caso de alineación de secuencias. Por ejemplo, el comando tview de samtool es necesario; Gbrowse2 también lo necesita al mostrar el gráfico de comparación de lectura. También es necesaria la comparación de monitores IGV.

Uso: samtools index & ltin.bam & gt[out.index]

Ejemplo:

Los dos comandos siguientes tienen el mismo resultado.

$ índice de herramientas abc.sort.bam

índice de herramientas $ Sam ABC. p> p>

Archivos fasta de índice, los archivos de índice generados comienzan con. sufijo fai. Este comando también puede extraer rápidamente una (sub)secuencia de un archivo fasta basado en un archivo de índice.

Uso: samtools faidx & ltin.bam & gt[ [...]]

Indexar archivos del genoma para facilitar la extracción de secuencias.

Ejemplo: $ samtools faidx genome.fasta

Debido al archivo de índice, las subsecuencias en formato fasta se pueden extraer rápidamente del genoma usando el siguiente comando.

$ Sam tools faidx genome .fasta SCF fold _ 10 & gt;scaffold_10.fasta

6 .View

Tview puede mostrar visualmente cómo se relacionan las lecturas con el genoma A modo de comparación, similar al navegador de genoma.

Debe utilizar los comandos de clasificación e indexación mencionados anteriormente de antemano y luego utilizar los siguientes comandos.

Uso: samtools tview & ltaln.bam & gt[reference fasta]

Cuando se muestra el genoma de referencia, la secuencia del genoma de referencia se mostrará en la primera línea; de lo contrario, , la primera línea Que todos estén representados por n.

Presiona G y se te pedirá que ingreses a un sitio determinado en el genoma. El ejemplo "scaffold_10:1000" representa la base número 1000 del primer andamio 1000.

Utilice h (izquierda) j (arriba) k (abajo) l (derecha) para mover la interfaz de visualización. Las letras mayúsculas se mueven más rápido, las minúsculas se mueven más lento.

Utiliza la barra espaciadora para moverte rápidamente hacia la izquierda (similar a L) y la tecla de retroceso para moverte rápidamente hacia la izquierda (similar a H).

Ctrl+H mueve 1kb de distancia base hacia la izquierda; Ctrl+L mueve 1kb de distancia base hacia la derecha.

Se pueden utilizar colores para marcar calidad de comparación, calidad de base, nucleótidos, etc. La calidad básica o calidad comparativa de 30~40 está representada por el blanco; 20~30 es amarillo; 10~20 es verde;

Utilice el punto "." para cambiar la visualización de bases y puntos; utilice el interruptor r para mostrar los nombres leídos, etc.

Hay muchas otras instrucciones, según las específicas? Ver clave.

7.flagstat

proporciona los resultados de la comparación de archivos BAM.

Uso: samtools flagstat & ltin.bam & gt

Ejemplo de $ samtools flagstat bam

Un total de 11945742+0 (lecturas de control de calidad aprobadas + control de calidad fallido). lecturas)

#Total * * *Número de lecturas

0 + 0 copias

7536364 + 0 asignadas (63,09%:-nan%)

p>

#Tasa general de lecturas coincidentes

11945742 + 0 secuenciación emparejada

#¿Cuántas lecturas son lecturas emparejadas?

5972871 + 0 lectura1

#reads1 número de lecturas

5972871 + 0 lectura2

#reads2 número de lecturas

p>

6412042 + 0 emparejamiento correcto (53,68%: nan%)

#Número de lecturas que coinciden exactamente: coincide con la misma secuencia de referencia y la distancia entre las dos lecturas cumple con el establecer umbral.

6899708 + 0 y sus asignaciones propias y emparejadas

#El número de lecturas emparejadas que coinciden con la secuencia de referencia.

636656 + 0 monomórfico (5,33%:-nan%)

#La lectura única que coincide con la secuencia de referencia, más la lectura anterior, es el número total de lecturas coincidentes.

469868 + 0, la relación de pareja se asigna a diferentes chr

#lecturas emparejadas, el número de lecturas en las que dos pares coinciden con dos secuencias de referencia diferentes respectivamente.

243047 + 0, relaciona mapas con diferentes chr (mapQ>=5)

#Igual que arriba, solo compara la calidad>; >8. Depth

Obtiene la profundidad de secuenciación de cada sitio base y la genera en la salida estándar, así que agréguelo a un archivo con un signo mayor que.

Uso: bam 2 profundidad[-r reg][-Q baseQthres][-Q mapq thres][-b in .bed]<in 1. bam>[...]

-r seguido del número de cromosoma (región)

-q: el valor de calidad mínimo de la calidad de secuenciación básica requerida al calcular la profundidad.

-Q: El valor de masa mínimo que debe compararse al calcular la profundidad.

Nota: el índice de samtools; debe completarse antes de la profundidad

Ejemplo

profundidad de samtools aceptar.bam & gt profundidad

9. Otros pedidos

Reencabezado: Reemplaza el encabezado del archivo bam.

$ samtools header & ltin. Sam> Número de lecturas comparadas, número de lecturas no asignadas". La cuarta columna debe ser el número de lecturas en las que un extremo de las lecturas emparejadas puede coincidir con el andamio y el otro extremo no coincide con ningún andamio.

$ samtools idxstats & ltaln.bam & gt

Rmdup: Elimina las lecturas obtenidas por duplicados de PCR, reteniendo solo las lecturas con la mayor calidad de alineación.

Uso:? samtools rmdup [-sS]?

-s para lecturas de un solo extremo. De forma predeterminada, solo se usa para lecturas de extremos emparejados.

-S trata las lecturas de extremos emparejados como lecturas de un solo extremo.

10. Convertir archivos bam a archivos fastq

A veces, necesitamos extraer lecturas que coincidan con la secuencia de referencia y realizar análisis a pequeña escala para facilitar la depuración, etc. En este momento, debe convertir el archivo bam o sam al formato fastq.

Este sitio web proporciona un programa para convertir bam a fastq: http://www.hudsonalpha.org/GSL/information/software/bam2fatq.

$ wget http://www . hudsonalpha .org/GSL/static/software/bam 2 fastq-1.1.0 .

$ tar zxf bam 2 fastq-1.1. 0 . tgz

$ cd bam2fastq-1.1.0

$ make

$ ./ba m2 fastq & lt.bam & gt

11.mpileup

Samtools también tiene un comando muy importante compilar, anteriormente llamado pileup. Este comando se usa para generar archivos bcf y luego usar bcftools para analizar SNP e Indels. Bcftools es el software incluido en samtools y se puede encontrar en la carpeta de instalación de SamTools.

Los parámetros más utilizados son 2:

-f ingresa la secuencia de referencia fasta del archivo de índice; -g genera el formato bcf.

El uso y el ejemplo más simple son los siguientes

Uso: Sam tools mpileup[-EBug][-C capQcoef][-r reg][-f in ][-l list][-M. capMapQ] [-Q minBaseQ][-Q minMapQ]in . bam[in2 . bam[...]]$ Sam herramientas mpileup -f fasta ABC . herramientas mpileup-gSDf genoma . fasta ABC . bcftools view -cvNg->abc.vcf

Cuando compilar no utiliza los parámetros -u o -g, no genera un archivo binario bcf, sino un archivo de texto (salida estándar). El archivo de texto cuenta la alineación de cada sitio base en la secuencia de referencia; cada línea del archivo representa el resultado de la comparación de un sitio base en la secuencia de referencia. Por ejemplo:

Andamio_1?2841 A? 11 ?,,,...,....BHIGDGIJ? Extinción de incendios

Andamio_1?2842 C? 12 ?,$,,...,....^i. cfgeggcff+

andamios_1?2843 G? 11 ?,,...,....FDDDDCD? DD+

Andamios_1?2844 G? 11 ?,,...,... ley? AAAA<AA+

Andamios_1?2845 G? 11 ?,,...,....f 65666166 *

Andamio_1?2846 A? 11 ?,,...,....(1.1111)11*

Andamios_1?2847 A? 11?,,+9acggtgaag. +9ACGGTGAAT. +9ACGGTGAAG. +9ACGGTGAAG, +9acggtgaag. +9ACGGTGAAG. +9ACGGTGAAG. +9ACGGTGAAG. +9ACGGTGAAG. +9ACGGTGAAG? %.+....-..)

Andamio_1?2848 N? 11 ?agGGGgGGGGG! ! $!! ! ! ! ! ! !

Andamio_1?2849 A? ¿11 dólares canadienses,...,...! 0000000000

andamio_1?2850 A? 10?,...,...? 353333333

El resultado generado por compilación contiene 6 líneas: nombre de la secuencia de referencia; base de referencia; número de lectura al comparar; La quinta columna es más complicada y se explica de la siguiente manera:

1'. representa una coincidencia de cadena positiva con la secuencia de referencia.

2 ',' representa la coincidencia de cadena negativa con la secuencia de referencia.

3 ‘atcgn’ representa un desajuste en la hebra positiva.

4 "atcgn" indica una falta de coincidencia en la hebra negativa.

5' * ' representa la base difusa.

"6" ​​​​indica que la base coincidente es el comienzo de la lectura; el código ascii seguido de "" menos 33 indica la calidad de la comparación, estos dos símbolos modifican la base a continuación y la base (; ., ATCGatcgNn) Indica la primera lectura de base.

7' $ ' representa el final de una lectura, este símbolo modifica la base anterior a ella.

8 La fórmula regular "\+[0-9]+[ACGTNacgtn]+" representa la base insertada después de este sitio; por ejemplo, en el ejemplo anterior, la base insertada después de 2847 de scaffold_1 9 bases; acggtgaag. Indica que lo más probable es que esto esté en del.

La expresión regular 9 '-[0-9]+[ACGTNacgtn]+' representa una base eliminada después de este sitio

Utilice bcftools.

Bcftools es similar a samtools y se utiliza para procesar archivos vcf (formato de llamada variante) y archivos bcf (formato de llamada binaria). El primero es un archivo de texto y el segundo es su archivo binario.

Bcftools es fácil de usar. El comando principal es el comando ver, seguido de comandos como index y cat. Los comandos index y cat son similares a los de samtools. Aquí, el hablante usa el comando de vista para realizar llamadas SNP e Indel. El uso y los ejemplos de este comando son los siguientes:

$ bcftools view[-AbFGNQSucgv][-D seq dict][-l list loci][-s list sample]

? [-I brechaSNPratio][-t mu tasa][-P varThres][-P anterior]

? [-1n grupo 1][-d min frac][-U nPerm][-X umbral permanente]

? [-T trioType]in . BCF[región]

$ BCF tools view-cvNg ABC .snp_indel.vcf

El archivo de resultados generado está en formato vcf. 10 columnas, es decir: 1 nombre de la secuencia de referencia; 2 posiciones más a la izquierda de Varianti; 3 ID de variable (no establecidas de forma predeterminada, representadas por "."); 4 genes de variante de referencia, etc. (separados por '); ,' si hay múltiples alelos); 6 variantes/masa de referencia; 7 filtros aplicados; 8 información de variante, separada por punto y coma; 9 tipos de campos de genotipo, separados por dos puntos (opcional); opcional).

Por ejemplo:

scaffold_1?2847.? ¿respuesta? AACGGTGAAG? 194.? Indel; DP = 11; VDB = 0,0401; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 8, 3; GQ 1/1:235,33,0:63

scaffold_1?3908.? ¿gramo? ¿respuesta? 111.? DP = 13; VDB = 0,0085; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 5, 7; GT:PL:GQ 1/1:144,36,0:69

Andamios_1?4500.? ¿respuesta? ¿gramo? 31.5.? DP = 8; VDB = 0,0034; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 1, 3; MQ = 42FQ = -39 GT:PL:GQ 1/1:6412 p>

Andamios_1?4581.? ¿TGGGGG? TGG145.

? Indel; DP = 8; VDB = 0,0308; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 0, 8; :45

Andamio_1?4644.? ¿gramo? ¿respuesta? 195.? DP = 21; VDB = 0,0198; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 10, 10;

Andamio_1?4827.? ¿NACAAAGA NA? 4.42.? indel; DP = 1; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 1, 0; GT:PL:GQ 0/1:40,3,0:3

Andamios_1?4854.? ¿respuesta? ¿gramo? 48?. ? DP = 6; VDB = 0,0085; FQ = -36 GT:PL:GQ 1/1:80, 9, 0 :16

Andamio_1?5120.? ¿respuesta? ¿gramo? 85?. ? DP = 8; VDB = 0,0355; af 1 = 1; AC 1 = 2; DP4 = 0, 0, 5, 3; MQ = 42FQ = -51 GT:PL:GQ 1/1:11824, 0:45 p>

La columna 8 muestra la descripción informativa de la variante, que es más importante. La descripción de la etiqueta es la siguiente:

Descripción del formato de etiqueta

Estimación de máxima verosimilitud doble de la frecuencia del alelo del sitio (AF) del primer alelo ALT de AF1

DP int profundidad de lectura sin procesar (sin filtrado de calidad)

DP4 int[4] #Línea de base directa de referencia de alta calidad, línea de base inversa de referencia, línea de base inversa alternativa y alternativa

Consenso internacional de FQ Calidad. Positivo: las muestras tienen diferentes genotipos; Negativo: de lo contrario

MQ int calidad del mapeo cuadrático medio de las lecturas cubiertas

PC2 int[2] La probabilidad de fibrilación auricular de la muestra del grupo 1 de Phred. es mayor que (menor que) Grupo 2

PCHI2 valor de chi^2 ponderado posterior doble entre las muestras del Grupo 1 y el Grupo 2

Sesgo de cadena PV4, sesgo de base, sesgo de mapQ y doble [4] Valor P para el sesgo de distancia de la cola

PCHI2 de escala entera

La permutación RP int # produce un PCHI2 más pequeño

Con y sin restricciones de triple/par CLR int Relación logarítmica de Phred de probabilidad genotípica

Configuración genotípica más probable de cadena UGT sin restricción triple

Configuración de cadena CGT más probable con restricción trío

Después de usar bcftools para obtener el resultado de la llamada variante. Los resultados deben filtrarse nuevamente. Basado principalmente en la información de la octava columna de los resultados de la comparación. La línea DP4 es particularmente importante y proporciona 4 datos: 1. El resultado de la comparación es consistente con la cadena positiva, 2. El resultado de la comparación es consistente con la cadena negativa, 3. El resultado de la comparación está en la variante de la cadena positiva, y 4 El resultado de la comparación es superior a la variante de cadena negativa. Puede establecer (valor 3 + valor 4) para que sea mayor que un cierto umbral, que se trata como una variable. Por ejemplo:

$ perl -ne 'imprime $_ si /DP4=(\d+), (\d+), (\d+)/& amp;& amp($3+$4)>= 10 & & amp($3+$4)/($1+$2+$3+$4)>= 0.8 ' snp _ indel.vcf & gtsnp_indel.final.vcf