Aula Mapeamento

O mapeamento tem diferentes objetivos.

Aqui seguiremos por dois caminhos:

Identificação de variante em SARS-CoV-2

Para identificar as variantes, iremos mapear as sequências ao genoma referência, identificar os pontos de mutação e caracterizar a variante.

Primeiro vamos ativar o ambiente que tem os programas instalados:

micromamba activate mapping

Escolha sua amostra

Amostras Amostras Amostras Amostras Amostras
IRR1820_S1 IRR1821_S2 IRR1822_S3 IRR1823_S4 IRR1824_S5
IRR1825_S6 IRR1826_S7 IRR1827_S8 IRR1828_S9 IRR1829_S10
IRR1830_S11 IRR1831_S12 IRR1832_S13 IRR1833_S14 IRR1834_S15
IRR1835_S16 IRR1836_S17 IRR1837_S18 IRR1838_S19 IRR1839_S20
IRR1840_S21 IRR1841_S22 IRR1842_S23 IRR1843_S24 IRR1844_S25
IRR1845_S26 IRR1846_S27 IRR1847_S28 IRR1848_S29 IRR1849_S30

Controle de qualidade

Trimmomatic vai retirar as regiões de baixa qualidade, remover os adaptadores e filtrar sequências curtas. Dando um ls repare como usando $AMOSTRA vc completa o nome do arquivo R1 e R2 com a variável exportada que tem o código da sua.

trimmomatic PE -threads 4 $AMOSTRA\_L001_R1_001.fastq.gz $AMOSTRA\_L001_R2_001.fastq.gz Paired_R1.fastq.gz Unpaired_R1.fastq.gz Paired_R2.fastq.gz Unpaired_R2.fastq.gz ILLUMINACLIP:/home/bioufmg/big/ref/adapters.fasta:2:30:10:2:keepBothReads LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:50

Mapeamento ao genoma referência

Antes de mapear, precisamos criar um índice que o programas BWA consiga utilizar. Usamos o genoma de referência para criar o índice.

Em seguida, mapeamos os reads ao índice, criando um alinhamento no formato binário "bam" ("sam" daria pra ler, "bam" só computador lê). Vamos usar somente as saídas "Paired" do Trimmomatic, os reads R1 e R2 que produziram forward e reverse.

#Cria o índice usando o arquivo FASTA de referência
bwa index /home/bioufmg/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa
#Mapeia os reads ao indice
bwa mem /home/bioufmg/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa Paired_R1.fastq.gz Paired_R2.fastq.gz > map.bam

Ordenar o mapeamento e empilhar as sequências

Para facilitar a identificação de variantes, é importante ordenarmos os reads e empilharmos para olharmos a informação que cobre cada base. Para isso, usaremos a ferramenta samtools, que trabalha com um arquivo BAM.

A função sort serve para ordenar. E a função mpileup serve empilhar as sequências em relação ao genoma de referência.

samtools sort map.bam -o map.sorted.bam
# -o = nome do arquivo de saída

samtools mpileup -d 50000 --reference /home/bioufmg/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa -a -Q 30 map.sorted.bam > pile.bam
# -d = número máximo de sequências empilhadas em uma mesma posição.
# --reference = referência utilizada no mapeamento
# -a = reportar até as regiões que não tiveram nenhuma sequência cobrindo.
# -Q = qualidade mínima do mapeamento
# > = "Salvar como" pile.bam

Identificação das variantes

Compara as sequências mapeadas com a referência, e aplica teste estatístico para eliminar a hipótese de acaso. Para realizar esse teste, usamos a ferramenta ivar, embora existam diversas outras que fazer abordagens semelhantes.

cat pile.bam | ivar variants -p variant -q 30 -t 0.1 -r /home/bioufmg/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa -g /home/bioufmg/big/ref/Sars_cov_2.ASM985889v3.101.gff3

# cat pile.bam | = usar o conteúdo do arquivo como entrada para o próximo comando
# variants = função para chamada de SNV
# -p = prefixo do nome do seu arquivo de saída
# -q = qualidade mínima do mapeamento para ser considerada na chamada de variante
# -t = proporção mínima de uma variante para ser chamada de "SNV"
# -r = FASTA do genoma de referência
# -g = GFF (arquivo de anotações) para predição do impacto na tradução dos genes

cat pile.bam | ivar consensus -p variant -q 30 -t 0.5 -m 30 -n N

# consensus = função para construção de um FASTA
# -p = prefixo do nome do arquivo de saída
# -q = qualidade mínima do mapeamento para o cálculo do consenso
# -t = proporção mínima de uma variante para ser incluída na sequência consenso.
# -m = quantidade mínima de reads na posição para definir qual base será escrita
# -n = caracter da base a ser incluída quandos os critérios mínimos não forem satisfeitos

Identificar a variante de SARS-CoV-2

Uma das utilizadades da identificação de SNV é a caracterização de variantes de virus. No exemplo, usaremos o pangolin que é uma ferramenta que agrupa diversos modelos para identificação das variantes de SARS-CoV-2 a partir da sequência de nucleotídeos deles.

micromamba activate pangolin 
pangolin variant.fa
micromamba deactivate

Mapeamento para quantificação de RNA (Genoma referência)

Outra importância do mapeamento de sequências é para fins de quantificação.

Neste exemplo iremos quantificar RNA sequenciados ao mapeá-los a um genoma de referência.

Controle de qualidade

O controle de qualidade de reads a serem mapeados para quantificação não precisa ser tão restritivo quanto o para identificação de variantes, uma vez que o mapeamento serve apenas para contar o número de evidências sobre a expressão de um gene. Para isso, usaremos o trimmomatic mas um parâmetro de qualidade mais permissivo.

Vamos criar o ambiente para trabalharmos, com estrutura de diretórios.

#Ative o ambiente que contém os programas
micromamba activate mapping

#Crie e entre no diretório Map
mkdir Map
cd Map

#Crie diretórios para escrever os resultados
mkdir trim logs map

#Veja os arquivos de entrada
ls /home/bioufmg/big/input/data/yeast/reads

Neste caso, temos uma lista com 9 arquivos, e vamos criar um comando para automatizar as análises.

for i in `cat /home/bioufmg/big/input/data/yeast/reads/list.txt`; do trimmomatic PE -threads 4 /home/bioufmg/big/input/data/yeast/reads/$i\_1.fq.gz /home/bioufmg/big/input/data/yeast/reads/$i\_2.fq.gz trim/$i.Paired_R1.fastq.gz trim/$i.Unpaired_R1.fastq.gz trim/$i.Paired_R2.fastq.gz trim/$i.Unpaired_R2.fastq.gz ILLUMINACLIP:/home/bioufmg/big/ref/adapters.fasta:2:30:10:2:keepBothReads LEADING:5 TRAILING:5 SLIDINGWINDOW:4:15 MINLEN:50 2> logs/$i.trim.log; done

#Mudamos o critério de qualidade média de 20 para 15.

Criação do índice para mapeamento

Para este mapeamento, não usaremos o BWA. O BWA é uma ferramenta que faz mapeamento contínuo e muito preciso. Precisamos de um mapeamento "splice aware", capaz de separar a sequência para escapar a região intrônica. Usaremos a ferramenta STAR que é rápido e é capaz de mapear sequências que contém introns.

A primeira etapa é criar um índice que o STAR seja capaz de ler.

STAR --runMode genomeGenerate --runThreadN 4 --genomeDir map/STAR_ref --genomeFastaFiles /home/bioufmg/big/ref/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa --sjdbGTFfile /home/bioufmg/big/ref/Saccharomyces_cerevisiae.R64-1-1.112.gtf

# --runMode = programa a ser usado. Neste caso é o genomeGenerate para criar o índice
# --runThreadN = número de processadores
# --genomeDir = Nome da pasta onde irá escrever o índice
# --genomeFastaFiles = Caminho para o arquivo FASTA que contém o genoma referência
# --sjdbGTFfile = Caminho para o arquivo de anotação, que contém as posições de cada gene

Mapeamento

Após cria o índice, iremos mapear as sequências limpas ao genoma de referência

for i in `cat /home/bioufmg/big/input/data/yeast/reads/list.txt`; do STAR --runMode alignReads --runThreadN 4 --genomeDir map/STAR_ref/ --readFilesIn trim/$i.Paired_R1.fastq.gz trim/$i.Paired_R2.fastq.gz --readFilesCommand zcat --outFilterMultimapNmax 1 --outSAMtype BAM SortedByCoordinate --outFileNamePrefix map/$i. --quantMode GeneCounts; done

# Para cada arquivo na nossa lista, irá rodar o STAR
# --runMode = programa a ser rodado, agora é o alignReads
# --runThreadN = Número de processadores
# --genomeDir = Diretório que contém o índice que vai ler
# --readFilesIn = Localização dos arquivos que contém as leituras
# --readFilesCommand = Caso os arquivos estejam compactados (.gz) diz qual o comando que descompacta
# --outFilterMultimapNmax = Máximo de loci que cada sequência pode ser mapeada
# --outSAMtype = Tipo do formato de saída, no caso, será o formato BAM já ordenado
# --outFileNamePrefix = Padrão dos nomes dos arquivos de saída
# --quantMode = Informa se irá quantificar os reads mapeados, no caso por contagem dos genes (GeneCounts)

Isso vai gerar arquivos de quantificação para cada amostra. A quantificação do STAR é bem restritiva, e muitos reads acabam não sendo quantificados.

Quantificação

Vamos então quantificar usando uma ferramenta do bedtools que serve para quantificar utilizando como referência um arquivo de anotações. Filtrei somente as anotações de CDS. E então rodamos o comando:

# Primeira etapa é gerar índices para os arquivos BAM
for i in `ls *.bam`; do samtools index $i; done
# Para cada arquivo BAM que temos, rode o indexador do samtools

# Então quantificamos
multiBamCov -bams etoh60_1.Aligned.sortedByCoord.out.bam etoh60_2.Aligned.sortedByCoord.out.bam etoh60_3.Aligned.sortedByCoord.out.bam temp33_1.Aligned.sortedByCoord.out.bam temp33_2.Aligned.sortedByCoord.out.bam temp33_3.Aligned.sortedByCoord.out.bam -bed /home/bioufmg/big/ref/Saccharomyces_cerevisiae.R64-1-1.cds.gff > Quantification.tsv

# -bams = lista de arquivos BAM
# -bed = arquivo que tem a posição de cada CDS no genoma

Quantificação sem mapeamento

Outra forma de mapear é utilizando ferramentas que fazem a contagem dos trancritos sem mapear. Para isso não usamos o genoma referência, mas sim um arquivo contendo todos os transcritos (CDS) do organismo.

Neste exemplo, usaremos o salmon, que é uma ferramenta que quantifica e normaliza usando sequências de kmer para inferir o transcrito alvo.

A primeira etapa (como sempre) é criar um índice com a informação das CDS.

salmon index -t /home/bioufmg/big/ref/Saccharomyces_cerevisiae.R64-1-1.cdna.all.fa.gz -i salmon/index/Sce -p 4

# index = indica o programa a ser executado
# -t = arquivo que contém os transcritos em formato FASTA
# -i = nome da pasta que vai guardar os índices
# -p = número de processadores

Após cria o índice, podemos quantificar nossos reads limpos utilizando o transcritoma referência.

for i in `cat /home/bioufmg/big/input/data/yeast/reads/list.txt`; do salmon quant -i salmon/index/Sce/ -l A -1 trim/$i.Paired_R1.fastq.gz -2 trim/$i.Paired_R2.fastq.gz -p 4 -o salmon/$i; done

# Para cada amostra da nossa lista, vai rodar o salmon
# quant = programa do salmon a ser rodado
# -i = pasta que contém o index
# -l = Orientação da biblioteca, no caso usamos o A, que refere-se a Automático
# -1 = Arquivo com as sequências diretas
# -2 = Arquivo com as sequências reversas
# -p = Número de processadores
# -o = nome da pasta em que escreverá os arquivos de saída

O salmon vai gerar uma pasta para cada amostra, e cada uma tem uma quantificação. Vamos juntar todas as quantificações em uma única tabela para as futuras análises estatísticas.

salmon quantmerge --quants salmon/etoh60_1/ salmon/etoh60_2/ salmon/etoh60_3/ salmon/ref1/ salmon/ref2/ salmon/ref3/ salmon/temp33_1/ salmon/temp33_2/ salmon/temp33_3/ --names etoh1 etoh2 etoh3 ref1 ref2 ref3 temp1 temp2 temp3 -o salmon/quant.txt

# --quants = Pastas que tem as quantificações
# --names = Nomes que as colunas da sua tabela terão (na mesma ordem que os nomes das pastas)
# -o = Nome do arquivo de saída