Manipulação com samtools - lmigueel/Bioinformatica GitHub Wiki

1. Sobre

O formato SAM (Sequence Alignment/Map) é um formato genérico para armazenar grandes alinhamentos de sequências de nucleotídeos. O SAM pretende ser um formato que:

 1. É flexível o suficiente para armazenar todas as informações de alinhamento geradas por vários programas de alinhamento;
 2. É simples o suficiente para ser facilmente gerado por programas de alinhamento ou convertido de formatos de alinhamento existentes;
 3. É compacto em tamanho de arquivo;
 4. Permite que a maioria das operações no alinhamento funcionem em um fluxo sem carregar todo o alinhamento na memória;
 5. Permite que o arquivo seja indexado por posição genômica para recuperar com eficiência todas as leituras alinhadas a um locus.

O samtools fornece vários utilitários para manipular alinhamentos no formato SAM, incluindo classificação, mesclagem, indexação e geração de alinhamentos em um formato por posição, arquivos BAM, entre outros.

2. Instalação

A instalação pode ser realizada pelo conda:

 conda install -c bioconda samtools 

3. Exemplos de uso

1. Verificar o alinhamento de BAM

samtools view sample.bam | head

2. Ordenar arquivo BAM

samtools sort sample.bam -o sample.sorted.bam

3. Extrair genes mapeados de uma região

Primeiro criamos o index do BAM ordenado:

samtools index sample.sorted.bam

E extaimos os reads de uma região, por exemplo, cromossomo 1 (chr1) da posição 1000 até 3000. Cuidado com o header.

samtools view sample.sorted.bam chr1:1000-3000 > geneX.sorted.bam

Caso queiramos contar quantos reads estão naquela região:

samtools view sample.sorted.bam chr1:1000-3000 | wc -l

4. SAM para BAM

samtools view -bS teste.sam > teste.bam

5. SAM para BAM já ordenado

samtools view -bS file.sam | samtools sort - file_sorted

6. BAM para SAM

samtools view -h file.bam > file.sam

7. Sumário do alinhamento

samtools flagstat file.sorted.bam

7. BAM para FASTQ

Caso seja paired-end:

samtools fastq -1 file_1.fastq -2 file_2.fastq file.bam

Caso seja single-end:

samtools fastq file.bam > file.fastq

8. Reads mapeados e não mapeados

Para você extrair os reads totais que mapearam utilizamos a flag -F 4, e para os não mapeados utilizamos a flag -f 4 :

samtools view -h -F 4 file.paired.bam > file_only_mapped.sam

samtools view -h -F 4 file.paired.bam > file_unmapped.sam

Caso você deseje saber os reads mapeados e não mapeados, divididos por arquivo de entrada (R1 e R2), você pode executar os comandos abaixo.

1. R1 não mapeado e R2 mapeado

samtools view -u -f 4 -F 264 file.paired.bam > file_unmap_map.bam

2. R1 mapeado e R2 não mapeado

samtools view -u -f 8 -F 260 file.paired.bam > file_map_unmap.bam

3. R1 e R2 mapeados

samtools view -u -f 12 -F 256 file.paired.bam > file_map_map.bam

9. Ver alinhamento no terminal

Você pode observar os alinhamentos pelo terminal, e não por um software, através do tview, passando como argumento o BAM e o FASTA do mapeamento. Quando abrir a janela, você precisa clicar no 'g' do teclado e digitar as coordenadas, por exemplo: chr1:1000-3000.

samtools tview file.bam file.fasta

Citação

Li H., Handsaker B., Wysoker A., Fennell T., Ruan J., Homer N., Marth G., Abecasis G., Durbin R. and 1000 Genome Project Data Processing Subgroup (2009) The Sequence alignment/map (SAM) format and SAMtools. Bioinformatics, 25, 2078-9. [PMID: 19505943]