diff --git a/elab/index.html b/elab/index.html index cdf9db7..93fe1af 100755 --- a/elab/index.html +++ b/elab/index.html @@ -1,3192 +1,3192 @@ - - - - - Conteúdo ELAB - - - - - - + + + + + - - - - - - -
-

Conteúdo ELAB

-

Aula 1 - Aleatoreidade

-

Aula 2 - Anotação

-

Tutorial - Montagem e anotação

-

Download do material da aula no - Google Drive

-
-
- - - - +} + + + + + + + +
+

Conteúdo ELAB

+

Aula 1 - Aleatoreidade

+

Aula 2 - Anotação

+

Tutorial - Montagem e anotação

+

Tutorial - Mapeamento e quantificação

+

Download do material da aula no - Google Drive

+
+
+ + + + diff --git a/elab/tutorial/index.html b/elab/tutorial/index.html index ca16b52..162e4b2 100755 --- a/elab/tutorial/index.html +++ b/elab/tutorial/index.html @@ -1,3193 +1,3193 @@ - - - - - Montagem e anotação funcional - - - - - - - + + + + + + - - - - - - -
-
-
-

Montagem e anotação funcional

-
-

Objetivo: Montagem e anotação funcional de um genoma bacteriano.

-
-
-

Origem dos dados

-

Os dados utilizados são provenientes de um sequenciamento feito na plataforma Illumina. Uma biblioteca de DNA foi utilizada para duas diferentes abordagens de sequenciamento, porém com a mesma profundidade esperada de 60x.

- -

Conhecendo o conteúdo da pasta

- -
pwd #path of working directory
-
-ls #lista o conteudo da pasta
-
- -
mkdir Aula_Funcional
-cd Aula_Funcional
-ls
-
-

Baixe o arquivo com o conteúdo da aula e descompacte dentro da pasta.

-
tar -zxf elab.tar.gz
-
- -
mkdir 0.Logs 1.Trim 2.Assembly 3.Genes 4.Function
-ls
-
-conda activate assembly
-
-
-

Controle de qualidade das sequências

-

O sequenciamento Illumina tem uma taxa de erro de substituição de bases de cerca de 1%. À medida em que a sequência vai sendo sintetizada essa taxa de erro pode aumentar. Por isso, é importante retirarmos as regiões com baixa qualidade de bases.

-

A ferramenta Trimmomatic é uma das muitas utilizadas para esse fim. Nela, podemos definir a qualidade média das bases em uma janela deslizante, além da remoção de adaptadores e exclusão de sequências muito curtas. -O comando utilizado é:

-
trimmomatic PE -threads 4 data/lib1_150bp_R1.fq.gz data/lib1_150bp_R2.fq.gz 1.Trim/lib1.Paired_R1.fastq.gz 1.Trim/lib1.Unpaired_R1.fastq.gz 1.Trim/lib1.Paired_R2.fastq.gz 1.Trim/lib1.Unpaired_R2.fastq.gz ILLUMINACLIP:data/adapters.fasta:2:30:10:2:keepBothReads LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:50 2> 0.Logs/lib1.trimmomatic.log
-
-

A linha de comando recebe parâmetros que chamamos de flags. Seccionando as linhas de comando, vemos a função de cada parâmetro.

-
trimmomatic #inicia o programa
-PE #informa que o sequenciamento é paired-end
--threads 4 #conta que vai usar 4 processadores
-data/lib1_150bp_R1.fq.gz data/lib1_150bp_R2.fq.gz #indicas os arquivos que contem as sequências forward e reversas a serem corrigidas
-1.Trim/lib1.Paired_R1.fastq.gz #Indica onde vai escrever as sequências Forward (R1) cujos pares sobreviveram
-1.Trim/lib1.Unpaired_R1.fastq.gz #Sequências forward cujos pares reversos não passaram no critério de qualidade
-1.Trim/lib1.Paired_R2.fastq.gz 1.Trim/lib1.Unpaired_R2.fastq.gz #Os mesmos que os anteriores, só que com as sequências reversas: R2
-ILLUMINACLIP:data/adapters.fasta:2:30:10:2:keepBothReads #Cita o nome do arquivo que contém a sequência dos adaptadores : número de mismatches aceitos : tamanho mínimo de um palíndromo (em caso de pares de reads) : tamanho mínimo do alinhamento simples : tamanho mínimo da sequência do adaptador identifica.
-LEADING:5 TRAILING:5 #Valores mínimos de qualidade para começar a análise no início (LEADING) e no final da sequência (TRAILING)
-SLIDINGWINDOW:4:20 #Tamanho da janela de bases : Qualidade média dentro da janela
-MINLEN:50 #Tamanho mínimo da sequência restante para que não seja descartada
-2> 0.Logs/lib1.trimmomatic.log #Escreve o relatório da análise no arquivo mencionado
-
-

Assim, dentro da pasta 1.Trim teremos os arquivos limpos. E no 0.Logs podemos ver o relatório.

-
cat 0.Logs/lib1.trimmomatic.log
-
-

Montagem do genoma

-

Para a montagem do genoma, utilizaremos apenas as sequências cujos pares "sobreviveram" ao controle de qualidade. Assim, teremos dois arquivos de entrada: as sequências forward e as reverse. -Dentre os vários programas para a montagem de um genoma, vamos utilizar o Spades.

-
spades.py --isolate -o 2.Assembly/lib1 -1 1.Trim/lib1.Paired_R1.fastq.gz -2 1.Trim/lib1.Paired_R2.fastq.gz -t 4 -m 16
-
-

Esse comando contém os seguintes parâmetros:

-
spades.py #Roda o programa
---isolate #Diz que o genoma é de um isolado. Ou seja, espera-se que todas as sequências tenham vindo de um mesmo organismo
--o 2.Assembly/lib1 #Diretório onde escreverá os resultados
--1 1.Trim/lib1.Paired_R1.fastq.gz #Arquivo forward
--2 1.Trim/lib1.Paired_R2.fastq.gz #Arquivo reverse
--t 4 #Número de processadores
--m 16 #Máximo de memória
-
-
-

Avaliando a qualidade da montagem

-

Existem métrica que informam sobre o quão completa está a montagem do seu genoma. Uma das ferramentas é o Quast, que fornece métricas como:

- -

Para rodar, podemos usar o comando com os seguintes parâmtros:

-
quast -o 2.Assembly/lib1/quast/ -m 1000 -t 4 -1 1.Trim/lib1.Paired_R1.fastq.gz -2 1.Trim/lib1.Paired_R2.fastq.gz 2.Assembly/lib1/scaffolds.fasta
-
-

Os parâmetros são:

-
quast #Chama o programa
--o 2.Assembly/lib1/quast/ #Diz onde vão escrever os resultados
--m 1000 #Tamano mínimo de um contig
--t 4 #Número de processadores
--1 1.Trim/lib1.Paired_R1.fastq.gz #[OPCIONAL] Sequências forward utilizadas na montagem
--2 1.Trim/lib1.Paired_R2.fastq.gz #[OPCIONAL] Sequências reverse utilizadas na montagem
-2.Assembly/lib1/scaffolds.fasta #Arquivo com sua montagem
-
-

Observando o conteúdo do arquivo, vemos as métricas da montagem

-
cat 2.Assembly/lib1/quast/report.txt
-
-
Assembly                    scaffolds
-# contigs (>= 0 bp)         212      
-# contigs (>= 1000 bp)      81       
-# contigs (>= 5000 bp)      62       
-# contigs (>= 10000 bp)     56       
-# contigs (>= 25000 bp)     46       
-# contigs (>= 50000 bp)     29       
-Total length (>= 0 bp)      4502140
-Total length (>= 1000 bp)   4473025
-Total length (>= 5000 bp)   4431666
-Total length (>= 10000 bp)  4386681
-Total length (>= 25000 bp)  4227455
-Total length (>= 50000 bp)  3592230
-# contigs                   81       
-Largest contig              329901
-Total length                4473025
-GC (%)                      50.79
-N50                         115885
-N90                         34661
-auN                         148131.5
-L50                         12
-L90                         40
-# total reads               1998966  
-# left                      999483   
-# right                     999483   
-Mapped (%)                  99.28
-Properly paired (%)         98.67
-Avg. coverage depth         66
-Coverage >= 1x (%)          100.0
-# N's per 100 kbp           4.47 
-
-

Predição de Genes

-

A predição de genes vai procurar por janelas abertas de leitura (ORF) e aprender sobre o uso de codons dentro dessa janela para predizer os genes possivelmente codificantes.

-

Dentre as diversas ferramentas, usaremos a ferramenta prodigal

-
prodigal -a 3.Genes/lib1.aa -d 3.Genes/lib1.nt -f gff -i 2.Assembly/lib1/scaffolds.fasta -o 3.Genes/lib1.gff
-
-

O comando contém as seguintes flags

-
prodigal #Chama o programa
--a 3.Genes/lib1.aa #Escreve um arquivo contendo as sequências de aminoácidos
--d 3.Genes/lib1.nt #Escreve um arquivo contendo os nucleotídeos
--f gff #Determina o formato GFF para a anotação
--i 2.Assembly/lib1/scaffolds.fasta #Conta qual o arquivo contém o genoma para a predição
--o 3.Genes/lib1.gff #Escreve um arquivo contendo a anotação
-
-

Podemos observar as primeiras linha do arquivo:

-
head 3.Genes/lib1.aa
-
-

No cabeçalho temos algumas informações: ->NODE_1_length_329901_cov_29.783351_1 # 19 # 486 # 1 # ID=1_1;partial=00;start_type=ATG;rbs_motif=AGGA;rbs_spacer=5-10bp;gc_cont=0.538.

-

A informação partial=00 informa se o gene foi predito completamente - Metionina no início e codon de parada ao final. partial=10 indica que não havia Metionina no início, mas havia codon de parada. partial=01 indica que o codon de parada não estava presente. E partial=11 indica que não havia metionina e nem codon de parada, porém o preditor encontrou uma sequência de possível codons que condiziam com codon usage do organismo.

-

Podemos contar quantos genes codificantes de proteínas foram preditos:

-
grep ">" -c 3.Genes/lib1.aa
-
-

O que esse comando faz?

-
grep #Chama o programa que faz buscas
-">" #Aponta o que ele precisa buscar
--c #Retorna a contagem de linhas que tinham aquele item da busca ao invés de retornar a linha inteira
-3.Genes/lib1.aa #Arquivo em que a busca será feita
-
-

Podemos usar o grep para contar quantos genes estavam completos:

-
grep ">" 3.Genes/lib1.aa | cut -d "#" -f5 | cut -d ";" -f2 | sort | uniq -c
-
-

Explicando...

-
grep #Chama o programa que faz as buscas
-">" #Aponta o que ele precisa buscar
-3.Genes/lib1.aa #Arquivo em que será feita a busca
-| #Conecta o resultado do comando anterior ao próximo comando
-cut #Corta o arquivo em colunas
--d "#" #Determina que o caracter delimitador de colunas é o #
--f5 #Mostra a coluna 5
-| #Conecta o resultado do comando anterior ao próximo comando
-cut -d ";" -f2 | #Corta usando o ; como delimitador e mostra a segunda coluna
-sort | #Ordena em ordem alfabética
-uniq -c #Conta quanta vezes cada item único apareceu
-
-

Classificação funcional dos genes

-

Como vimos, temos milhares de possíveis genes, e não sabemos suas funções. -Uma das formas de identificar as funções é através de alinhamentos locais com bancos de dados secundários que contenham alguma informação funcional. -Para tal tarefa, o BLAST é a ferramenta mais popular. Usaremos o diamond que faz um alinhamento local, assim como o BLAST, porém de maneira muito mais rápida.

-

Utilizaremos o UniProt como banco de dados para as buscas. A primeira etapa é transformar as sequências do UniProt em um formato que o diamond possa ler.

-
diamond makedb --threads 4 --in data/uniprotkb_proteome_UP000000625_2024_08_29.fasta.gz --db ecoli
-
-

As flags significam:

-
diamond #Chama o programa
-makedb #Conta qual a tarefa que ele vai fazer
---threads 4 #Número de processadores
---in data/uniprotkb_proteome_UP000000625_2024_08_29.fasta.gz #Onde está o arquivo do UniProt com as sequências
---db ecoli #Nome que quer dar ao seu banco de dados
-
-

Após formatarmos o banco de dados, podemos chamar o alinhamento. Temos um banco de dados de sequências de aminoácidos (UniProt), e arquivo de entrada de aminoácidos ou nucleotídeos. Poderiam então escolher entre duas opções:

- -

Para otimizar o tempo - dado que temos as sequências de aminoácidos - rodaremos o blastp:

-
diamond blastp -d ecoli.dmnd -q 3.Genes/lib1.aa -o 4.Function/lib1.diamond -e 1e-20 -k 1 -p 4 --outfmt 6
-
-

As flags são:

-
diamond #Chama o programa
-blastp #Diz a opção de alinhamento
--d ecoli.dmnd #Aponta o banco de dados
--q 3.Genes/lib1.aa #Aponta o arquivo de entrada 
--o 4.Function/lib1.diamond #Define o nome do arquivo de saída
--e 1e-20 #Estabelece um mínimo de e-value
--k 1 #Pede para mostrar só o melhor resultado para cada um dos meus genes
--p 4 #Usa 4 processadores
---outfmt 6 #Mostra o resultado em um formato tabular com 12 colunas
-
-

Podemos ver como é o resultado:

-
head 4.Function/lib1.diamond
-
-#Descrição das colunas
- # 1.  qseqid      query or source (gene) sequence id
- # 2.  sseqid      subject or target (reference genome) sequence id
- # 3.  pident      percentage of identical positions
- # 4.  length      alignment length (sequence overlap)
- # 5.  mismatch    number of mismatches
- # 6.  gapopen     number of gap openings
- # 7.  qstart      start of alignment in query
- # 8.  qend        end of alignment in query
- # 9.  sstart      start of alignment in subject
- #10.  send        end of alignment in subject
- #11.  evalue      expect value
- #12.  bitscore    bit score
-
-

Vemos que a segunda coluna contém o nome da proteína que melhor descreve nosso gene no banco de dados. -Agora queremos adicionar informação funcional à ela. -Vamos adicionar a informação do KEGG, e saber qual KO e via metabólica nosso gene pertence.

-

Para isso vamos executar uma sério de comandos.

-
#Entrar na pasta que contém os resultados
-cd 4.Function
-
-#Pega o nome de cada proteína que teve hit
-awk -F'|' '{print $2}' lib1.diamond > protein_ids.txt
-head protein_ids.txt
-
-#Adicionar o nome dessa proteína do UniProte à nossa tabela do resultado do Diamond
-paste lib1.diamond protein_ids.txt > diamond_with_protein_ids.txt
-
-# Colocar em ordem alfabética as colunas que contem as informações que vamos precisar.
-## A coluna que contém o nome das proteínas é a 13.
-sort -k13 diamond_with_protein_ids.txt > sorted_diamond.txt
-## A coluna que contém os nomes das mesmas proteínas no arquivo de conversão, é a 1
-sort -k1 data/Uniprot2KEGG.eco.txt > sorted_Uniprot2KEGG.txt
-## Juntamos as duas tabelas, levando em consideração o ponto de junção como sendo a 13 coluna do resultado do Diamond, com a 1 coluna da tabela de conversão.
-join -t$'\t' -1 13 -2 1 sorted_diamond.txt sorted_Uniprot2KEGG.txt > diamond_with_KEGG.txt
-
-# Colocar em ordem alfabética as colunas que contem as informações que vamos precisar.
-## A coluna que contém o nome da entrada no KEGG é a 14.
-sort -k14 diamond_with_KEGG.txt > sorted_diamond_with_KEGG.txt
-## A coluna que contém os nomes das mesmas proteínas no arquivo de conversão, é a 1
-sort -k1 data/KEGG2KO.eco.txt > sorted_KEGG2KO.txt
-## Juntamos as duas tabelas, levando em consideração o ponto de junção como sendo a 14 coluna do resultado do Diamond, com a 1 coluna da tabela de conversão.
-join -t$'\t' -1 14 -2 1 sorted_diamond_with_KEGG.txt sorted_KEGG2KO.txt > diamond_with_KO.txt
-
-# Temos uma tabela com 15 colunas, sendo que a última delas é o identificador do KEGG Orthology
-## Podemos ver qual é o KO mais representado
-cat diamond_with_KO.txt | cut -f15 | sort | uniq -c | sort -k1 -nr | head
-## Podemos gerar uma lista para visualizar no mapa metabólico
-cat diamond_with_KO.txt | cut -f15 | sort -u | sed -e 's/ko://g' > KO_list.txt
-
-

Visualização em vias metabólicas

-

A partir da lista de KOs, é possível ver suas funções no site do KEGG ou criar uma via metabólica interativa com o iPath3.

-
-
- - - - +} + + + + + + + +
+
+
+

Montagem e anotação funcional

+
+

Objetivo: Montagem e anotação funcional de um genoma bacteriano.

+
+
+

Origem dos dados

+

Os dados utilizados são provenientes de um sequenciamento feito na plataforma Illumina. Uma biblioteca de DNA foi utilizada para duas diferentes abordagens de sequenciamento, porém com a mesma profundidade esperada de 60x.

+ +

Conhecendo o conteúdo da pasta

+ +
pwd #path of working directory
+
+ls #lista o conteudo da pasta
+
+ +
mkdir ELAB
+cd ELAB
+ls
+
+

Baixe o arquivo com o conteúdo da aula e descompacte dentro da pasta.

+
tar -zxf elab.tar.gz
+
+ +
mkdir 0.Logs 1.Trim 2.Assembly 3.Genes 4.Function
+ls
+
+####### Caso não tenha o ambiente com os programas instalados, pode instalar com: #####
+conda create -n assembly -c bioconda trimmomatic spades quast prodigal diamond
+#######################################################################################
+
+#Ativar ambiente com programas instalados
+conda activate assembly
+
+
+

Controle de qualidade das sequências

+

O sequenciamento Illumina tem uma taxa de erro de substituição de bases de cerca de 1%. À medida em que a sequência vai sendo sintetizada essa taxa de erro pode aumentar. Por isso, é importante retirarmos as regiões com baixa qualidade de bases.

+

A ferramenta Trimmomatic é uma das muitas utilizadas para esse fim. Nela, podemos definir a qualidade média das bases em uma janela deslizante, além da remoção de adaptadores e exclusão de sequências muito curtas. +O comando utilizado é:

+
trimmomatic PE -threads 4 /databases/big/input/lib1_150bp_R1.fq.gz /databases/big/input/lib1_150bp_R2.fq.gz 1.Trim/lib1.Paired_R1.fastq.gz 1.Trim/lib1.Unpaired_R1.fastq.gz 1.Trim/lib1.Paired_R2.fastq.gz 1.Trim/lib1.Unpaired_R2.fastq.gz ILLUMINACLIP:/databases/big/ref/adapters.fasta:2:30:10:2:keepBothReads LEADING:5 TRAILING:5 SLIDINGWINDOW:4:20 MINLEN:50 2> 0.Logs/lib1.trimmomatic.log
+
+

A linha de comando recebe parâmetros que chamamos de flags. Seccionando as linhas de comando, vemos a função de cada parâmetro.

+
trimmomatic #inicia o programa
+PE #informa que o sequenciamento é paired-end
+-threads 4 #conta que vai usar 4 processadores
+/databases/big/input/lib1_150bp_R1.fq.gz /databases/big/input/lib1_150bp_R2.fq.gz #indicas os arquivos que contem as sequências forward e reversas a serem corrigidas
+1.Trim/lib1.Paired_R1.fastq.gz #Indica onde vai escrever as sequências Forward (R1) cujos pares sobreviveram
+1.Trim/lib1.Unpaired_R1.fastq.gz #Sequências forward cujos pares reversos não passaram no critério de qualidade
+1.Trim/lib1.Paired_R2.fastq.gz 1.Trim/lib1.Unpaired_R2.fastq.gz #Os mesmos que os anteriores, só que com as sequências reversas: R2
+ILLUMINACLIP:/databases/big/ref/adapters.fasta:2:30:10:2:keepBothReads #Cita o nome do arquivo que contém a sequência dos adaptadores : número de mismatches aceitos : tamanho mínimo de um palíndromo (em caso de pares de reads) : tamanho mínimo do alinhamento simples : tamanho mínimo da sequência do adaptador identifica.
+LEADING:5 TRAILING:5 #Valores mínimos de qualidade para começar a análise no início (LEADING) e no final da sequência (TRAILING)
+SLIDINGWINDOW:4:20 #Tamanho da janela de bases : Qualidade média dentro da janela
+MINLEN:50 #Tamanho mínimo da sequência restante para que não seja descartada
+2> 0.Logs/lib1.trimmomatic.log #Escreve o relatório da análise no arquivo mencionado
+
+

Assim, dentro da pasta 1.Trim teremos os arquivos limpos. E no 0.Logs podemos ver o relatório.

+
cat 0.Logs/lib1.trimmomatic.log
+
+

Montagem do genoma

+

Para a montagem do genoma, utilizaremos apenas as sequências cujos pares "sobreviveram" ao controle de qualidade. Assim, teremos dois arquivos de entrada: as sequências forward e as reverse. +Dentre os vários programas para a montagem de um genoma, vamos utilizar o Spades.

+
spades.py --isolate -o 2.Assembly/lib1 -1 1.Trim/lib1.Paired_R1.fastq.gz -2 1.Trim/lib1.Paired_R2.fastq.gz -t 4 -m 16
+
+

Esse comando contém os seguintes parâmetros:

+
spades.py #Roda o programa
+--isolate #Diz que o genoma é de um isolado. Ou seja, espera-se que todas as sequências tenham vindo de um mesmo organismo
+-o 2.Assembly/lib1 #Diretório onde escreverá os resultados
+-1 1.Trim/lib1.Paired_R1.fastq.gz #Arquivo forward
+-2 1.Trim/lib1.Paired_R2.fastq.gz #Arquivo reverse
+-t 4 #Número de processadores
+-m 16 #Máximo de memória
+
+
+

Avaliando a qualidade da montagem

+

Existem métrica que informam sobre o quão completa está a montagem do seu genoma. Uma das ferramentas é o Quast, que fornece métricas como:

+ +

Para rodar, podemos usar o comando com os seguintes parâmtros:

+
quast -o 2.Assembly/lib1/quast/ -m 1000 -t 4 -1 1.Trim/lib1.Paired_R1.fastq.gz -2 1.Trim/lib1.Paired_R2.fastq.gz 2.Assembly/lib1/scaffolds.fasta
+
+

Os parâmetros são:

+
quast #Chama o programa
+-o 2.Assembly/lib1/quast/ #Diz onde vão escrever os resultados
+-m 1000 #Tamano mínimo de um contig
+-t 4 #Número de processadores
+-1 1.Trim/lib1.Paired_R1.fastq.gz #[OPCIONAL] Sequências forward utilizadas na montagem
+-2 1.Trim/lib1.Paired_R2.fastq.gz #[OPCIONAL] Sequências reverse utilizadas na montagem
+2.Assembly/lib1/scaffolds.fasta #Arquivo com sua montagem
+
+

Observando o conteúdo do arquivo, vemos as métricas da montagem

+
cat 2.Assembly/lib1/quast/report.txt
+
+
Assembly                    scaffolds
+# contigs (>= 0 bp)         212      
+# contigs (>= 1000 bp)      81       
+# contigs (>= 5000 bp)      62       
+# contigs (>= 10000 bp)     56       
+# contigs (>= 25000 bp)     46       
+# contigs (>= 50000 bp)     29       
+Total length (>= 0 bp)      4502140
+Total length (>= 1000 bp)   4473025
+Total length (>= 5000 bp)   4431666
+Total length (>= 10000 bp)  4386681
+Total length (>= 25000 bp)  4227455
+Total length (>= 50000 bp)  3592230
+# contigs                   81       
+Largest contig              329901
+Total length                4473025
+GC (%)                      50.79
+N50                         115885
+N90                         34661
+auN                         148131.5
+L50                         12
+L90                         40
+# total reads               1998966  
+# left                      999483   
+# right                     999483   
+Mapped (%)                  99.28
+Properly paired (%)         98.67
+Avg. coverage depth         66
+Coverage >= 1x (%)          100.0
+# N's per 100 kbp           4.47 
+
+

Predição de Genes

+

A predição de genes vai procurar por janelas abertas de leitura (ORF) e aprender sobre o uso de codons dentro dessa janela para predizer os genes possivelmente codificantes.

+

Dentre as diversas ferramentas, usaremos a ferramenta prodigal

+
prodigal -a 3.Genes/lib1.aa -d 3.Genes/lib1.nt -f gff -i 2.Assembly/lib1/scaffolds.fasta -o 3.Genes/lib1.gff
+
+

O comando contém as seguintes flags

+
prodigal #Chama o programa
+-a 3.Genes/lib1.aa #Escreve um arquivo contendo as sequências de aminoácidos
+-d 3.Genes/lib1.nt #Escreve um arquivo contendo os nucleotídeos
+-f gff #Determina o formato GFF para a anotação
+-i 2.Assembly/lib1/scaffolds.fasta #Conta qual o arquivo contém o genoma para a predição
+-o 3.Genes/lib1.gff #Escreve um arquivo contendo a anotação
+
+

Podemos observar as primeiras linha do arquivo:

+
head 3.Genes/lib1.aa
+
+

No cabeçalho temos algumas informações: +>NODE_1_length_329901_cov_29.783351_1 # 19 # 486 # 1 # ID=1_1;partial=00;start_type=ATG;rbs_motif=AGGA;rbs_spacer=5-10bp;gc_cont=0.538.

+

A informação partial=00 informa se o gene foi predito completamente - Metionina no início e codon de parada ao final. partial=10 indica que não havia Metionina no início, mas havia codon de parada. partial=01 indica que o codon de parada não estava presente. E partial=11 indica que não havia metionina e nem codon de parada, porém o preditor encontrou uma sequência de possível codons que condiziam com codon usage do organismo.

+

Podemos contar quantos genes codificantes de proteínas foram preditos:

+
grep ">" -c 3.Genes/lib1.aa
+
+

O que esse comando faz?

+
grep #Chama o programa que faz buscas
+">" #Aponta o que ele precisa buscar
+-c #Retorna a contagem de linhas que tinham aquele item da busca ao invés de retornar a linha inteira
+3.Genes/lib1.aa #Arquivo em que a busca será feita
+
+

Podemos usar o grep para contar quantos genes estavam completos:

+
grep ">" 3.Genes/lib1.aa | cut -d "#" -f5 | cut -d ";" -f2 | sort | uniq -c
+
+

Explicando...

+
grep #Chama o programa que faz as buscas
+">" #Aponta o que ele precisa buscar
+3.Genes/lib1.aa #Arquivo em que será feita a busca
+| #Conecta o resultado do comando anterior ao próximo comando
+cut #Corta o arquivo em colunas
+-d "#" #Determina que o caracter delimitador de colunas é o #
+-f5 #Mostra a coluna 5
+| #Conecta o resultado do comando anterior ao próximo comando
+cut -d ";" -f2 | #Corta usando o ; como delimitador e mostra a segunda coluna
+sort | #Ordena em ordem alfabética
+uniq -c #Conta quanta vezes cada item único apareceu
+
+

Classificação funcional dos genes

+

Como vimos, temos milhares de possíveis genes, e não sabemos suas funções. +Uma das formas de identificar as funções é através de alinhamentos locais com bancos de dados secundários que contenham alguma informação funcional. +Para tal tarefa, o BLAST é a ferramenta mais popular. Usaremos o diamond que faz um alinhamento local, assim como o BLAST, porém de maneira muito mais rápida.

+

Utilizaremos o UniProt como banco de dados para as buscas. A primeira etapa é transformar as sequências do UniProt em um formato que o diamond possa ler.

+
diamond makedb --threads 4 --in /databases/big/ref/uniprotkb_proteome_UP000000625_2024_08_29.fasta.gz --db ecoli
+
+

As flags significam:

+
diamond #Chama o programa
+makedb #Conta qual a tarefa que ele vai fazer
+--threads 4 #Número de processadores
+--in /databases/big/ref/uniprotkb_proteome_UP000000625_2024_08_29.fasta.gz #Onde está o arquivo do UniProt com as sequências
+--db ecoli #Nome que quer dar ao seu banco de dados
+
+

Após formatarmos o banco de dados, podemos chamar o alinhamento. Temos um banco de dados de sequências de aminoácidos (UniProt), e arquivo de entrada de aminoácidos ou nucleotídeos. Poderiam então escolher entre duas opções:

+ +

Para otimizar o tempo - dado que temos as sequências de aminoácidos - rodaremos o blastp:

+
diamond blastp -d ecoli.dmnd -q 3.Genes/lib1.aa -o 4.Function/lib1.diamond -e 1e-20 -k 1 -p 4 --outfmt 6
+
+

As flags são:

+
diamond #Chama o programa
+blastp #Diz a opção de alinhamento
+-d ecoli.dmnd #Aponta o banco de dados
+-q 3.Genes/lib1.aa #Aponta o arquivo de entrada 
+-o 4.Function/lib1.diamond #Define o nome do arquivo de saída
+-e 1e-20 #Estabelece um mínimo de e-value
+-k 1 #Pede para mostrar só o melhor resultado para cada um dos meus genes
+-p 4 #Usa 4 processadores
+--outfmt 6 #Mostra o resultado em um formato tabular com 12 colunas
+
+

Podemos ver como é o resultado:

+
head 4.Function/lib1.diamond
+
+#Descrição das colunas
+ # 1.  qseqid      query or source (gene) sequence id
+ # 2.  sseqid      subject or target (reference genome) sequence id
+ # 3.  pident      percentage of identical positions
+ # 4.  length      alignment length (sequence overlap)
+ # 5.  mismatch    number of mismatches
+ # 6.  gapopen     number of gap openings
+ # 7.  qstart      start of alignment in query
+ # 8.  qend        end of alignment in query
+ # 9.  sstart      start of alignment in subject
+ #10.  send        end of alignment in subject
+ #11.  evalue      expect value
+ #12.  bitscore    bit score
+
+

Vemos que a segunda coluna contém o nome da proteína que melhor descreve nosso gene no banco de dados. +Agora queremos adicionar informação funcional à ela. +Vamos adicionar a informação do KEGG, e saber qual KO e via metabólica nosso gene pertence.

+

Para isso vamos executar uma sério de comandos.

+
#Entrar na pasta que contém os resultados
+cd 4.Function
+
+#Pega o nome de cada proteína que teve hit
+awk -F'|' '{print $2}' lib1.diamond > protein_ids.txt
+head protein_ids.txt
+
+#Adicionar o nome dessa proteína do UniProte à nossa tabela do resultado do Diamond
+paste lib1.diamond protein_ids.txt > diamond_with_protein_ids.txt
+
+# Colocar em ordem alfabética as colunas que contem as informações que vamos precisar.
+## A coluna que contém o nome das proteínas é a 13.
+sort -k13 diamond_with_protein_ids.txt > sorted_diamond.txt
+## A coluna que contém os nomes das mesmas proteínas no arquivo de conversão, é a 1
+sort -k1 /databases/big/ref/Uniprot2KEGG.eco.txt > sorted_Uniprot2KEGG.txt
+## Juntamos as duas tabelas, levando em consideração o ponto de junção como sendo a 13 coluna do resultado do Diamond, com a 1 coluna da tabela de conversão.
+join -t$'\t' -1 13 -2 1 sorted_diamond.txt sorted_Uniprot2KEGG.txt > diamond_with_KEGG.txt
+
+# Colocar em ordem alfabética as colunas que contem as informações que vamos precisar.
+## A coluna que contém o nome da entrada no KEGG é a 14.
+sort -k14 diamond_with_KEGG.txt > sorted_diamond_with_KEGG.txt
+## A coluna que contém os nomes das mesmas proteínas no arquivo de conversão, é a 1
+sort -k1 /databases/big/ref/KEGG2KO.eco.txt > sorted_KEGG2KO.txt
+## Juntamos as duas tabelas, levando em consideração o ponto de junção como sendo a 14 coluna do resultado do Diamond, com a 1 coluna da tabela de conversão.
+join -t$'\t' -1 14 -2 1 sorted_diamond_with_KEGG.txt sorted_KEGG2KO.txt > diamond_with_KO.txt
+
+# Temos uma tabela com 15 colunas, sendo que a última delas é o identificador do KEGG Orthology
+## Podemos ver qual é o KO mais representado
+cat diamond_with_KO.txt | cut -f15 | sort | uniq -c | sort -k1 -nr | head
+## Podemos gerar uma lista para visualizar no mapa metabólico
+cat diamond_with_KO.txt | cut -f15 | sort -u | sed -e 's/ko://g' > KO_list.txt
+
+

Visualização em vias metabólicas

+

A partir da lista de KOs, é possível ver suas funções no site do KEGG ou criar uma via metabólica interativa com o iPath3.

+
+
+ + + + diff --git a/elab/tutorial/map.html b/elab/tutorial/map.html new file mode 100644 index 0000000..b7a3f9b --- /dev/null +++ b/elab/tutorial/map.html @@ -0,0 +1,406 @@ + + + + + Aula Mapeamento + + + + + + + + + + + + + + +
+

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:

+
conda activate mapping
+
+
+

Escolha sua amostra

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
AmostrasAmostrasAmostrasAmostrasAmostras
IRR1820_S1IRR1821_S2IRR1822_S3IRR1823_S4IRR1824_S5
IRR1825_S6IRR1826_S7IRR1827_S8IRR1828_S9IRR1829_S10
IRR1830_S11IRR1831_S12IRR1832_S13IRR1833_S14IRR1834_S15
IRR1835_S16IRR1836_S17IRR1837_S18IRR1838_S19IRR1839_S20
IRR1840_S21IRR1841_S22IRR1842_S23IRR1843_S24IRR1844_S25
IRR1845_S26IRR1846_S27IRR1847_S28IRR1848_S29IRR1849_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:/databases/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 /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
+
+

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 /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
+
+

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 /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
+
+
+

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.

+
conda activate pangolin 
+pangolin variant.fa
+conda 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
+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.
+
+

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 /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
+
+

Mapeamento

+

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.

+

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 /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
+
+
+

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 /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
+
+ + + + + \ No newline at end of file