Mapeamento de reads - lmigueel/Bioinformatica GitHub Wiki
1. Sobre
O sequenciamento de DNA é obviamente primordial entre essas técnicas, mas o sequenciamento tem uma grande limitação: mesmo com a tecnologia mais sofisticada, raramente é possível obter uma sequência de mais de cerca de 750 bp em um único experimento. Isso significa que a sequência de uma longa molécula de DNA deve ser construída a partir de uma série de sequências mais curtas. Isso é feito quebrando a molécula em fragmentos, determinando a sequência de cada um e usando um computador para pesquisar sobreposições e construir a sequência consenso. Este método shotgun é a abordagem padrão para sequenciar pequenos genomas procarióticos, mas é muito mais difícil com genomas maiores porque a análise de dados necessária torna-se desproporcionalmente mais complexa à medida que o número de fragmentos aumenta (para n fragmentos, o número de sobreposições possíveis é dado por 2n^2 - 2n). Um segundo problema com o método shotgun é que ele pode levar a erros quando regiões repetitivas de um genoma são analisadas. Quando uma sequência repetitiva é quebrada em fragmentos, muitas das peças resultantes contêm os mesmos motivos de sequência ou muito semelhantes. Seria muito fácil remontar essas sequências de modo que uma parte de uma região repetitiva fosse deixada de fora, ou mesmo conectar duas peças bem separadas do mesmo cromossomo ou de cromossomos diferentes.
Para verificar as posições dos fragmentos (reads) sequenciados em relação a um genoma de interesse executamos um MAPEAMENTO. Por exemplo, abaixo está um mapeamento executado sobre um genoma. Note que os reads das duas strains abaixo (em laranja e azul) mapeiam em lugares diferentes. No final, caso esteja executando um projeto de RNA-Seq, pode-se extrair os counts (quantificação) de cada gene. Mas isso é outro assunto, mas deixo a dica de que você pode estar mapeando errado seus reads.
Acesse aqui a lista com mais de 50 programas de alinhamento: AQUI.
2. Instalação
Para mapear reads iremos utilizar dois programas: BWA e BOWTIE2. Para isso, instale-os via conda:
conda install -c bioconda bwa
conda install -c bioconda bowtie2
Acesse também aos manuais:
3. Exemplos
Vamos analisar dados de RNA-Seq em E. Coli. Se você não viu ainda o manual de como baixar dados do SRA, acesse AQUI. Acesse os dados aqui.
Download dos dados e conversão em paired-end:
fastq-dump --split-files -A SRR030257
Esse comando irá gerar os arquivos SRR030257_1.fastq e SRR030257_2.fastq.
Download do genoma:
Site: https://www.ncbi.nlm.nih.gov/genome/?term=txid413997[Organism:exp]
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
gunzip *.gz
BWA
Primeiro criamos o index pelo BWA
bwa index GCF_000005845.2_ASM584v2_genomic.fna
Nesse momento ele irá criar quatro arquivos com extensões .amb, .ann, .bwt, .pac e .sa. Agora podemos rodar os mapeamentos com 10 processadores (-t 10), e gerar assim o arquivo .sam. Caso você esteja lidando com arquivos single-end, você deve apenas passar ele como argumento.
bwa mem -t 10 GCF_000005845.2_ASM584v2_genomic.fna SRR030257_1.fastq SRR030257_2.fastq > SRR030257.paired.sam
Vamos utilizar o samtools para visualizar as informações do alinhamento. Instalação:
conda install -c bioconda samtools
Vamos converter para .bam com 4 processadores:
samtools view -@ 4 -b SRR030257.paired.sam > SRR030257.paired.bam
Agora nós ordenamos os arquivos por posição de alinhamento para que sejam fáceis de juntar (merge) mais tarde:
samtools sort -@ 4 -o SRR030257.paired.sorted.bam SRR030257.paired.bam
Este arquivo pode ser utilizado em diversas ferramentas. Mas vamos continuar com algumas análises interessantes:
Vamos criar o index do nosso .bam para que possamos buscar informações sobre ele:
samtools index SRR030257.paired.sorted.bam
Agora vamos usar o samtools flagstat para observar um sumário do alinhamento:
samtools flagstat SRR030257.paired.sorted.bam
Caso deseje extrair apenas uma região (cromosso: star-end) de interesse, utilize o samtools view. Você pode abrir este alinhamento no IGV (Acesse AQUI.
samtools view -b SRR030257.paired.sorted.bam NC_000913.3:5400-5500 > regiao_geneX.bam
BOWTIE2
Após gerar o arquivo .sam, tanto pelo Bowtie2 ou pelo BWA, as análises posteriores são as mesmas, logo, apenas irei mostrar a funcionalidade do Bowtie2 e seus comandos.
Vamos primeiro criar o index. Irei criar um output chamado Ecoli_index. Serão criados arquivos .bt2 :
bowtie2-build GCF_000005845.2_ASM584v2_genomic.fna Ecoli_index
Agora podemos mapear:
bowtie2 -t -x Ecoli_index -1 SRR030257_1.fastq -2 SRR030257_2.fastq -S SRR030257.sam
As próximas etapas de conversão para .bam, ordenaçao e visualização são as mesmas.
4. Extraindo números
Contar os reads
Para saber o total de reads presentes no arquivo .bam, caso não tenha feito o flagstat, você pode rodar:
samtools view -c SRR030257.paired.bam
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 SRR030257.paired.bam > SRR030257_only_mapped.sam
samtools view -h -F 4 SRR030257.paired.bam > SRR030257_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.
- R1 não mapeado e R2 mapeado
samtools view -u -f 4 -F 264 SRR030257.paired.bam > SRR030257_unmap_map.bam
- R1 mapeado e R2 não mapeado
samtools view -u -f 8 -F 260 SRR030257.paired.bam > SRR030257_map_unmap.bam
- R1 e R2 mapeados
samtools view -u -f 12 -F 256 SRR030257.paired.bam > SRR030257_map_map.bam
Lembrando que para cada tipo de análise deve ser utilizado corretamente os argumentos de cada mapeador. Por exemplo, análises de EXOMA possuem algumas restrições de match ou mismatches, que são diferentes de análises CNV/SNP.
5. Citação
Brown TA. Genomes. 2nd edition. Oxford: Wiley-Liss; 2002. Available from: https://www.ncbi.nlm.nih.gov/books/NBK21128/
https://www.rna-seqblog.com/youve-been-aligning-your-rna-seq-reads-all-wrong/