Relatório de Bioinformática aplicada II: Análise de transcriptomas

# Relatório de Bioinformática aplicada II: Análise de transcriptomas ### Docente: Dr. Daniel Guariz Pinheiro ### Discente: Daiane Silva Bonaldi # 1. Simulação A simulação a seguir, necessita do pacote E-Direct com as ferramentas E-utilities. ## 1.1 Requisitos Seguir as instruções para instalação do e-direct. ``` cd ~ /bin/bash perl -MNet::FTP -e \ '$ftp = new Net::FTP("ftp.ncbi.nlm.nih.gov", Passive => 1); $ftp->login; $ftp->binary; $ftp->get("/entrez/entrezdirect/edirect.tar.gz");' gunzip -c edirect.tar.gz | tar xf - rm edirect.tar.gz builtin exit export PATH=${PATH}:$HOME/edirect >& /dev/null || setenv PATH "${PATH}:$HOME/edirect" ./edirect/setup.sh ``` Adicionar o caminho no $PATH caso não tenha solicitado que seja feito automaticamente no final do passo anterior ``` echo "export PATH=\$PATH:\$HOME/edirect" >> $HOME/.bash_profile ``` ## 1.2. Construção do transcriptoma para simulação Teremos que utilizar um transcriptoma com um genoma de referência. O relatório utilizará genes de Apis mellifera. A seguir vamos utilizar as ferramentas esearch/efetch para fazer o download e construir o arquivo fasta contendo os transcritos referência para a simulação dos experimentos de investigação do transcriptoma. Para conhecer os IDs de interesse, devemos consultar a página do banco de dados Gene do NCBI. #### Melt melittin [ Apis mellifera (honey bee) ] Gene ID: 406130 mRNA and Protein(s) NM_001011607.2 :::success Protein: Melittin Gene: MELT Melitina: Principal toxina do veneno de abelha com forte atividade hemolítica e atividade antimicrobiana. Produz dor ativando direta e indiretamente as células nociceptoras primárias, devido à sua capacidade de ativar a fosfolipase A2 da membrana plasmática e sua atividade de formação de poros. ::: #### Hbg3 alpha-glucosidase [ Apis mellifera (honey bee) ] Gene ID: 406131 mRNA and Protein(s) NM_001011608.1 :::success Protein Submitted name: Alpha-glucosidase Gene: hbg3 Funções moleculares: atividade catalítica. Processos biológicos: processo metabólico de carboidratos. ::: #### Lop1 long wavelength sensitive opsin 1 [ Apis mellifera (honey bee) ] Gene ID: 413961 mRNA and Protein(s) NM_001011639.2 :::success Protein Submitted name: Long wavelength sensitive opsin 1 Gene: Lop1 Funções moleculares: atividade do receptor acoplado à proteína. atividade de fotorreceptores. Processos biológicos: fonte de fototransdução. ligação proteína-cromóforo. ::: #### Fem feminizer [ Apis mellifera (honey bee) ] Gene ID: 724970 mRNA and Protein(s) NM_001134828.1 :::success Protein Submitted name: Feminizer Gene: fem Acredita-se que o gene fem esteja envolvido na regulação do processamento dos transcritos do gene amdsx. Em insetos; o pré-mRNA de dsx sofre processamento alternativo nos sexos que determina variações na estrutura do fator de transcrição DSX e consequente regulação diferencial em machos e fêmeas. ::: #### Ac3 adenylate cyclase 3 [Apis mellifera (honey bee) ] Gene ID: 408363 mRNA and Protein(s) NM_001077808.1 :::success Gene: Ac3 o gene pertence a família de proteínas adenilato ciclase, está relacionado com a homeostase energética. ::: ## 1.3. Criando os arquivos de abundância Os arquivos abundance_A.txt e abundance_B.txt devem ser criados manualmente, respectivamente para a condição biológica A e condição biológica B. ``` abundance_A.txt NM_001011607.2 0.50 NM_001011608.1 0.05 NM_001011639.2 0.05 NM_001134828.1 0.25 NM_001077808.1 0.15 ``` ``` abundance_B.txt NM_001011607.2 0.15 NM_001011608.1 0.05 NM_001011639.2 0.45 NM_001134828.1 0.25 NM_001077808.1 0.10 ``` ## 1.4 Gerando os arquivos .fastq Gerando a primeira réplica da simulação da corrida para a amostra A utilizando o (abundance_A.txt). ### Passo 1 Gerar “25000” sequências de fragmentos aleatórios em média de “300 pb” com desvio padrão de “30 pb” a partir das sequências do transcriptoma montado anteriormente “transcriptoma.fa”, considerando a abundância contida em “abundance_A.txt” ``` generate_fragments.py -r transcriptoma.fa \ -a ./abundance_A.txt \ -o ./tmp.frags.r1 \ -t 25000 \ -i 300 \ -s 30 ``` ### Passo 2 Renomear as sequências para “SAMPLEA” juntamente com um número hexadecimal sequencial começando por “0000000001” e posicionar todas as bases em uma única linha (até 1000 pb). A linha de SED irá substituir a descrição adicional. ``` cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \ -if FASTA \ -of FASTA \ -p SAMPLEA1 \ -w 1000 | \ sed 's/^>$\S\+$.*/>\1/' \ > ./frags_A1.fa ``` ### Passo 3 Fazer a simulação considerando leituras “paired” 2x “151” bp , com parâmetros configurados usando simNGS: ``` cat ./frags_A1.fa | simNGS -a \ AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ -p paired \ /usr/local/bioinfo/simNGS/data/s_4_0099.runfile \ -n 151 > ./SAMPLEA1.fastq 2> SAMPLEA1.err.txt ``` ### Passo 4 Desintercalar as leituras em dois arquivos SAMPLEA_R1.fastq e SAMPLEA_R2.fastq no diretório …/raw ``` mkdir -p ./raw deinterleave_pairs SAMPLEA1.fastq \ -o ./raw/SAMPLEA1_R1.fastq \ ./raw/SAMPLEA1_R2.fastq ``` ### Passo 5 Remover arquivos desnecessários: `rm -f ./tmp.frags.r1.1.fasta ./frags_A1.fa ./SAMPLEA1.fastq ./SAMPLEA1.err.txt` ### Passo 6 Verificar o número de sequências geradas ``` wc -l ./raw/SAMPLEA1_R1.fastq wc -l ./raw/SAMPLEA1_R2.fastq ``` ### Repetir o mesmo processo para a réplica 2 da condição biológica A (SAMPLEA2) utilizando o arquivo abundance_A.txt ``` generate_fragments.py -r transcriptoma.fa \ -a ./abundance_A.txt \ -o ./tmp.frags.r1 \ -t 25000 \ -i 300 \ -s 30 ``` ``` cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \ -if FASTA \ -of FASTA \ -p SAMPLEA2 \ -w 1000 | \ sed 's/^>$\S\+$.*/>\1/' \ > ./frags_A2.fa ``` ``` cat ./frags_A2.fa | simNGS -a \ AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ -p paired \ /usr/local/bioinfo/simNGS/data/s_4_0099.runfile \ -n 151 > ./SAMPLEA2.fastq 2> SAMPLEA2.err.txt ``` ``` mkdir -p ./raw deinterleave_pairs SAMPLEA1.fastq \ -o ./raw/SAMPLEA2_R1.fastq \ ./raw/SAMPLEA2_R2.fastq ``` ``` rm -f ./tmp.frags.r1.1.fasta ./frags_A1.fa ./SAMPLEA2.fastq ./SAMPLEA2.err.txt ``` ``` wc -l ./raw/SAMPLEA2_R1.fastq wc -l ./raw/SAMPLEA2_R2.fastq ``` ### Repetir o mesmo processo para a réplica 1 e réplica 2 para a condição biológica B (SAMPLEB) ## réplica 1 ``` generate_fragments.py -r transcriptoma.fa \ -a ./abundance_B.txt \ -o ./tmp.frags.r1 \ -t 25000 \ -i 300 \ -s 30 ``` ``` cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \ -if FASTA \ -of FASTA \ -p SAMPLEB1 \ -w 1000 | \ sed 's/^>$\S\+$.*/>\1/' \ > ./frags_B1.fa ``` ``` cat ./frags_B1.fa | simNGS -a \ AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ -p paired \ /usr/local/bioinfo/simNGS/data/s_4_0099.runfile \ -n 151 > ./SAMPLEB1.fastq 2> SAMPLEB1.err.txt ``` ``` mkdir -p ./raw deinterleave_pairs SAMPLEB1.fastq \ -o ./raw/SAMPLEB1_R1.fastq \ ./raw/SAMPLEB1_R2.fastq ``` ``` rm -f ./tmp.frags.r1.1.fasta ./frags_B1.fa ./SAMPLEB1.fastq ./SAMPLEB1.err.txt ``` ``` wc -l ./raw/SAMPLEB1_R1.fastq wc -l ./raw/SAMPLEB1_R2.fastq ``` ## réplica 2 ``` generate_fragments.py -r transcriptoma.fa \ -a ./abundance_B.txt \ -o ./tmp.frags.r1 \ -t 25000 \ -i 300 \ -s 30 ``` ``` cat ./tmp.frags.r1.1.fasta | renameSeqs.pl \ -if FASTA \ -of FASTA \ -p SAMPLEB2 \ -w 1000 | \ sed 's/^>$\S\+$.*/>\1/' \ > ./frags_B2.fa ``` ``` cat ./frags_B2.fa | simNGS -a \ AGATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCTCGTATGCCGTCTTCTGCTTG:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT \ -p paired \ /usr/local/bioinfo/simNGS/data/s_4_0099.runfile \ -n 151 > ./SAMPLEB2.fastq 2> SAMPLEB2.err.txt ``` ``` mkdir -p ./raw deinterleave_pairs SAMPLEB2.fastq \ -o ./raw/SAMPLEB2_R1.fastq \ ./raw/SAMPLEB2_R2.fastq ``` `rm -f ./tmp.frags.r1.1.fasta ./frags_B2.fa ./SAMPLEB2.fastq ./SAMPLEB2.err.txt` ``` wc -l ./raw/SAMPLEB2_R1.fastq wc -l ./raw/SAMPLEB2_R2.fastq ``` ## 2. Obtenção das referências para análises ### 2.1.Sequências e Anotações #### A) Obtendo as sequências no formato FASTA e anotações gênicas no formato GFF. mkdir ./refs cd ./refs #### B) Obtendo arquivo fna `wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/254/395/GCF_003254395.2_Amel_HAv3.1/GCF_003254395.2_Amel_HAv3.1_genomic.fna.gz` #### C) Obtendo arquivo gff `wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/254/395/GCF_003254395.2_Amel_HAv3.1/GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz` #### D) Para descompactar `gunzip GCF_003254395.2_Amel_HAv3.1_genomic.fna.gz` `gunzip GCF_003254395.2_Amel_HAv3.1_genomic.gff.gz` ## 3. Limpeza dos arquivos referência Voltar para diretório ./refs: `cd ../refs/` ### 3.1. Limpeza dos cabeçalhos FASTA Criar arquivo com o script (cleanfasta.sh) abaixo: ``` nano cleanfasta.sh ``` ``` #!/bin/bash infile=$1 if [ ! ${infile} ]; then echo "Missing input fasta file" exit fi if [ ! -e ${infile} ]; then echo "Not found input fasta file (${infile})" exit fi sed 's/^>$\S\+$.*/>\1/' ${infile} ``` Dar permissão: `chmod a+x cleanfasta.sh` Executar: ``` ./cleanfasta.sh GCF_003254395.2_Amel_HAv3.1_genomic.fna > genome.fa ``` ## 3.2. Ajustes do arquivo GFF: Executar o script fixNCBIgff.sh do repositório Github: `fixNCBIgff.sh GCF_003254395.2_Amel_HAv3.1_genomic.gff genome.gff` ## Testes com GFF #### Para a conversão de GFF versão 3 para um arquivo GTF: `gffread genome.gff -g genome.fa -T -o genome.gtf` #### Para obtenção de um arquivo com a sequência FASTA do trasncriptoma: `gffread genome.gff -g genome.fa -w transcriptome.fa` #### Para obtenção de um arquivo com a sequência FASTA do proteoma: `gffread genome.gff -g genome.fa -y proteome.fa` ## 3.3. Testes de alinhamento com Bowtie2 :::info O Bowtie 2 é uma ferramenta para o alinhamento de leituras de sequenciamento relativamente curtas para genomas longos. Dito isso, ele lida com sequências de referência arbitrariamente pequenas (por exemplo, amplicons) e leituras muito longas (ou seja, mais de 10s ou 100s de kilobases), embora seja mais lento nessas configurações "Alinhamento" é o processo pelo qual descobrimos como e onde as seqüências de leitura são semelhantes à sequência de referência. Um "alinhamento" é o resultado desse processo, especificamente: um alinhamento é uma maneira de "alinhar" alguns ou todos os caracteres na leitura com alguns caracteres da referência, de maneira a revelar como eles são semelhantes. Usamos o alinhamento para adivinhar qual a origem da leitura em relação ao genoma de referência, mas nem sempre é possível determinar isso com certeza. O Bowtie 2 gera alinhamentos no formato SAM, permitindo a interoperação com um grande número de outras ferramentas (por exemplo, SAMtools, GATK) que usam SAM. O Bowtie 2 é distribuído sob a licença GPLv3 e é executado na linha de comando no Windows, Mac OS X e Linux. ::: #### Indexação do genoma: ``` cd ./refs bowtie2-build -f genome.fa genome ``` :::warning `---bowtie2-build \` cria um índice Bowtie a partir de um conjunto de sequências de DNA.Esses arquivos juntos constituem o índice: são tudo o que é necessário para alinhar as leituras a essa referência. Os arquivos FASTA da sequência original não são mais usados pelo Bowtie 2 depois que o índice é criado. `--f \` Os arquivos de entrada de referência (especificados como <reference_in>) são arquivos FASTA (geralmente com extensão .fa, .mfa, .fna ou similar). ::: #### Alinhamento ``` cd ../ time bowtie2 -x ./refs/genome -q -1 ./raw/SAMPLEA1_R1.fastq -2 ./raw/SAMPLEA1_R2.fastq -S bowtie2-test.sam ``` #### Conversão de SAM para BAM #### A) Convertendo SAM para BAM com samtools "view" ``` samtools view -b bowtie2-test.sam > bowtie2-test.bam ``` :::info Para fazer algo significativo com os dados de alinhamento dos alinhadores (que produzem saída SAM com base em texto), precisamos primeiro converter o SAM em sua contraparte binária, o formato BAM. Para converter SAM em BAM, usamos o comando samtools view. Devemos especificar que nossa entrada está no formato SAM (por padrão, espera BAM) usando a opção -S. Também devemos dizer que queremos que a saída seja BAM (por padrão, ela produz BAM) com a opção -b. O Samtools segue a convenção do UNIX de enviar sua saída para o UNIX STDOUT, portanto, precisamos usar um operador de redirecionamento (">") para criar um arquivo BAM a partir da saída. ::: #### B) samtools “sort” ``` samtools sort bowtie2-test.bam -o bowtie2-test-sorted.bam ``` :::info Quando você alinha arquivos FASTQ com todos os alinhadores de sequência atuais, os alinhamentos produzidos são em ordem aleatória em relação à sua posição no genoma de referência. Em outras palavras, o arquivo BAM está na ordem em que as seqüências ocorreram nos arquivos FASTQ de entrada. ::: #### C) samtools “index” ``` samtools index bowtie2-test-sorted.bam ``` :::info A indexação de um arquivo BAM classificado por genoma permite extrair rapidamente alinhamentos sobrepondo regiões genômicas específicas. Além disso, a indexação é exigida pelos visualizadores de genoma, como IGV, para que os visualizadores possam exibir rapidamente alinhamentos em cada região genômica para a qual você navega. ::: ### Observação: Devem ser copiados os seguintes arquivos para serem abertos no IGV: bowtie2-test-sorted.bam bowtie2-test-sorted.bam.bai Mapeamento de RNA-Seq reads x Genoma: :::success ### Resultados do alinhamento com Bowtie2: ::: A visualização do resultado dos alinhamentos foi realizada com o programa Integrative Genomics Viewer (IGV). :::info O Integrative Genomics Viewer (IGV) é uma ferramenta de visualização de alto desempenho para exploração interativa de grandes conjuntos de dados genômicos integrados. Ele suporta uma ampla variedade de tipos de dados, incluindo dados de sequência baseados em array e de próxima geração, e anotações genômicas. ::: ### Melt melittin Gene ID: 406130 mRNA and Protein(s) NM_001011607.2 ![](https://i.imgur.com/Z3axCux.png) ### Hbg3 alpha-glucosidase Gene ID: 406131 mRNA and Protein(s) NM_001011608.1 ![](https://i.imgur.com/aceBKFl.png) ### Lop1 long wavelength sensitive opsin 1 Gene ID: 413961 mRNA and Protein(s) NM_001011639.2 ![](https://i.imgur.com/7Yve6im.png) ### Fem feminizer Gene ID: 724970 mRNA and Protein(s) NM_001134828.1 ![](https://i.imgur.com/Mtvx7jd.png) ### Ac3 adenylate cyclase 3 Gene ID: 408363 mRNA and Protein(s) NM_001077808.1 ![](https://i.imgur.com/muUVhQb.png) ## Montando o script de pré-processamento final A ordem de execução será a seguinte: FastQC (pré) -> ATROPOS (insert) -> ATROPOS (adater) -> resultados concatanados ATROPOS (insert) & ATROPOS (adapter) -> PrinSeq -> FastQC (pós) ## 4. Montando o script de análise RNA-Seq com referência 1. O script rnaseq-ref.sh foi criado utilizando o script preprocess3.sh para fazer o pré-processamento e posteriormente executar o STAR com as leituras Paired-End e as leituras Single-End (singletons), além de combinar os resultados, ordenar e executar o script SAM_nameSorted_to_uniq_count_stats.pl para coletar as estatísticas de alinhamento. ## 4.1.preprocess3.sh Link: [Script preprocess3.sh](https://github.com/dbonaldi/RelatorioBio/blob/master/preprocess3.sh) ``` wget https://github.com/dbonaldi/RelatorioBio/blob/master/preprocess3.sh ``` ## 4.2.rnaseq-ref.sh Link: [Script rnaseq-ref.sh ](https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-ref.sh) ``` wget https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-ref.sh ``` 3. Para verificar que o script foi executado com sucesso, confira com a listagem dos arquivos resultantes dos alinhamentos e que foram posteriormente ordenados com samtools. `ls -lh ./output/star_out_final/*/Aligned.out.sorted.bam` 4. Verifique a porcentagem geral de alinhamento em pares de modo adequado (proper pairs): `find ./output/star_out_final/ -name 'Aligned.stats.txt' -exec bash -c 'echo $(basename $(dirname $0)) ":" $(grep "^Overall" $0)' {} \;` 5. Execução do pipeline O script rnaseq-ref.sh deve ser executado da seguinte forma: `./rnaseq-ref.sh ./raw ./output ./refs/genome.gtf ./refs/genome.fa` :::success #### RESULTADO DOS ALINHAMENTOS ::: ### Número total de reads dos arquivos simulados ``` find ./ -name ‘*_R[12].fastq’ -exec bash -c ‘echo -e “0\t(echo “$(cat $0 | wc -l)/4” | bc)”’ {} ; ``` ./SAMPLEA2_R1.fastq 25000 ./SAMPLEA1_R1.fastq 25000 ./SAMPLEB2_R1.fastq 25000 ./SAMPLEB2_R2.fastq 25000 ./SAMPLEB1_R1.fastq 25000 ./SAMPLEA2_R2.fastq 25000 ./SAMPLEB1_R2.fastq 25000 ./SAMPLEA1_R2.fastq 25000 ### Número de reads que foram alinhadas, mas com redundância `samtools view All.sorted.bam | cut -f1,3| sed ‘/*/d’ | wc -l` Foram alinhadas: 197210 reads. ### Número de reads que foram alinhadas, sem redundância `samtools view All.sorted.bam | cut -f1,3| sed ‘/*/d’ | sort -u| wc -l` Foram alinhadas: 97948 reads. ## 4.3. Acrescentando informações aos resultados de expressão gênica diferencial #### A) Para a consulta ao NCBI, vamos utilizar diretamente via linha de comando, as ferramentas E-utilities do NCBI. Será realizada a busca pelos genes de Apis mellifera utilizando o Taxonomy ID específico desta linhagem. Taxonomy ID: 7460 ``` esearch -db gene -query "7460" \ | efetch -format tabular > ./refs/gene_result.txt ``` :::success #### Resultado: ::: tax_id Org_name GeneID CurrentID Status Symbol Aliases description other_designations map_location chromosome genomic_nucleotide_accession.version start_position_on_the_genomic_accession end_position_on_the_genomic_accession orientation exon_count OMIM 9606 Homo sapiens 7460 7461 secondary WBSCR3 WSCR3 Williams-Beuren syndrome chromosome region 3 7q11.23 7 #### B) Selecionaremos somente as colunas 3, 6 e 8, respectivamente, GeneID, Symbol e description, para um novo arquivo ./refs/gene_info.txt. `cut -f 3,6,8 ./refs/gene_result.txt > ./refs/gene_info.txt` :::success #### Resultado: GeneID: 7460 Symbol: WBSCR3 description: Williams-Beuren syndrome chromosome region 3 ::: #### C) O cufflinks/cuffmerge, durante o seu processamento, pode fundir genes considerados inicialmente distintos, ou seja, que possuem IDs distintos (exemplo: M431DRAFT_12071 e M431DRAFT_502167). Isso pode ser observado no arquivo gene_exp.diff. Neste caso, o script mergeR.R, criado para realizar a tarefa final de fundir os dados, não irá funcionar, afinal, na linha correspondente, os identificadores estarão unidos por vírgula (exemplo: M431DRAFT_12071,M431DRAFT_502167) impossibilitando o relacionamento. Essa impossibilidade se dá pelo fato de que no arquivo das descrições os dois identificadores estarão cada um em uma linha distinta no arquivo ./refs/gene_info.txt. ``` splitteR.R --x="./output/cuffdiff/gene_exp.diff" \ --col.x="gene" \ --by.x="," \ --out="./output/cuffdiff/gene_exp.diff.splitted.txt" ``` :::info ### SplitteR.R Dividir sequência (s) em seqüências menores. O divisor divide uma ou mais seqüências de entrada em subsequências menores, opcionalmente sobrepostas. O tamanho da subsequência e a sobreposição (se houver) podem ser especificados. ::: #### D) Após gerar o novo arquivo, executar a tarefa de fusão dos dados tabulares. Note que utilizaremos o parâmetro –all.x para que todas as linhas do arquivo ./output/cuffdiff/gene_exp.diff.spplitted.txt estejam no resultado final, no arquivo ./output/cuffdiff/gene_exp.diff.spplitted.merged.txt, mesmo que não tenha uma linha correspondente em ./refs/gene_info.txt. Essa última situação pode acontecer quando o cufflinks cria um novo mapeamento gênico, que não corresponde a nenhum outro contido no GFF/GTF original, usado como guia para alinhamento e montagem do transcritoma. ``` mergeR.R --x="./output/cuffdiff/gene_exp.diff.spplitted.txt" \ --all.x \ --y="./refs/gene_info.txt" \ --by.x="gene" \ --by.y="Symbol" \ --out="./output/cuffdiff/gene_exp.diff.spplitted.merged.txt" ``` ::: info ### mergeR.R Tem como função mesclar duas seqüências sobrepostas. A fusão lê duas seqüências de entrada sobrepostas do mesmo tipo (normalmente nucleotídeo) e usa um algoritmo de alinhamento global (Needleman & Wunsch) para alinhar as seqüências de maneira ideal. Uma sequência mesclada é gerada a partir do alinhamento e gravada no arquivo de saída. O alinhamento da sequência também é gravado. ::: :::success #### Resultado dos genes diferentemente expressos que foram significativos: ::: MSTRG.1 MSTRG.1 Fem NC_037640.1:11792750-11800256 SAMPLEA SAMPLEB OK 58987.6 192804 1.70865 48.3841 5e-05 5e-05 yes MSTRG.10 MSTRG.10 Ac3 NC_037648.1:12956996-12973602 SAMPLEA SAMPLEB OK 12421.4 24053.2 0.953396 21.0613 5e-05 5e-05 yes MSTRG.11 MSTRG.11 Lop1 NC_037652.1:7414821-7417248 SAMPLEA SAMPLEB OK 8554.44 248052 4.85782 100.816 5e-05 5e-05 yes MSTRG.4 MSTRG.4 Melt NC_037641.1:12394665-12411860 SAMPLEA SAMPLEB OK 3.88417e+06 1.7656e+06 -1.13745 -32.6396 5e-05 5e-05 yes MSTRG.8 MSTRG.8 Hbg3 NC_037643.1:4479593-4655231 SAMPLEA SAMPLEB OK 15534.2 34983 1.1712 17.9086 5e-05 5e-05 yes :::danger ## Descrição dos softwares utilizados para a montagem dos scripts para a análise do RNA-Seq com referência ::: ## FastqC :::info O FastQC é um aplicativo de controle de qualidade que permite aos usuários executar inúmeras verificações de controle de qualidade em dados brutos de sequência gerados por pipelines de sequenciamento de alto rendimento, como as plataformas Illumina e ABI SOLiD no formato FASTQ. Ele gera como saída um relatório abrangente de várias páginas sobre a composição e a qualidade das leituras no formato HTML, com uma página para cada uma das leituras (por exemplo, extremidade única, extremidade emparelhada: frente e extremidade emparelhada: reversa). ::: ## Atropos :::info O Atropos é uma ferramenta de corte e leitura que possui recursos para opções de corte específicas para Methyl-Seq, detecção automatizada de adaptadores, estimativa de erro de sequenciamento, cálculo de métricas de controle de qualidade antes e depois do corte e suporte para dados gerados por muitos métodos de sequenciamento. O software inclui um comando que fornece uma estimativa da taxa de erro em cada arquivo de entrada e opções para coletar métricas de controle de qualidade (CQ) antes e / ou após o corte. ::: :::warning ### Parâmetros utilizados com o Atropos `--e ` O nível de tolerância a erros é ajustado especificando uma taxa de erro máxima, que é 0,1 (= 10%) por padrão. Usa-se a opção -e para definir um valor diferente. Para determinar o número de erros permitidos, a taxa máxima de erros é multiplicada pela duração da correspondência. `--aligner insert ` Algorítmo de alinhamento baseado em inserção (mais preciso) `--aligner adapter ` Algorítmo de alinhamento baseado em adaptadores `--e ` O nível de tolerância a erros é ajustado especificando uma taxa de erro máxima, que é 0,1 (= 10%) por padrão. Usa-se a opção -e para definir um valor diferente. Para determinar o número de erros permitidos, a taxa máxima de erros é multiplicada pela duração da correspondência. `--n ` Por padrão, no máximo uma sequência do adaptador é removida de cada leitura, mesmo que várias seqüências do adaptador tenham sido fornecidas. Isso pode ser alterado usando a opção --times (ou sua forma abreviada -n). O Atropos procurará todas as seqüências de adaptadores fornecidas repetidamente, até que nenhuma correspondência de adaptador tenha sido encontrada ou até que o número especificado de rodadas seja atingido. `--m ` Exclui as leituras processadas com menos de N bases. `--op-order GAWCQ ` Pode-se ter algum controle sobre a ordem em que as operações são aplicadas. A opção --op-order recebe uma sequência de caracteres que controla a ordem na qual as quatro primeiras operações são aplicadas. A = adapter trimming; C = cutting (unconditional); G = NextSeq trimming; Q = quality trimming; W = overwrite poor quality reads. `--match-read-wildcards ` Use --match-read-wildcards para ativar curingas também em reads. Todos os códigos de nucleotídeos (caracteres curinga) são suportados. Por exemplo, use um N na sequência do adaptador para corresponder a qualquer nucleotídeo na leitura ou use -a YACGT para um adaptador que corresponda a CACGT e TACGT. O caractere curinga N é útil para aparar adaptadores com um código de barras variável incorporado. `-O ` Usado para controlar correspondências aleatórias especificando uma sobreposição mínima entre o adaptador e a leitura. O comprimento mínimo de sobreposição pode ser alterado com o parâmetro --overlap (ou sua versão curta -O). Correspondências mais curtas são simplesmente ignoradas e as bases não são cortadas. `-q ` O parâmetro -q (ou --quality-cutoff) pode ser usado para aparar extremidades de baixa qualidade das leituras antes da remoção do adaptador. Para que isso funcione corretamente, os valores de qualidade devem ser codificados como ascii (qualidade phred + 33). Se eles são codificados como ascii (qualidade phred + 64), você precisa adicionar --quality-base = 64 à linha de comando. `-T ` A opção -T (ou --threads) permite especificar o número de threads a serem usados para ler entradas, processar leituras e gravar saídas `--correct-mismatches ` Os pares de leitura que se sobrepõem são derivados da mesma sequência e, portanto, espera-se que sejam idênticos. No entanto, todas as tecnologias de seqüenciamento estão associadas a algum nível de erro que pode levar a incompatibilidades. O alinhador de insertos pode corrigir incompatibilidades entre leituras sobrepostas usando uma das três estratégias . `--pair-filter any ` Parâmetro de filtragem das reads `-a ` A sequência do adaptador é fornecida com a opção -a. Obviamente, você precisa substituir o AACCGGTT pela sua sequência real do adaptador. As leituras são lidas no arquivo de entrada input.fastq e gravado no arquivo de saída output.fastq `-A ` Parâmetro usado para remoção de adaptadores. `-o /` Especificar o arquivo de saída. `-p ` Para leituras de Paired-end, a segunda leitura em um par é sempre gravada no arquivo especificado pela opção -p. `-pe1 e -pe2 ` Os arquivos Paired-end FASTQ são suportados usando -pe1 e -pe2 para especificar os arquivos de read 1 e read 2, respectivamente. `--untrimmed-output ` Escreve todas as reads sem adaptadores em FILE (no formato FASTA / FASTQ) em vez de gravá-las no arquivo de saída regular. `--untrimmed-paired-output ` Usado junto com --untrimmed-output. A segunda leitura em um par é gravada nesse arquivo quando o par processado não foi aparado. ::: ## Prinseq :::info O PRINSEQ é uma ferramenta que gera estatísticas resumidas dos dados de sequência e qualidade e é usada para filtrar, reformatar e aparar dados de sequenciamento de nova geração. Foi especialmente desenvolvido para dados 454 / Roche, mas também pode ser usado para outros tipos de dados de sequência. ::: :::warning ### Parâmetros utilizados com o Prinseq `--fastq ` Arquivo de entrada no formato FASTQ que contém os dados de sequência e qualidade. Use stdin em vez de um nome de arquivo para ler em STDIN (-fasta stdin). Isso pode ser útil para processar arquivos compactados usando pipes Unix. `--fastq2 ` Apenas para dados de extremidade emparelhada. Arquivo de entrada no formato FASTQ que contém os dados de sequência e qualidade. Os identificadores de sequência para duas sequências de extremidade emparelhadas correspondentes em arquivos separados podem ser marcados por / 1 e / 2, ou _L e _R, ou _left e _right, ou devem ter o mesmo identificador exato nos dois arquivos de entrada. As sequências de entrada devem ser classificadas por seus identificadores de sequência. Singletons são permitidos nos arquivos de entrada. `--out_format 3 ` Formato de saída 1 (apenas FASTA) 2 (FASTA e QUAL) 3 (FASTQ) 4 (FASTQ e FASTA) 5 (FASTQ, FASTA e QUAL) ` --trim_qual_window 3 ` O tamanho da “janela deslizante” usada para calcular o índice de qualidade por tipo. Para parar na primeira base que falha na regra definida, use o tamanho da janela 1. `--trim_qual_step 1 ` Tamanho do passo usado para mover a “janela deslizante”. Para mover a janela sobre todos os índices de qualidade sem perder nenhum, o tamanho da etapa deve ser menor ou igual ao tamanho da janela. `--trim_qual_right 30 ` Corte a sequência pelo índice de qualidade do final 3 'com este limiar. `--trim_qual_type mean ` Tipo de cálculo do índice de qualidade a ser usado. `--trim_qual_rule lt ` Regra a ser usada para comparar o índice de qualidade com o valor calculado. As opções permitidas são lt (menor que), gt (maior que) e et (igual a). `--out_good ` Por padrão, os arquivos de saída são criados no mesmo diretório que o arquivo de entrada que contém os dados da sequência com um "_prinseq_good_XXXX" adicional em seu nome. Para alterar o nome e o local da saída, especifique o nome do arquivo usando esta opção. `--out_bad null ` Por padrão, os arquivos de saída são criados no mesmo diretório que o arquivo de entrada que contém os dados da sequência com um "_prinseq_bad_XXXX" adicional em seu nome (onde XXXX é substituído por caracteres aleatórios para evitar a substituição de arquivos anteriores). Para alterar o nome e o local da saída, especifique o nome do arquivo usando esta opção. ` --lc_method dust ` Método para filtrar sequências de baixa complexidade. `--lc_threshold 30 ` O valor limite usado para filtrar sequências por complexidade de sequência. O método de poeira usa isso como pontuação máxima permitida e o método de entropia como valor mínimo permitido. `--min_len 20 ` Sequência de filtro menor que min_len. `--trim_tail_right 5 ` Apare a cauda de poli-A / T com um comprimento mínimo de trim_tail_right na extremidade 3'. `--trim_tail_left 5 ` Apare a cauda de poli-A / T com um comprimento mínimo de trim_tail_left na extremidade 5 '. ` --ns_max_p 80 ` Sequência de filtros com mais de ns_max_p porcentagem de Ns. ` --noniupac ` Filtrar sequência com caracteres diferentes de A, C, G, T ou N. ::: ## STAR :::info O alinhamento ou mapeamento de leitura é realizado para determinar de onde no genoma as leituras se originaram. O processo de alinhamento consiste em escolher um genoma de referência apropriado para mapear nossas leituras e executar o alinhamento de leitura usando uma das várias ferramentas de alinhamento com reconhecimento de emenda, como STAR ou HISAT2. A escolha do alinhador geralmente é uma preferência pessoal e também depende dos recursos computacionais disponíveis. O pacote de software STAR executa essa tarefa com altos níveis de precisão e velocidade. Além de detectar junções de emenda anotadas e novas, o STAR é capaz de descobrir arranjos de sequência de RNA mais complexos, como RNA quimérico e circular. O STAR pode alinhar seqüências emendadas de qualquer tamanho com taxas de erro moderadas, fornecendo escalabilidade para tecnologias emergentes de sequenciamento. O programa também gera arquivos de saída que podem ser usados para muitas análises posteriores, como quantificação de transcrição / expressão de genes, expressão diferencial de genes, nova reconstrução de isoformas e visualização de sinais. ::: :::warning ### Parâmetros utilizados com o STAR `--genomeSAindexNbases 12 ` Comprimento (bases) da sequência de pré-indexação do SA. Normalmente entre 10 e 15. (Seqüências mais longas usarão muito mais memória, mas permitirão pesquisas mais rápidas.) Para pequenos genomas, o parâmetro --genomeSAindexNbases deve ser reduzido, com um valor típico de min (14, log2 (GenomeLength) / 2 - 1). `--sjdbOverhang 149 ` Especifica o comprimento da sequência genômica em torno da junção anotada a ser usada na construção do banco de dados de junções de emenda. Idealmente, esse comprimento deve ser igual ao ReadLength-1, em que ReadLength é o comprimento das leituras. `--runThreadN ` A opção define o número de encadeamentos a serem utilizados para a geração do genoma; ele deve ser definido como o número de núcleos disponíveis no nó do servidor. `--runMode ` A opção genomeGenerate direciona a STAR para executar o trabalho de geração de índices de genoma. `--genomeFastaFiles ` Especifica um ou mais arquivos FASTA com as sequências de referência do genoma. Múltiplas seqüências de referência são permitidas para cada arquivo fasta. `--genomeDir ` Especifica o caminho para o diretório (denominado "diretório genoma" onde os índices do genoma são armazenados. Esse diretório deve ser criado (com mkdir) antes da execução do STAR e precisa de permissões de gravação. `--sjdbGTFfile ` Especifica o caminho para o arquivo com transcrições anotadas no formato GTF padrão. O STAR extrairá junções de emenda deste arquivo e as utilizará para melhorar significativamente a precisão do mapeamento. Embora isso seja opcional e o STAR possa ser executado sem anotações, o uso de anotações é altamente recomendado sempre que elas estiverem disponíveis. ` --outSAMstrandField intronMotif ` Para dados de RNA-seq sem cadeia, Cufflinks / Cuffdiff exigem alinhamentos emendados com o atributo XS strand, que o STAR gerará com a opção introutMotif --outSAMstrandField. Conforme necessário, o atributo da cadeia XS será gerado para todos os alinhamentos que contêm junções de emenda. `--outFilterIntronMotifs ` Usado para remover as junções não canônicas das execuções de botão de punho. `--readFilesIn ` O (s) nome (s) (com caminho) dos arquivos que contêm as seqüências a serem mapeadas (por exemplo, arquivos FASTQ RNA-seq). Se estiver usando leituras de extremidade pareada Illumina, os arquivos read1 e read2 deverão ser fornecidos. O STAR pode processar arquivos FASTA e FASTQ. `--outFilterMultimapNmax 20 ` É o número máximo de alinhamentos múltiplos permitidos para uma leitura: se excedida, a leitura é considerada não mapeada `--outFileNamePrefix ` O STAR produz vários arquivos de saída. Todos os arquivos têm nome padrão, no entanto, você pode alterar os prefixos de arquivo usando --outFileNamePrefix /path/to/output/dir/prefix. Por padrão, esse parâmetro é ./, ou seja, todos os arquivos de saída são gravados no diretório atual. `--outSAMtype BAM Unsorted ` As extremidades emparelhadas de um alinhamento são sempre adjacentes, e vários alinhamentos de uma leitura também são adjacentes. Esse arquivo "não classificado" pode ser usado diretamente com software downstream, como o HTseq, sem a necessidade de classificação de nomes. A ordem das leituras corresponderá à dos arquivos FASTQ (A) de entrada apenas se um encadeamento for usado: --runThread 1 e --outFilterType --BySJout não é usado. ::: ## Cufflinks :::info O conjunto de ferramentas Cufflinks pode ser usado para realizar vários tipos diferentes de análises para experimentos de RNA-Seq. O conjunto de Cufflinks inclui vários programas diferentes que trabalham juntos para executar essas análises. Compreende os programas: 1.Cuffcompare 2.Cuffmerge 3.Cuffquant 4.Cuffdiff 5.Cuffnorm ::: :::warning ### Parâmetros utilizados com o Cufflinks `--g/–-GTF-guide ` Diz aos Cufflinks para usar a anotação de referência fornecida em um arquivo GFF para orientar a montagem do RABT. As transcrições de referência serão lado a lado com leituras falsas para fornecer informações adicionais na montagem. A saída incluirá todos os transcritos de referência, bem como quaisquer novos genes e isoformas reunidos. `--b/ --frag-bias-correct ` O fornecimento de cufflinks com um arquivo multifasta por meio dessa opção instrui-o a executar o algoritmo de detecção e correção de viés. `--u/ --multi-read-correct ` Diz aos Cufflinks para fazer um procedimento de estimativa inicial para ponderar com mais precisão as leituras de mapeamento para vários locais no genoma. ` --library-type fr-unstranded ` Reads da extremidade mais à esquerda do fragmento (em coordenadas de transcrição) mapeia para a cadeia de transcrição e a extremidade mais à direita mapeia para a cadeia oposta. `--m/ --frag-len-mean 300 ` Esse é o comprimento esperado do fragmento. (O padrão é 200bp.) `--s/ --frag-len-std-dev 50 ` O desvio padrão para a distribuição nos comprimentos dos fragmentos. (O padrão é 80bp.) ` --total-hits-norm ` Com esta opção, os cufflinks contam todos os fragmentos, incluindo aqueles não compatíveis com nenhuma transcrição de referência, em relação ao número de ocorrências mapeadas usadas no denominador FPKM. ` --max-frag-multihits 20 ` Este comando permite ignorar fragmentos que são mapeados para o genoma mais de um número especificado de vezes. ` --F/ --min-isoform-fraction 0.20 ` Depois de calcular a abundância de isoformas para um gene, os cufflinks filtram as transcrições que elas acreditam serem de abundância muito baixa, porque as isoformas expressas em níveis extremamente baixos geralmente não podem ser montadas com segurança e podem até ser artefatos de precursores de transcrições processadas incompletamente unidos. Esse parâmetro também é usado para filtrar íntrons que possuem muito menos alinhamentos emendados que os suportam. [O padrão é 0,1 ou 10% da isoforma mais abundante (a principal isoforma) do gene.] ` --I/ --max-intron-length 10000 ` O comprimento máximo do íntron. Os botões de punho não reportarão transcrições com íntrons maiores que isso e ignorarão os alinhamentos de SAM com as operações REF_SKIP CIGAR por mais tempo. (O padrão é 300.000.) `--min-intron-length 100 ` Tamanho mínimo de íntron permitido no genoma. (O padrão é 50 pb.) ` --overhang-tolerance 4 ` O número de bp permitido exceder a extremidade 3 'de uma transcrição de referência ao determinar se uma transcrição montada deve ser mesclada com ela (ou seja, a transcrição montada não é nova). (O padrão é 600 bp.) ` --max-bundle-frags 999999 ` Define o número máximo de fragmentos que um locus pode ter antes de ser ignorado. Os locais ignorados estão listados em skipped.gtf. (Padrão: 1000000) ` --max-multiread-fraction 0.45 ` A fração de leituras de suporte de um transfrag que pode ser multiplicada por mapeamento para o genoma. Uma transcrição composta por mais do que essa fração não será relatada pelo montador. [Padrão: 0,75 (75% de múltiplas linhas ou mais é suprimido).] `--overlap-radius 10 ` Transfrags separados por menos que essa distância são mesclados e a lacuna é preenchida. (Padrão: 50, em bp). `--3-overhang-tolerance 300 ` O número de bp permitido exceder a extremidade 3 'de uma transcrição de referência ao determinar se uma transcrição montada deve ser mesclada com ela (ou seja, a transcrição montada não é nova). (O padrão é 600 bp.) ::: :::info ### 1. Cuffcompare Depois de montar um transcriptoma de uma ou mais amostras, você provavelmente desejará comparar seu assembly com transcrições conhecidas. Mesmo se não houver transcriptoma de "referência" para o organismo que você está estudando, convém comparar os transcriptomas montados em diferentes bibliotecas de RNA-Seq. O Cuffcompare ajuda a realizar essas comparações e avaliar a qualidade da sua montagem. ::: :::warning ### Parâmetros utilizados com o Cuffcompare `-r ` An optional “reference” annotation GFF file. Cada amostra é comparada com esse arquivo e as isoformas da amostra são marcadas como sobrepostas, correspondentes ou novas, quando apropriado. ` -s ` Faz com que o cuffcompare procure arquivos fasta com as seqüências genômicas subjacentes (um arquivo por contig) com as quais suas leituras foram alinhadas para algumas funções de classificação opcionais. `-o ` Todos os arquivos de saída criados pelo Cuffcompare terão esse prefixo (por exemplo, .loci, .tracking etc.). ::: :::info ### 2. Cuffmerge Quando você possui várias bibliotecas RNA-Seq e reuniu transcriptomas de cada uma delas, recomendamos que você mescla esses assemblies em um transcriptoma mestre. Esta etapa é necessária para uma análise de expressão diferencial das novas transcrições que você montou. O Cuffmerge executa esta etapa de mesclagem. ::: :::warning ### Parâmetros utilizados com o Cuffmerge ` --o ` Escreva as estatísticas de resumo no arquivo de saída de texto <outprefix> . `--g/ --ref-gtf ` Uma anotação opcional de "referência" GTF. Os conjuntos de entrada são mesclados com o GTF de referência e incluídos na saída final. `--s/ --ref-sequence ` Este argumento deve apontar para as sequências de DNA genômico para a referência. Se um diretório, ele deve conter um arquivo fasta por contig. Se um arquivo multifasta, todos os contigs devem estar presentes. O script de mesclagem passará essa opção para comparar as algemas, que usarão as sequências para auxiliar na classificação de transfrags e na exclusão de artefatos (por exemplo, repetições). Por exemplo, as transcrições de abotoaduras que consistem principalmente em letras minúsculas são classificadas como repetições. ``--F/ --min-isoform-fraction 0.20 `` Depois de calcular a abundância de isoformas para um gene, as Abotoaduras filtram as transcrições que elas acreditam serem de abundância muito baixa, porque as isoformas expressas em níveis extremamente baixos geralmente não podem ser montadas com segurança e podem até ser artefatos de precursores de transcrições processadas incompletamente unidos. Esse parâmetro também é usado para filtrar íntrons que possuem muito menos alinhamentos emendados que os suportam. O padrão é 0,1 ou 10% da isoforma mais abundante (a principal isoforma) do gene. `--p/ --num-threads ` Use esse número de threads para alinhar as leituras. O padrão é 1. ::: :::info ### 3. Cuffquant A quantificação da expressão de genes e transcritos em amostras de RNA-Seq pode ser computacionalmente cara. O Cuffquant permite calcular os perfis de expressão de gene e transcrição e salvar esses perfis em arquivos que você pode analisar posteriormente com Cuffdiff ou Cuffnorm. Isso pode ajudá-lo a distribuir sua carga computacional em um cluster e é recomendado para análises que envolvem mais de um punhado de bibliotecas. ::: :::warning ### Parâmetros utilizados com o Cuffquant `--o/ --output-dir ` Define o nome do diretório no qual o Cuffquant gravará toda a sua saída. O padrão é "./". `--b/ --frag-bias-correct ` Fornece ao Cuffquant um arquivo multifasta por meio dessa opção instrui-o a executar o algoritmo de detecção e correção de viés, o que pode melhorar significativamente a precisão das estimativas de abundância de transcrição. `--u/ --multi-read-correct` Diz ao Cuffquant para realizar um procedimento de estimativa inicial para ponderar com mais precisão as leituras de mapeamento para vários locais no genoma. `--p/ --num-threads ` Use esse número de threads para alinhar as leituras. O padrão é 1. `--library-type fr-unstranded` Reads da extremidade mais à esquerda do fragmento (em coordenadas de transcrição) mapeia para a cadeia de transcrição e a extremidade mais à direita mapeia para a cadeia oposta. `--m/ --frag-len-mean 300 ` Esse é o comprimento (médio) esperado do fragmento. O padrão é 200bp. `--s/ --frag-len-std-dev 50 ` O desvio padrão para a distribuição nos comprimentos dos fragmentos. O padrão é 80bp. `--max-bundle-frags 9999999` Define o número máximo de fragmentos que um locus pode ter antes de ser ignorado. Padrão: 1000000 `--max-frag-multihits 20` Este comando permite ignorar fragmentos que são mapeados para o genoma mais de um número especificado de vezes. ::: :::info ### 4. Cuffdiff Comparar os níveis de expressão de genes e transcritos em experimentos de RNA-Seq é um problema difícil. O Cuffdiff é uma ferramenta altamente precisa para realizar essas comparações e pode dizer não apenas quais genes são regulados para cima ou para baixo entre duas ou mais condições, mas também quais genes estão emendados diferencialmente ou estão passando por outros tipos de regulação no nível de isoformas. ::: :::warning ### Parâmetros utilizados com o Cuffdiff `--o/ --output-dir ` Define o nome do diretório no qual o Cuffdiff gravará toda a sua saída. O padrão é "./". `--L/ --labels ` Especifica um rótulo para cada amostra, que será incluído em vários arquivos de saída produzidos pelo Cuffdiff. `--b/ --frag-bias-correct ` Fornece cufflinks com o arquivo multifasta para o qual suas leituras foram mapeadas por meio dessa opção instrui-o a executar nosso algoritmo de detecção e correção de viés, que pode melhorar significativamente a precisão das estimativas de abundância de transcrição. `--u/ --multi-read-correct ` Diz aos Cufflinks para fazer um procedimento de estimativa inicial para ponderar com mais precisão as leituras mapeadas para vários locais do genoma. `--p/ --num-threads` Use esse número de threads para alinhar as leituras. O padrão é 1. `--library-type fr-unstranded ` Reads da extremidade mais à esquerda do fragmento (em coordenadas de transcrição) mapeia para a cadeia de transcrição e a extremidade mais à direita mapeia para a cadeia oposta. `--m/ --frag-len-mean 300 ` Esse é o comprimento (médio) esperado do fragmento. O padrão é 200bp. `--frag-len-std-dev 50` O desvio padrão para a distribuição nos comprimentos dos fragmentos. O padrão é 80bp. `--max-bundle-frags 9999999 ` Define o número máximo de fragmentos que um locus pode ter antes de ser ignorado. Os locais ignorados são marcados com o status HIDATA. Padrão: 1000000 `--max-frag-multihits 20` Este comando permite ignorar fragmentos que são mapeados para o genoma mais de um número especificado de vezes. `--total-hits-norm ` Com esta opção, os cufflinks contam todos os fragmentos, incluindo aqueles não compatíveis com nenhuma transcrição de referência, em relação ao número de fragmentos mapeados usados no denominador FPKM. `--min-reps-for-js-test 2` Cuffdiff não testará genes para regulação diferencial, a menos que as condições em questão tenham pelo menos tantas repetições. Padrão: 3. `--library-norm-method geometric ` Os FPKMs e as contagens de fragmentos são dimensionados pela mediana das médias geométricas das contagens de fragmentos em todas as bibliotecas, conforme descrito em Anders e Huber (Genome Biology, 2010). `--dispersion-method per-condition ` Cada condição replicada recebe seu próprio modelo. Disponível apenas quando todas as condições tiverem réplicas. `--c/ --min-alignment-count 10 ` O número mínimo de alinhamentos em um locus necessário para realizar testes de significância nas alterações nesse locus observadas entre as amostras. Se nenhum teste for realizado, as alterações no locus são consideradas não significativas e as alterações observadas no locus não contribuem para a correção de vários testes. O padrão é 10 alinhamentos de fragmentos. :::: :::info ### 5. Cuffnorm Às vezes, tudo o que você quer fazer é normalizar os níveis de expressão de um conjunto de bibliotecas de RNA-Seq para que elas estejam na mesma escala, facilitando análises posteriores, como agrupamentos. Os níveis de expressão relatados pelos botões de punho nas unidades FPKM são geralmente comparáveis entre as amostras, mas em determinadas situações, a aplicação de um nível extra de normalização pode remover fontes de viés nos dados. O Cuffnorm normaliza um conjunto de amostras para estar na escala mais semelhante possível, o que pode melhorar os resultados obtidos com outras ferramentas a jusante. ::: :::warning ### Parâmetros utilizados com o Cuffnorm `--o/ --output-dir ` Define o nome do diretório no qual o Cuffdiff gravará toda a sua saída. O padrão é "./". `--L/ --labels ` Especifica um rótulo para cada amostra, que será incluído em vários arquivos de saída produzidos pelo Cuffdiff. `--p/ --num-threads` Use esse número de threads para alinhar as leituras. O padrão é 1. `--library-type fr-unstranded ` Reads da extremidade mais à esquerda do fragmento (em coordenadas de transcrição) mapeia para a cadeia de transcrição e a extremidade mais à direita mapeia para a cadeia oposta. `--library-norm-method geometric ` Os FPKMs e as contagens de fragmentos são dimensionados pela mediana das médias geométricas das contagens de fragmentos em todas as bibliotecas, conforme descrito em Anders e Huber (Genome Biology, 2010). ` --output-format simple-table` Por padrão, o Cuffnorm relata os níveis de expressão nos arquivos de texto delimitados por tabulação “tabela simples”. O programa também relata informações sobre suas amostras e sobre os genes, transcrições, grupos TSS e grupos CDS como arquivos de texto delimitados por tabulação. ::: :::info ## StringTie O StringTie é um montador rápido e altamente eficiente de alinhamentos RNA-Seq em possíveis transcrições. Ele usa um novo algoritmo de fluxo de rede, bem como uma etapa de montagem opcional de novo para montar e quantificar transcrições completas que representam múltiplas variantes de emenda para cada locus genético. Sua entrada pode incluir não apenas alinhamentos de leituras curtas que também podem ser usadas por outros montadores de transcrição, mas também alinhamentos de sequências mais longas que foram montadas a partir dessas leituras. Para identificar genes expressos diferencialmente entre experimentos, a produção do StringTie pode ser processada por software especializado como Ballgown, Cuffdiff ou outros programas (DESeq2, edgeR, etc.). ::: :::warning ### Parâmetros utilizados no StringTie `--merge` É o modo de mesclagem de transcrição. Neste modo de mesclagem, o StringTie recebe como entrada uma lista de arquivos GTF / GFF e mescla / monta essas transcrições em um conjunto não redundante de transcrições. Esse modo é usado na análise diferencial para gerar um conjunto global e unificado de transcritos (isoformas) em várias amostras de RNA-Seq. Quando a opção -G (anotação de referência) for fornecida, o StringTie montará os transfrags dos arquivos GTF de entrada com as transcrições de referência. ` --G` Use o arquivo de anotação de referência (no formato GTF ou GFF3) para orientar o processo de montagem. A saída incluirá transcrições de referência expressas, bem como quaisquer novas transcrições reunidas. ` --o ` Define o nome do arquivo GTF de saída em que StringTie gravará as transcrições montadas. Isso pode ser especificado como um caminho completo; nesse caso, os diretórios serão criados conforme necessário. Por padrão, o StringTie grava o GTF na saída padrão. ` --p ` Especifique o número de threads de processamento (CPUs) a serem usados para montagem de transcrição. O padrão é 1. `--f 0.20 ` Define a abundância isoforma mínima dos transcritos previstos como uma fração do transcrito mais abundante reunida em um determinado local. Transcrições de menor abundância são frequentemente artefatos de precursores de transcrições processadas incompletamente unidos. Padrão: 0.1 `--a 10 ` As junções que não têm leituras emendadas que se alinham entre elas com pelo menos essa quantidade de bases nos dois lados são filtradas. Padrão: 10 `--j 3 ` Deve haver pelo menos tantas leituras emendadas que se alinham em uma junção (ou seja, cobertura da junção). Esse número pode ser fracionário, pois algumas leituras se alinham em mais de um lugar. Uma leitura alinhada em n lugares contribuirá 1 / n para a cobertura da junção. Padrão: 1 `--c 2` Define a cobertura mínima de leitura permitida para as transcrições previstas. Uma transcrição com uma cobertura menor que esse valor não é mostrada na saída. Padrão: 2.5 `--g 10` Valor mínimo de separação do espaço entre pontos. As leituras mapeadas mais próximas que essa distância são mescladas no mesmo pacote de processamento. Padrão: 50 (bp) `--M 0.45 ` Define a fração máxima de leituras mapeadas em vários locais que podem estar presentes em um determinado local. Padrão: 0,95. `--A` A abundância de genes será relatada (formato delimitado por tabulações) no arquivo de saída com o nome fornecido. ::: ## 5. Montagem de novo de transcriptomas Os dados de entrada para este programa serão os mesmos utilizados com o protocolo baseado em referência usando Cufflinks (Tuxedo). Esses arquivos processados estão no diretório ./output/processed/prinseq. Para preparar os dados para montagem com Trinity além de construir a linha de comando, foi desenvolvido o script rnaseq-denovo.sh para realizar a montagem de novo. ## rnaseq-denovo.sh Link: [Script rnaseq-denovo.sh](https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-denovo.sh) ``` wget https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-denovo.sh ``` :::danger ## Descrição dos softwares utilizados para a montagem dos scripts para a análise do RNA-Seq sem referência - "Montagem de novo" ::: :::info ## Trinity O Trinity monta sequências de transcrição a partir de dados de Illumina RNA-Seq. O Trinity representa um novo método para a reconstrução eficiente e robusta de novo de transcriptomas a partir de dados de RNA-seq. O Trinity combina três módulos de software independentes: Inchworm, Chrysalis e Butterfly, aplicados sequencialmente para processar grandes volumes de leituras de RNA-seq. O Trinity particiona os dados da sequência em muitos gráficos individuais de Bruijn, cada um representando a complexidade da transcrição em um determinado gene ou locus, e depois processa cada gráfico independentemente para extrair isoformas de emenda de comprimento total e separar os transcritos derivados de genes paralógicos. Resumidamente, o processo funciona assim: ### Inchworm Reúne os dados de RNA-seq em seqüências únicas de transcrições, geralmente gerando transcrições completas para uma isoforma dominante, mas depois re:porta apenas as partes únicas de transcrições com splicing alternativo. ### Chrysalis Agrupa os contágios do Inchworm em agrupamentos e constrói gráficos completos de Bruijn para cada agrupamento. Cada cluster representa a complexidade transcriptonal completa de um determinado gene (ou conjuntos de genes que compartilham seqüências em comum). Crisálida, em seguida, particiona o conjunto completo de leitura entre esses gráficos disjuntos. ### Butterfly Processa os gráficos individuais em paralelo, traçando os caminhos que as leituras e os pares de leituras percorrem no gráfico, relatando, em última análise, transcrições completas para isoformas alternadas emendadas e separando as transcrições que correspondem a genes paralógicos. ::: :::warning ### Parâmetros utilizados no Trinity `--KMER_SIZE 27 ` Tamanho do Kmer `--output ` Nome do diretório de saída `--seqType fq ` Os tipos de reads, podem ser: fa ou fq. `--max_memory ` Memória máxima sugerida para uso pelo Trinity, onde a limitação pode ser ativada. O padrão é «--max_memory 10G» `--CPU` Número de CPUs a serem usadas, padrão: 2 `--min_per_id_same_path 95 ` `--max_diffs_same_path 5 ` `--path_reinforcement_distance 5` `--group_pairs_distance 500 ` Se você tiver vários tamanhos de fragmento de biblioteca de extremidade pareada, defina '--group_pairs_distance' de acordo com a biblioteca de inserção maior. Os pares que excederem essa distância serão tratados como se não fossem emparelhados pelo processo Butterfly. `--min_glue 5 ` `--min_contig_length 600` Comprimento mínimo contig montado para relatar, padrão = 200. `--min_kmer_cov 3 ` Tamanho mínimo do Kmer. `--genome_guided_bam ` Modo guiado por genoma, forneça o caminho para o arquivo bam classificado por coordenadas. ``--genome_guided_max_intron 10000 `` ## "Estimating abundances" `--transcripts` Arquivo de transcrição fasta. `--est_method` Método de estimativa da abundância. `--salmon_add_opts ` salmon opções: --salmon_idx_type <string> quasi|fmd (defalt: quasi) --salmon_add_opts <string> default: `--samples_file ` Arquivo de texto delimitado por tabulação indicando relacionamentos de replicação biológica. `--gene_trans_map` O arquivo de mapeamento de gene para transcrição. (Certifique-se de usar os parâmetros --gene_trans_map ou --trinity_mode para obter uma matriz de contagem de genes, além da matriz de contagem de isoformas.) `--prep_reference ` Referência de preparação (cria o índice de destino). `--thread_count ` Número de threads a serem usados (padrão = 4). `--seqType fq ` `--output_dir` Gravar todos os arquivos no diretório de saída. ## "Constructing abundance matrix" `--est_method salmon` `--gene_trans_map` Arquivo contendo identificadores 'gene (tab) transcript' por linha. `--name_sample_by_basedir` Nome da coluna de amostra pelo nome do diretório em vez do nome do arquivo `--cross_sample_norm none` TMM|UpperQuartile|none (padrão: TMM) `--quant_files` Arquivo contendo uma lista de todos os arquivos de destino. `--out_prefix` Padrão: valor para --est_method. ## "Calculating Differentially Expressed Genes" ### run-DESeq2.R O pacote DESeq2 fornece métodos para testar a expressão diferencial usando modelos lineares generalizados binomiais negativos; as estimativas de alterações de dispersão e dobras logarítmicas incorporam distribuições anteriores orientadas por dados. `-in` `--groups` `--out` ::: ## 5.1 Para avaliação das montagens #### A) Indexação do transcriptoma utilizado para a simulação: ``` makeblastdb -title "RefTrans" \ -dbtype nucl \ -out RefTrans \ -in transcriptoma.fa \ > makebleastdb.log.out.txt \ 2> makeblastdb.log.err.txt ``` Para conferir se esta indexação foi bem sucedida, neste caso, aparecerão os seguintes arquivos: - RefTrans.nhr - RefTrans.nin - RefTrans.nsq :::info ## Makeblastdb O makeblastdb é uma ferramenta de linha de comando para criar bancos de dados BLAST. Faz parte do novo pacote blast + do NCBI. O aplicativo produz bancos de dados BLAST a partir de arquivos FASTA. É possível usar linhas de definição FASTA completamente não estruturadas (ou mesmo em branco), mas este não é o procedimento recomendado. A designação de um identificador exclusivo para cada sequência no banco de dados permite recuperar a sequência por identificador e permite associar todas as sequências a um nó taxonômico (por meio do taxid da sequência). O identificador exclusivo pode ser uma sequência simples ou pode ser a acessão real da sequência se a sequência vier de um banco de dados público (por exemplo, GenBank). ::: :::warning ### Parâmetros utilizados no Makeblastdb `-title` Título do banco de dados BLAST. Se não estiver definido, o nome do arquivo de entrada será usado. ``-dbtype `` Tipo de molécula de entrada em um banco de dados BLAST, podem ser nucl ou prot. No script: ``-dbtype nucl`` => a molécula escolhida foi do tipo nucleotídeo `-out ` Nome do banco de dados BLAST a ser criado. O nome do arquivo de entrada é usado se nenhum for fornecido. Este campo é obrigatório se a entrada consistir em vários arquivos. No script, foi utilizado: `-out RefTrans` `-in` Arquivo de entrada / nome do banco de dados. No script, escolhemos: `-in transcriptoma.fa`. ::: #### B) Para executar a comparação vamos utilizar o blastn: Primeiro, para a montagem com Trinity de novo: ``` blastn -max_hsps 1 \ -max_target_seqs 1 \ -num_threads 8 \ -query ./output/trinity_assembled/Trinity.fasta \ -task blastn \ -db ./RefTrans \ -out ./Trinity_x_RefTrans.txt \ -evalue 1e-5 \ -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \ > ./Trinity_x_RefTrans.log.out.txt \ 2> ./Trinity_x_RefTrans.log.err.txt ``` #### C) Visualizar os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos: ``` cut -f 1,2 ./Trinity_x_RefTrans.txt ``` #### D) Contar quantas ocorrências de sequências da base de dados (subject) foram encontradas: `cut -f 2 ./Trinity_x_RefTrans.txt | sort | uniq -c` #### E) Depois, para a montagem com Trinity Genome Guided: ``` blastn -max_hsps 1 \ -max_target_seqs 1 \ -num_threads 8 \ -query ./output/trinity_GG_assembled/Trinity-GG.fasta \ -task blastn \ -db ./RefTrans \ -out ./Trinity-GG_x_RefTrans.txt \ -evalue 1e-5 \ -outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \ > ./Trinity-GG_x_RefTrans.log.out.txt \ 2> ./Trinity-GG_x_RefTrans.log.err.txt ``` #### F) Visualizar os nomes das sequências query (qseqid) e subject (sseqid) HSPs de todos HITs obtidos: ``` cut -f 1,2 ./Trinity-GG_x_RefTrans.txt ``` #### G) Contar quantas ocorrências de sequências da base de dados (subject) foram encontradas: ``` cut -f 2 ./Trinity-GG_x_RefTrans.txt | sort | uniq -c ``` :::info ## Blastn BLAST nucleotideo (BLASTn) O programa blastn é usado quanto temos uma sequência de DNA (query) e queremos buscar a sequência de DNA (subject) do banco de dados similar a nossa sequência, são chamados de alinhamentos nucleotídeo-nucleotídeo. As comparações nucleotídeo-nucleotídeo tem duas variantes dentro da família BLAST, o programa blastn e o programa MegaBLAST. O programa blastn é usado para fazer buscas e alinhamentos de segmentos de DNA mais gerais, usando regiões de genomas onde encontramos misturas de regiões codificantes e não codificantes, comparaçoes de rRNA, tRNA e mRNA. O MegaBLAST faz a comparação 10 vezes mais rápida que o blastn, mas ele busca sequências que são extremamente similares, permitindo mapeamento de transcritos no genoma de 3 bilhões de bases em segundos, e útil para processamento de um grande conjunto de sequências. ::: :::warning ### Parâmetros utilizados no Blastn `-max_hsps 1 ` Número máximo de HSPs (alinhamentos) a serem mantidos por um único par de assuntos de consulta. Os HSPs mostrados serão os melhores conforme julgado pelo valor esperado. Esse número deve ser um número inteiro que seja um ou mais. Se essa opção não estiver definida, o BLAST mostrará todos os HSPs que atendem aos critérios de valor esperado. `-max_target_seqs 1 ` Número de sequências alinhadas a serem mantidas. Use com formatos de relatório que não tenham linhas de definição e seções de alinhamento separadas, como tabular (todos os outfmt> 4). Não é compatível com num_descriptions ou num_alignments. Os laços são quebrados por ordem de sequências no banco de dados. `-num_threads 8 ` Número de threads (CPUs) a serem usados na pesquisa por blast. `-query ` Nome do arquivo de consulta. -task blastn O aplicativo blastn pesquisa uma consulta de nucleotídeo em relação a sequências sujeitas a nucleotídeos ou um banco de dados de nucleotídeos. Quatro tarefas diferentes são suportadas: 1.) megablast, para seqüências muito semelhantes, 2.) dc-megablast, normalmente usado para comparações entre espécies 3.) blastn, o tradicional programa usado para comparações interespécies 4.) blastn-short, otimizado para seqüências com menos de 30 nucleotídeos. `-db ` `-out ` `-evalue 1e-5 ` `-outfmt ` ::: ## 5.2. Para estimar a abundância e encontrar os genes/isoformas diferencialmente expressos #### A) Para executar esse protocolo de análise de expressão gênica diferencial de uma forma automatizada, foi criado o script rnaseq-trinity-abundance.sh. `cp /usr/local/bioinfo/bioaat/rnaseq-trinity-abundance.sh .` `./rnaseq-trinity-abundance.sh output/renamed/ output/trinity_assembled/ output/` ## rnaseq-trinity-abundance.sh. Link: [Script rnaseq-trinity-abundance.sh.](https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-trinity-abundance.sh) ``` wget https://github.com/dbonaldi/RelatorioBio/blob/master/rnaseq-trinity-abundance.sh ``` #### B) Para visualizar o resultado da anális com DESeq2 de modo alinhado: ``` cat output/abundance/DEG/SAMPLEA-SAMPLEB.txt | column ``` #### C) Para ordenar essas colunas: `cat ./output/abundance/DEG/SAMPLEA-SAMPLEB.txt | column -s$'\t' -t | less -S` | X | SAMPLEA2 | SAMPLEB1 | SAMPLEB2 | SAMPLEA1 | logFC | PValue | FDR | |----|----|----|----|----|----|----|----| | TRINITY_DN1_c0_g1| 1863.96003271285| 551.175741513274 | 649.090298609964| 1122.72072683247| 1.31404557077059 | 0.000113022839817923 | 0.000226045679635846| | TRINITY_DN1_c0_g1| 2016.10574896126| 7747.01265987757 |6398.44650323994| 2101.01637257929| -1.78061446741346 | 5.79241866606607e-10| 2.31696746642643e-09| | TRINITY_DN3_c0_g1| 1454.13257661881| 1113.20252091598 | 1317.75093735092| 1005.57507268683| 0.0163544028134741 | 0.959223154295421 | 0.959223154295421| | TRINITY_DN4_c0_g1| 6437.93238256796 | 6415.46184200393 | 11385.5387486897| 8712.82586688147 | -0.232385041676405 | 0.50124565561255 | 0.668327540816734| # 6. Acessando o RStudio para algumas análises dos resultados Na interface do RStudio inserir o script para selecionar os genes diferencialmente expressos: Link: [Script RStudio](https://github.com/dbonaldi/RelatorioBio/blob/master/RStudio) ``` wget https://github.com/dbonaldi/RelatorioBio/blob/master/RStudio ```` ::: info ## heatmap.2 Mapas de calor são comumente usados para visualizar resultados de RNA-Seq. Eles são úteis para visualizar a expressão de genes nas amostras. O mapa de calor também pode ser combinado com métodos de agrupamento que agrupam genes e / ou amostras com base na similaridade do seu padrão de expressão gênica. Isso pode ser útil para identificar genes comumente regulamentados ou assinaturas biológicas associadas a uma condição específica. Nos mapas de calor, os dados são exibidos em uma grade onde cada linha representa um gene e cada coluna representa uma amostra, sendo que a cor e a intensidade das caixas são usadas para representar alterações (não valores absolutos) da expressão do gene. ::: ::: success ### Resultado do heatmap.2 ::: ![](https://i.imgur.com/5zNJ6qB.png) O heatmap.2 nos permite a visualização dos genes TRINITY_DN1_cO_g1 e TRINITY_DN2_cO_g1 que estão com expressão up regulation.