O mapeamento tem diferentes objetivos.
+Aqui seguiremos por dois caminhos:
+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:
+conda 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 | +
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 /databases/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 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:/databases/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 /databases/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa
+#Mapeia os reads ao indice
+bwa mem /databases/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa Paired_R1.fastq.gz Paired_R2.fastq.gz -o 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
+# -o = nome do arquivo de saída
+
+samtools mpileup -d 50000 --reference /databases/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 /databases/big/ref/Sars_cov_2.ASM985889v3.dna.toplevel.fa -g /databases/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.
+conda activate pangolin
+pangolin variant.fa
+conda deactivate
+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.
+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
+conda 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 /databases/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 /databases/big/input/data/yeast/reads/list.txt`; do trimmomatic PE -threads 4 /databases/big/input/data/yeast/reads/$i\_1.fq.gz /databases/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:/databases/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.
+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 /databases/big/ref/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa --sjdbGTFfile /databases/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
+Após cria o índice, iremos mapear as sequências limpas ao genoma de referência
+for i in `cat /databases/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.
+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 /databases/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
+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 /databases/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 /databases/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
+