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
| 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 |
Entre na sua pasta
cd eusoujacu
Crie sua pasta e copie os arquivos
mkdir sars
cd sars
#Variável guarda o nome da sua amostra
export AMOSTRA="AMOSTRA_ESCOLHIDA" #Substitua o conteúdo entre aspas pelo código da amostra da tabela que escolheu
cp ~/big/input/map/$AMOSTRA*.fastq.gz .
Trimmomatic vai retirar as regiões de baixa qualidade, remover os adaptadores e filtrar sequências curtas. Dando um ls repare como usando $AMOSTRA você 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:~/big/ref/adapters.fasta:2:30:10:2:keepBothReads LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:50
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 ~/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa
#Mapeia os reads ao indice
bwa mem ~/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa Paired_R1.fastq.gz Paired_R2.fastq.gz > map.bam
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
# > = salva como arquivo de saída
samtools mpileup -d 50000 --reference ~/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
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 ~/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa -g ~/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
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.
pangolin variant.fa
## Para verificar de forma mais limpa os resultados das variantes
cat lineage_report.csv | awk -F "," '{print $2"\t"$4"\t"$5}'