(tem mais coisas antes. O que está descrito aqui é da aula do dia 01/11/23)
limpar o genome.fa, tirando os títulos dos >.
Seguir as instruções para instalação do e-direct.
Caso precise atualizar o conjunto de programas do e-direct, remova primeiro a pasta correspontende
```
```
Sequenciamento
é bem comum encontrar contaminação por adaptador na porção 3' de adaptador em R1 e R2, isso acontece porque não é sequenciado por completo o inserto. Isso é bem comum, um erro da máquina que tem que ser considerado.
avaliação do RNA -> eletroforese das bandas do 28S/18S (razão de 2:1) em eucarioto e 23S/16S. O bioanalyzer faz eletroforese capilar, que é bem mais precisa.Os picos de sinal são correspondentes às bandas do gel de agarose.
rRNA = problema para sequenciamento, pela quantidade.
Tem como enriquecer a amostra com RNAs mais raros usando a ideia de melting. Os mais abunantes vão encontrar seu parceiro mais rápido, então coloca uma enz de restrição para cortar os RNAs mais abundantes e os mais raros são enriquecidos.
POr questões da máquina e processo, a porção 3' tem qualidade menor e R2 tenha sua porção 3' ainda mais prejudicada. Mas isso pode ser arrumado com R1 e R2 sendo sobreposto até certo nivel. Quando a qualidade da base é igual e as bases são diferentes/divergentes é melhor deixar um N (indetermidado) a colocar uma base que não seja a real.
RiN - RNA integrity. Essa classificação é fruto de um algoritmo de aprendizado de máquinas que mostra o que é íntegro do que não é íntegro. Varia de 0 a 10, sendo que 7 é minimamente bom para o sequenciamento.
É feita uma qPCR para avaliar a quantidade de moléculas antes de ser sequenciada.
fasta:
@nome da máquina + detalahmentos do sequenciamento
bases nitrogenadas
+(informações sobre)
caracteres de qualidade (referidos na tabela ascii), que varia de 33 ou 64 de codificação. Hoje se usa só o 33, programas mais antigos era 64.
em que saber qual a versão da illumina phred que você está usando!! o fastq precisa de menos bytes para passar a mesma informação do fasta com o (qual) - índice de qualidade.
metadados - dados que explicam os dados. isso é bastante comum no depósito dos dados no ncbi ou outras bases de dados.
Onde ficam as sequências de RNA-seq => SRA no NCBI
PROCESSAMENTO E ANÁLISE DE DADOS
alinhamento = princípio fundamental da bioinformática
matrizes de substituíção de nucleotídeos: I)jukes-Kantor (mesma probabilidade de troca de nucleotpides, sem levar em consideração a natureza da base) ou II) Kimura (leva em consideração a natureza de transversão, da natureza das bases)
PAM = global e Blossum = local
Alinhamentos PAM com valores menores = menor divergência entre elas, mais próximo evolutivamente.
e-value - número de chances de encontrar um match igual ou melhor do que eu encontrei ao acaso, para o tamanho do banco de dados naquele momento.
MATRIX DE PONTOS - identifica o match de bases idênticas, análise qualitativa de estrutura. A diagonal mostra as regiões similares/idênticas. Indica assinaturas de inversão e repetições
ALINHAMENTOS CONSIDERANDO AS NOVAS GERAÇÕES DE SEQUENCIAMENTO
o Bowtie usa algoritmos de Árvores de prefixos e sufixos que busca a maior profundidade. Ele segmenta todas as letras e dá um valor para ela, por exemplo de matches e mismatches, e "soma" esses valores. No final, ele verifica quem foi mais profundo. Ele também inicia de trás pra frente na sequência.
Transformação de Burrows-Wheeler - ele maximiza a compressão de dados, por meio de ordernações diferenciais.
O bowtie necessita de um índice para trabalhar
Prática de alinhamento:
Fizemos uma nova pasta index pra trabalhar nela. Detalhe importante: criamos um link "de referência" acessar o genome.fa da outra pasta
```
mkdir index
cd index
ln -s ../ref/genome.fa
bowtie2-build -f genome.fa genome
```
o inspect volta o que era antes, sem o index (o less é só pra aparecer a última página)
```
bowtie2-inspect genome | less
```
Depois para fazer o alinhamento de uma sequência aleatória no arquivo genome e exportar para o arquivo test.sam
```
bowtie2 -c "GTTCCCAGCGAGAGCGTTACGCATGGATCTCGCATGAGCACGCCTCAG" -x genome -S test.sam
```
Para alinhar sequências pair-ended, isso com sequências feitas à mão. Para arquivos fasta é o próximo passo.
```
bowtie2 -c -1 "GTGTAACTACGAAACATTTGCCAGCATTTTATTATTCATCTCCTGATTTCTATTTTGGTAGAGACCAGGA" -2 "GCAGGGACTGCCAACCGTCGAAATCTGGTACTAAGCTTCAAAGTTTTTGAAAGCCTTTGACTTCTATCAA" -x genome -S test.sam --maxins 1000
```
> 1 reads; of these:
1 (100.00%) were paired; of these:
0 (0.00%) aligned concordantly 0 times
1 (100.00%) aligned concordantly exactly 1 time
0 (0.00%) aligned concordantly >1 times
----
0 pairs aligned concordantly 0 times; of these:
0 (0.00%) aligned discordantly 1 time
----
0 pairs aligned 0 times concordantly or discordantly; of these:
0 mates make up the pairs; of these:
0 (0.00%) aligned 0 times
0 (0.00%) aligned exactly 1 time
0 (0.00%) aligned >1 times
100.00% overall alignment rate
Alinhando arquivos inteiros:
> bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index$ bowtie2 -1 ../raw/SAMPLEC1_R1.fastq -2 ../raw/SAMPLEC1_R2.fastq -x genome -S test.sam --maxins 600
974 reads; of these:
974 (100.00%) were paired; of these:
878 (90.14%) aligned concordantly 0 times
96 (9.86%) aligned concordantly exactly 1 time
0 (0.00%) aligned concordantly >1 times
----
878 pairs aligned concordantly 0 times; of these:
327 (37.24%) aligned discordantly 1 time
----
551 pairs aligned 0 times concordantly or discordantly; of these:
1102 mates make up the pairs; of these:
681 (61.80%) aligned 0 times
421 (38.20%) aligned exactly 1 time
0 (0.00%) aligned >1 times
65.04% overall alignment rate
**O bowtie não entende splicing**, ele acha que é genoma com genoma por isso que deu só 96 reads que alinham 1 vez só. não funciona direito para eucariotos. Se for alinhamento rna-rna é de boa, é uma boa ferramenta. Ele não consegue alinhar de forma descontínua, por exemplo uma read que pega a região de união de exon. Ele não é capaz de "dividir" a fração de cada éxon. Daí, neste caso elas são classificadas como reads inicialmente não mapeadas. Ele também é capaz de predizer regiões de introns e "faz" o splicing e tenta alinhar essas inicialmente não mapeadas. Se der certo, boa, ele alinha. Caso não, a read é perdida.
Fazendo com o bwa
```
bwa index genome.fa
```
> [bwa_index] Pack FASTA... 0.37 sec
[bwa_index] Construct BWT for the packed sequence...
[BWTIncCreate] textLength=86550302, availableWord=18089684
[BWTIncConstructFromPacked] 10 iterations done. 29839102 characters processed.
[BWTIncConstructFromPacked] 20 iterations done. 55123726 characters processed.
[BWTIncConstructFromPacked] 30 iterations done. 77592734 characters processed.
[bwt_gen] Finished constructing BWT in 35 iterations.
[bwa_index] 24.96 seconds elapse.
[bwa_index] Update BWT... 0.22 sec
[bwa_index] Pack forward-only FASTA... 0.23 sec
[bwa_index] Construct SA from BWT and Occ... 10.10 sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa index genome.fa
[main] Real time: 35.982 sec; CPU: 35.872 sec
vendo com o bwa agora
```
bwa mem genome.fa ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq > bwa_sampleC1.sam
```
o **samtools** é um programa que te permite trabalhar com arquivos sam. Para organizar o arquivo sam por nome para ser comparável:
```
samtools sort -n -o bwa_SAMPLEC1_sorted.sam bwa_sampleC1.sam
samtools sort -n -o test_sorted.sam test.sam
```
Para ver os resultados
*o SAM_nameSorted_to_uniq_counts_stats.pl já está no path, então é só rodar
```
SAM_nameSorted_to_uniq_count_stats.pl bwa_SAMPLEC1_sorted.sam > bwa_SAMPLEC1_sorted_results.txt
SAM_nameSorted_to_uniq_count_stats.pl test_sorted.sam > test_sorted_results.txt
```
Resultado do bwa_SAMPLEC1_sorted_results.txt
>496 aligned fragments; of these:
496 were paired; of these:
1 aligned concordantly 0 times
197 aligned concordantly exactly 1 time
298 aligned concordantly >1 times
----
1 pairs aligned concordantly 0 times; of these:
0 aligned as improper pairs
1 pairs had only one fragment end align to one or more contigs; of these:
1 fragments had only the left /1 read aligned; of these:
1 left reads mapped uniquely
0 left reads mapped >1 times
0 fragments had only the right /2 read aligned; of these:
0 right reads mapped uniquely
0 right reads mapped >1 times
Overall, 99.80% of aligned fragments aligned as proper pairs
Resultado do test_sorted_results.txt
>461 aligned fragments; of these:
461 were paired; of these:
371 aligned concordantly 0 times
84 aligned concordantly exactly 1 time
6 aligned concordantly >1 times
----
371 pairs aligned concordantly 0 times; of these:
234 aligned as improper pairs
137 pairs had only one fragment end align to one or more contigs; of these:
73 fragments had only the left /1 read aligned; of these:
61 left reads mapped uniquely
12 left reads mapped >1 times
64 fragments had only the right /2 read aligned; of these:
54 right reads mapped uniquely
10 right reads mapped >1 times
Overall, 19.52% of aligned fragments aligned as proper pairs
no bowtie, não colocamos parametros de sensibilidade. Isso pode ter contribuido para que o alinhamento tenha dado pior
Fazendo com --very-sensitive não mudou nada
> 461 aligned fragments; of these:
461 were paired; of these:
371 aligned concordantly 0 times
84 aligned concordantly exactly 1 time
6 aligned concordantly >1 times
----
371 pairs aligned concordantly 0 times; of these:
237 aligned as improper pairs
134 pairs had only one fragment end align to one or more contigs; of these:
71 fragments had only the left /1 read aligned; of these:
59 left reads mapped uniquely
12 left reads mapped >1 times
63 fragments had only the right /2 read aligned; of these:
53 right reads mapped uniquely
10 right reads mapped >1 times
Overall, 19.52% of aligned fragments aligned as proper pairs
ALINHAMENTO DE SEQUÊNCIAS DE TRANSCRIPTOMAS (RNA-seq)
TopHat
Tem como avaliar onde estão as uniões/junções de éxons. È mais díficil, porque tem que ter evidências como parte da sequência dos dois lados dos exons e com múltiplas cópias para aumentar a confiabilidade de que ali é mesmo uma união de exon.
* ele está no path, agora
o tophat precisa do python 2, e agora ele está no 3. Então tem que entrar nesse ambiente
cd
```
bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index$ source /usr/local/bioinfo/anaconda3/etc/profile.d/conda.sh
bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index$ conda activate python2
(python2) bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index$
```
o (python2) mostra que deu certo!!
o tophat usa o mesmo indice do bowtie2
para rodar o tophat
```
tophat2 genome ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq
```
*para procariotos tem que colocar --no-novel-junc e --max-intron-length 5000000 (valores absurdos!! justamente para que não haja problemas e que não exista na realidade) -> tudo iso para garantir que ele não busque junções de introns.
ele cria um novo diretório tophat_out com arquivo .bam. Permite identificação de inserção e deleção nos arquivos deletions.bed ou insertions.bed.
junctions.bed mostra onde ele prediz que tenha união de éxons.
Ele só mostra o que alinhou no formato bam. Como é um formato 'ruim', vamos transformar em sam
```
samtools view -h accepted_hits.bam > accepted_hits.sam
```
> (python2) bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index/tophat_out$ samtools view -h accepted_hits.bam > accepted_hits.sam
(python2) bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/index/tophat_out$ ls
accepted_hits.bam align_summary.txt insertions.bed logs unmapped.bam
accepted_hits.sam deletions.bed junctions.bed prep_reads.info
```
more accepted_hits.sam
```
> @HD VN:1.0 SO:coordinate
@SQ SN:NC_023750.1 LN:43275151
@PG ID:TopHat VN:2.1.1 CL:/usr/local/bioinfo/tophat-2.1.1.Linux_x86_64/tophat ge
nome ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq
@PG ID:samtools PN:samtools PP:TopHat VN:1.16.1 CL:samtools view -h accep
ted_hits.bam
Frag_1 161 NC_023750.1 41452501 50 151M = 41454868 2929 G
CTTCTTCTTGTTAACAGATTCGTTCACAGAACTCCCCACCAACTTTCTAGATGCTGCCAAGAAGCCTGAGGTATCTGACTGGATGGTCAGAATCAGA
AGGAAGATTCACCAGAATCCAGAATTGGGTTATGAGGAATTTGAGACCAGTAA FECGGCCG=GFGGEF@<GGF/GFBGEGG?9FGF@FFCFGG@
G2EFGEGG9GAGFDEGCEA@AFGBG4GG1GFAFG'B=F4@FEDEFDD;/1/EC?EF>CCFD>FCD><ACB0FD=7B@$E@43B=&E:0<>D=2?/&2
A<:66252;?8,& AS:i:-2 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:75T75 YT:Z:UU NH:i:1
Frag_9 161 NC_023750.1 41452502 50 151M = 41452766 2439 C
TTCTTCTTGTTAACAGATTCGTTCACAGAACTCCCCACCAACTTTCTAGATGCTGCCAAGAAGCCTGAGGTATTTGACTGGATGGTCAGAATCAGAA
GGAAGATTCACCAGAATCCAGAATTGGGTTATGAGGAATTGGAGCCCAGTAAA >DFFFG&FGFEG@EE8GGGFGGGGGFGEGGGGGEFGGFG@?
GEEGEEGGG9F<FFEGG@FG<CA31GCAGFF4GFG(@<=GGBGG>FF<GCGF-FFC86G>FFBD8DE7.=+E5D@@@>:?>8FD7DDA6DA793=%/
%&=A(?3:=7..% AS:i:-4 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:138T3A8 YT:Z:UU NH:i:1
Frag_4 161 NC_023750.1 41452525 50 151M = 41455439 3065 T
CACAGAACTCCCCACCAACTTTCTAGATGCTGCCAAGAAGCCTGAGGTATTTGACTGGATGGTCAGAATCAGAAGGAAGATTCACCAGAATCCAGAA
TTGGGTTATGAGGAATTTGAGACCAGTAAAGTCATTAGAGAGAAATTGGATCA EFGGGGGGGFCGGEGGG;GGFD?FGDGGGCDGGFG?GGGFE
GG8GF@EGFG;GGFGGGGGF7GGGGGGEGGDFFCGFGG-EG0GFEGFEGGD4GFGFGAFDFC/ECFGDFGBD:;=GFFFF%DEE@?CBA?>BF?E0D
4F<D4EF1>;D70 AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:151 YT:Z:UU NH:i:1
Frag_6 161 NC_023750.1 41452633 50 151M = 41457270 4788 A
GGAATTTGAGACCAGTAAAGTCATTAGAGAGAAATTGGATCAACTGGGCATACCCTATAAACACCCTGTTGCAATCACTGGTGTAATTGGCTTCATT
GGGACTGGAAGATCGCCTTTTGTTGCTATAAGAGCTGATATGGATGCTCTTCC @GFGGFB?FG*GEGGGG<BGFG@FFGG>GGFGGDGGGGGGG
FDD>GGGG(EGG=GGCGFFBGFFG>GGGGGBADGGG@;@GGGB@GGF:FBF=GGACEGFGFDGGBEFFCG>-%:F@@FG?;0/F@F@BC&@A>EAE4
3DD<AE399A0=. AS:i:-2 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:49T101 YT:Z:UU NH:i:1
Frag_9 81 NC_023750.1 41452766 50 25M2024N126M = 41452502 -
2439 TGATATGGATGCTCTTCCCATGCAGGAAATGGTGGAGTGGGAGCACAAGAGTAAAGTACCTGGAAAGATGCATGCTTGTGGTCATGATG
CTCATGTGGCTATGCTTCTTGGTGCTGCGAAGATTCTTAAACAGCATGAAAAAGAGATACAA %2>D5?AAE8A,@@7DGEFEF6F3CCG>E5EC7
F:C:EF@<GDEDC*FF?DGFGEGDFFGEGFBFFFGF;GGBGG+GGD?GFCG<GG=G4GGEFFEGG@EBGFGGFEGGGGG?CEGGEGGGGGFFGGGFG
GGDGGBGGGDGAGGCGF$GDG AS:i:-8 XM:i:0 XO:i:0 XG:i:0 MD:Z:151 NM:i:0 XS:A:+ NH:i:1
Frag_1 81 NC_023750.1 41454868 50 73M411N78M = 41452501 -
2929 TGGTCATGATGCTCATGTGGCTATGCTTCTTGGTGCTGCGAAGATTCTTAAACAGCATGAAAAAGAGATACAAGGAACTGTAGTTCTTG
TTTTCCAACCAGCAGAGGAAGGAGGTGGGGGAGCCAAGAAAATTTTAGATGCTGGAGCCTTA 2>A=<*??B=C57FCD,CABD?C=BFBBCEAB)
B9F:FGFGF9FCGFCGGGAGFAEFGGGF@GEGAFG*GFECGGFGG0EGGEGGG?G.FGFEGFGFGGEGEFGGGFGFGFFGGEDF6GEFGGGBGGGGG
GFEFFGG1BGEGGGAGDFGAG AS:i:-8 XM:i:0 XO:i:0 XG:i:0 MD:Z:151 NM:i:0 XS:A:+ NH:i:1
q
Isso é pra usar as referências do gft que já foram mapeadas dos introns para facilitar a identificação dos introns
```
tophat2 --GTF ../../ref/genome.gtf genome ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq
tophat2 --b2-very-sensitive --transcriptome-index ./transcriptome --GTF ../ref/genome.gtf genome ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq
tophat2 --b2-very-sensitive --transcriptome-index ./transcriptome --library-type fr-firststrand --GTF ../ref/genome.gtf genome ../raw/SAMPLEC1_R1.fastq ../raw/SAMPLEC1_R2.fastq -> aqui não fez muita diferença. Ver de novo, rodar para ver o que dá.
o fr-firststrand só dá certo se for a biblioteca TruSeq da Illumina, senão ele lê errado.
Quando é o tophat, o alinhamento é em relação à fita senso, quando comparado no transcriptoma. No genoma isso não funciona, porque nao tem senso e antisenso.
rf no bowtie <- ->
fr -> <- (eu vou usar essa aqui, que é pair-ended da illumina)
o bedtools é interessante para selecionar coisas específicas em arquivos bed, por exemplo a expressão de genes específicos usando as coordenadas do arquivo gff ou gtf
O tophat verifica a predição de introns quando uma parte da read se alinha parcialmente e a outra se alinha em dois pedaços diferentes
introntab.pl -> script para ver introns. Ele não é tããão bom
o outro que vamos usar -> checkMinMaxIntronSize.sh
ele dá algumas informações do genoma como saída também.
```
checkMinMaxIntronSize.sh genome.gff .
```
> (python2) bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/ref$ checkMinMaxIntronSize.sh genome.gff .
* Gene structure statistics from current genome reference annotation (introntab.pl)
Running introntab.pl based on gff
13,7102
No caso o intron mínimo é 13 e o tamanho máximo é 7102 nt. Não existe certo e errado.
COm esse script dá para fazer uma descrição básica de genoma.
> (python2) bmoreno@spiderman:/data/bioaat/bmoreno/pvulgaris/ref$ checkMinMaxIntronSize.sh genome.gff . 0.05 0.9
* Gene structure statistics from current genome reference annotation (introntab.pl)
79,1108
Dá pra mudar os critérios o 0.05 e 0.9 são os percentis do tamanho dos introns encontrados no meu gff e espécie.
Para tirar a redundância a lista tem que estar em ordem SEMPRE!! por isso sempre usar o sort numa lista para ver valores únicos. Na lista é pra ter no máximo dois fragmentos que se alinham, quando o seq é pair-ended.
# FORMATO SAM
inicia por @ também e tem um cabeçalho
samtools view -h aparece o cabeçalho
samtools view -H aparece SÓ o cabeçalho
R1 e R2 tem o mesmo nome.
4ª coluna – onde inicia o pareamento com a referência, tipo 7, é o inicio do pareamento da read com a 7ª base da sequência de referência.
Existe um phred de alinhamento/mapeamento. Tem a mesma ideia do sequenciamento. Valor de 30 é um numero bom, cada programa tem o seu valor bom, sempre ver isso.
Soft clipping – se desconsidera uma parte das bases, não as exclui
Hard clipping – exclui as bases não pareadas.
N – splicing. Foram ignoradas no mapeamento por parecerem introns, sempre aparecem no meio da read.
O tamanho do fragmento é da primeira base do R1 até a última base do R2. A 12ª coluna é opcional
No plieup é que a gente tira as informações de SNPs.
Bitwise flag - por meio de numeros binários a gente tem informações de como foi feita o alinhamento, tipo se alinhou, se alinhou corretamente, se é forward ou reverse, alinhamentos primários (aqueles com maior score, ou mais confiável), duplicata PCR. 0 indica falso, não é verdadeiro e 1 é verdadeiro.
Tem um site que decifra isso .
Pacote do trinity
Programa IGV ajuda a visualizar o alinhamento feito.
o bam é a versão binária do Sam
Mexendo no IGV
carregar o genoma com genoma.fa na aba genome
depois colocar o bam e o index bam. Ainda para ver certinho no nosso caso tem que pegar o nome do gene para "filtrar"
![Uploading file..._ay6w8ze5h]()
![Uploading file..._5l9tch0w1]()
(falta um pedaço no libre office)
Distribuição dos dados
rna-seq normalmente é uma variável discontínua
expressão gênica diferencial (DEG genes diferentemente expressos)
Repetição biológica!!! Te dá margem para calcular a variância e média com mais precisão.
O tamanho do transcrito importa - viés de tamanho
NOrmalização dos dados
fator de normalização dos dados tem problemas com amostras com diferenças de profundidade muito significativas. N/T (N- fator de normalização tipo a soma total deu 25, e eu quero transformar num total de 100 ou 1000).
Ajuste também deve ser feito pelo tamanho do read, porque um transcrito maior vai precisar de mais reads para juntar e formar 1. então o RPKM é uma união do F/T + # de reads/tamanho do mRNA.
FPKM = RPKM a diferença é que o FPKM é para single end, 1 leitura só do fragmento. Isso muda a amostragem do mRNA.
Cufflink - ele já não é o melhor programa. Existe o Cuftie? (não entendi o nome)
ele pode fazer tudo do zero
usar o gff e o que tem de dados do experimento
restringir só a gff
o star é muito bom. hisat2 tem é bom e parece o tophat.
grafos de sobreposição - cada nó representa uma read e a aresta conectando ela indica que elas são alinhadas (indica que elas são mutualmente exclusivas. Os diferentes caminhos que pode se formar são as isoformas.
As reads em comum são distribuídas de forma estatística (modelo de máxima verossimilhança) para as isoformas mais abundantes. Porque é o esperado, que as mais abundantes tenham mais desses éxons.
Ele também pode fazer a montagem de novo - isso é possível porque ele cria falsas reads a partir do gff para ajudar a criar profundidade suficiente para que seja identificados os exons de forma correta.
rPKM é diferente do TMP porque ele cria a medida por milhão primeiro, depois ele normaliza por tamanho. O TPM, por sua vez, divide o tamanho do gene por reads.
o DESeq e o edgeR fazem a normalização Log Geometric Mean (calcular a média dos loges, que não é igual a tirar a média dos valores normais, sem o loge) que não leva em consideração o tamanho das reads, porque não faz muito sentido uma vez que eu comparo tamanho de reads e sim genes. O DESeq assume que os dados não estão normalizados, é importante ver isso antes de fazer qualquer análise. É uma normalização por gene não total de genes de um sample.
https://www.youtube.com/@statquest -> assistir mais vídeos dele.
--library-type é igual ao tophat!!
o cuff norm normaliza os dados de maneira default. Tem como desnormalizar para usar no DESeq
cuff merge = une diferentes transcriptomas ele cria uma referência unindo as reads das diferentes replicas.isso para um gene.
cuffdiff - é um teste parecido com o teste -t e tem a distribuição do log2fold change
-> false discovery rate - múltiplos testes de hipótese. correção por Bonferroni (aula de biometria)
ele mostra os diferentes TSSs das diferentes isoformas. Ele dá também a expressão do cds, uso de isoformas.
o p-ajustado é a aplicação da correção de falsos positivos igual a correção de bonferroni nos SNPs. A diferença é que o FDR é um método mais brando que o bonferroni. Na tabela da Nicole a quantidade de genes diferentemente expressos na verdade é aqueles que o p-ajds é menor ou igual à 0,05.
saídas do cufflink
o transcriptoma tem um código para cada isoforma:
class code = códigos das isoformas, mostra o quanto de match deu de introns em quantidade e qualidade e outras características.
colocar as bases em minúsculo mostra regiões repetitivas (isso também se chama de soft-masking).
string tie = igual ao cufflink, o que muda é o algoritmo que monta as isoformas.
Illumina iGenome - tem espécies boas com bastante informações sobre o genoma de referência. Deixa tudo bonitinho, como nós deixamos os dados arrumados para trabalhar.
**Pre-processamento dos dados
sempre fazer uma análise exploratória dos dados para ver se a facility não trocou as amostras. Principalmente uma análise hierárquica.
estrutura de diretórios
- [ ] -> colocar numero na sequencia dos programas que foram usados e as pastas com os nomes dos programas utilizados para cada etapa.
raw = dados exatamente do jeito que veio da facility
input = colocar todos os dados preparados para rodar.
dados alfanuméricos (nome igual é um grupo de amostras e o número a réplica biológica) para rodar o pipeline do professor. Ele também só roda paired-end. SEMPRE salvar os arquivos intermediários até a publicação pelo menos, para evitar problemas no momento da publicação.
output -> processed (melhor nome seria pre-process) e analysis.
para tirar adaptador - atropos (mais recente)
qualidade de reads - prinseq
-> repositórios bons
SRA - NCBI
ENA - EBI
DAR - DDBJ
para transformar sra em fastq fastq-dump ou sff-dump (quando é dado ROche 454)
`fastq-dump --slit3 nomedoaquivo.sra - separa os arquivos em R1 e R2 e transforma .sra em .fastq`
tem como fazer reamostragem, criar um subconjunto de dados reduzido
análise de qualidade no fastq - tem vários programas que fazem isso.
FASTX-Toolkit
FastQC
Prinseq
assemblystats
k-mer -> subsequencias
o fastqc é bom para sequenciamento genômicos, o RNA-seq normalmente não segue a distribuição normal da porcentagem de GC. Então GC a gente não leva em consideração. A análise de k-mers indica se há ou não a presença de adaptadores. Vendo os k-mers se espera que sua distribuição seja regular, se não for tem algo errado.
remoção de cauda poli-A -> prinseq-lite.pl
fita 1 tem o poli-T, porque é complementar do DNA
se faz trimagem na read do lado direito porque normalmente é o lugar que as bases tem valor menor.
a complexidade é inversamente proporcional. Valores altos, baixa complexidade, ou seja, só um nucleotídeo (TTTTT, AAAAA). Quanto mais complexo, menor o valor (TAGCCCGAT) tem valores baixos tipo 30, 20. Acima de 7 nt é de baixa complexidade.
Algoritmos de trimagem de adaptador
Scythe - trima só 3'.ele vai classificar do match com base na qualidade de bases para decidir se é adaptador ou não.
Cutadapt parece o atropos. tira adaptador de qualquer lugar da read. Mais sensível.
https://pt.slideshare.net/dgpinheiro/illuminatruseqpairedend
(gunzip -c é um descompactamento temporário para fazer uma análise rápida)
o atropos é um pouco diferente. ele funciona igual ao trimmomatic. Correção de bases no atropos dá um pouco de problema
atropos -> fastqc
trimagem -> avaliação de qualidade.
**ctrl +a esc - deixa a gente caminhar pra cima e pra baixo**
o problema do nextseq é que por usar duas cores de fluoroforo o G, que não tem cor, ser erroneamente medido, porque lugares mal mapeado também ficam sem cor. Sendo que pode ser erro mesmo ou Gs que não foram mapeados.
Por fim, vamos começar tudo de novo pra poder rodar o script de maneira completa!!
levamos tudo o que fizemos pro test e comecei uma nova pasta de pvulgaris
criamos o raw e copiamos tudo de abundância e script de pegar as sequências do ncbi.
todo script tem que colocar uma validação, tipo mostrar uma mensagem se tem erro, ou se tá faltando alguma coisa para o script.
o readlink parece o pwd, ele mostra por meio do link relativo (. ou ..) o absoluto (todo o caminho que leva até lá)
git clone link do github -> como faz para pegar um rpograma no github
#! indicador para quem vai executar esse script
a validação ajuda bastante para gente quando depois de um tempo se esquecemos alguma coisa.
./rnaseq-ref.sh star ./input ./output 3 ./ref/genome.gtf ./ref/genome.fa stringtie
Para ver as mensagens de erro, os logs em geral do processamento.
```
cat ./output/star_index/STAR.genomeGenerate.log.*txt
```
checkminandmaxintronsize = rodar ele também para saber
min anchor= minimo de bases que vai ancorar no exon à esquerda e o resto tem que alinhar à direita
o segundo tophat fica unstranded, porque é o alinhamento dos singletons que pode ser R1 ou R2. Por isso não pode ter sentido correto.
alt + u = desfazer (undo)
alt + e = refazer
tag XS é a última coluna no arquivo sam, aquelas quesão opcionais.
Caso os alinhamentos não funcionem direito, é importante ir opção por opção verificar o que pode ser adicionado ou retirado dependendo da situação.
--alignSJDBoverhangMin = esse aqui funciona quando se conhece os introns e tem que alinhar pelo menos 1 nem um exon e o resto no outro exon.
--pre-mrna-fraction - o quanto de reads que alinham dentro dos introns que pode ser "deixada" lá.
o que é transfrag ????
--intron-overhang-tolerance não considera que haja introns no meio dos exons no mRNA.
cuffcmp.combined.gtf - mostra o código e classe das isoformas. e tem um script getTransFragInfo.pl
--st=${outdir}/${aligner}_${assembler}_cuffnorm/samples.table - é aqui que tem o fator de normalização, para poder fazer a desnormalização.
deu erro, daí ele precisa ser compilado de novo. Ele foi e tem que usar python2
```
conda run -n python2 no cuffnorm, cuffdiff, tophat (nos dois) -> atualização para dar certo a "rodagem"
./rnaseq-ref.sh star ./input/ ./output/ 3 ./ref/genome.gtf ./ref/genome.fa stringtie
```
samtools view arquivo.
depois retiramos o conda porque o professor fez um arquivo tipo script que ativava o conda e usava os parametros da linha de comando
gene_exp.diff é o arquivo que tem a diferença de expressão.
o q_value é o p_value corrigido pelo FDR , o padrão é 0,05 global.
No cuffdiff há um problema de isoforma para calcular o fpkr. No exemplo do professor deu errado.
programa tipo iDR (programa que te ajuda a trabalhar com uma linguagem específica).
para acessar o Rstudio do servidor
http://200.145.102.33:8787/
modo slave do R no terminal do linux ele executa apenas o comando, se ter que abrir o R e deixar ele aberto no terminal
modo slave --slave.
para chamar um script em R no linux
#!/usr/bin/env Rscript no heading do script, igual o do bash!!
objetos podem ter atributos ou ações associadas a eles.
complexos 1. matrix só armazena um tipo de variável, 2. data.frame armazena mais de um tipo de variável 3. factor - são vetores categóricos 4. array tem mais dimensões do que linhas e colunas 5. function - é uma função
simples/básicos - 1. vetor atômico (1 valor - escalares ou vetores - vários valores em uma só coluna, ou de verdadeiro e falso) 2. lista - tem propriedade de ter uma chave e um valor e onsegue misturar tipos de dados, tipo caracter e números. Indexado por caracetres e números 3.Null - tem valor nulo.
`l <- list ("Daniel" =43, "Luis" = 50)`
<<- atribui não só naquele script como em outros
Tem como atribuir nomes aos elementos de uma matrix ou de vetores
não é possível transformar um factor em numeric, porque ele funciona de maneira diferente.
table() - faz a contagem de cada nível no R.
as.matrix - também existe, transforma tabelas em matrizes.
df[[]] é igual ao df$col. Existem essas duas maneiras de acessar as colunas.
o sep é diferente do collapse porque o sep são vetores escalares. Para juntar vetores normais é collapse.
O scan também serve para abrir files no R.
`gene.df <- read.delim(file = "/data/bioaat/bmoreno/pvulgaris/output/star_stringtie_cuffdiff/gene_exp.diff",
header = TRUE,
sep= "\t",
dec =".",
stringsAsFactors = FALSE)
dim (gene.df)
`
> > gene.df <- read.delim(file = "/data/bioaat/bmoreno/pvulgaris/output/star_stringtie_cuffdiff/gene_exp.diff",
+ header = TRUE,
+ sep= "\t",
+ dec =".",
+ stringsAsFactors = FALSE)
> View(gene.df)
> dim (gene.df)
[1] 1721 14
`colnames(gene.df)
signif.gene.df <- subset(gene.df, q_value <= 0.05)`
> > colnames(gene.df)
[1] "test_id" "gene_id" "gene" "locus" "sample_1"
[6] "sample_2" "status" "value_1" "value_2" "log2.fold_change."
[11] "test_stat" "p_value" "q_value" "significant"
> signif.gene.df <- subset(gene.df, q_value <= 0.05)
`signif.gene.df[,c("gene", "log2.fold_change.", "q_value")]`
> > signif.gene.df[,c("gene", "log2.fold_change.", "q_value")]
gene log2.fold_change. q_value
1 - 5.912350 7.5e-05
3 - 1.679970 7.5e-05
4 - -0.594616 6.1e-03
`iris.df <- read.delim(file = "/tmp/iris.data",
header = FALSE,
sep= ",",
dec =".",
stringsAsFactors = FALSE)
colnames(iris.df) <- c("sepal.length", "sepal.width", "petal.length", "petal.width", "class")
iris.df$class <- as.factor(iris.df$class)
table(iris.df$class)`
> > iris.df <- read.delim(file = "/tmp/iris.data",
+ header = FALSE,
+ sep= ",",
+ dec =".",
+ stringsAsFactors = FALSE)
> colnames(iris.df) <- c("sepal.length", "sepal.width", "petal.length", "petal.width", "class")
> View(iris.df)
> iris.df$class <- as.factor(iris.df$class)
> table(iris.df$class)
Iris-setosa Iris-versicolor Iris-virginica
50 50 50
AS DUAS BARRAS NO CAMINHO É PRA PROTEGER AS BARRAS!!! NÃO SABIA DISSO.
`library("xlsx")
write.xlsx(gene.df, "/data/bioaat/bmoreno/pvulgaris/gene_exp_diff.xlsx", "gene_exp_diff.xlsx")`
>
delim, subset e write.xlsx = mais importante.!!!
`write.table(signif.gene.df, /data/bioaat/bmoreno/pvulgaris/signif_gene_exp_diff.txt, sep = "\t", quote=FALSE)`
> ##### Testando sozinha
dados.df <- read.delim(file = "/data/bioaat/bmoreno/pvulgaris/output/star_stringtie_cuffdiff/gene_exp.diff",
header = TRUE,
sep= "\t",
dec =".",
stringsAsFactors = FALSE)
diff.exp <- subset (dados.df, dados.df$log2.fold_change. >= 2)
write.xlsx (diff.exp, "/data/bioaat/bmoreno/pvulgaris/teste_gene_exp.xlsx", "teste_gene_exp.xlsx")
>
lalala
deu certo:
[planilha certa](https://hackmd.io/_uploads/r1X1t_qEa.png)
apply() - função que quase não exite em outras linguagens. Aplicação de uma função em um conjunto de dados.
apply (iris.df[,c("sepal.width", "sepal.length", "petal.width", "petal.length")], 2, mean)
> apply (iris.df[,c("sepal.width", "sepal.length", "petal.width", "petal.length")], 2, mean)
sepal.width sepal.length petal.width petal.length
3.054000 5.843333 1.198667 3.758667
###exercício aula slides
golub <- read.delim("http://www.bioinf.ucd.ie/people/ian/Golub.txt", sep="\t")
head(golub)
colnames(golub)
library(limma)
golub[,c(samples.AML,samples.ALL)] <- normalizeQuantiles(golub[,c(samples.AML,samples.ALL)])
samples.AML <- colnames(golub)[ grep("^AML\\.", colnames(golub)) ]
samples.ALL <- colnames(golub)[ grep("^ALL\\.", colnames(golub)) ]
golub$ALL.mean <- apply(golub[, samples.ALL], 1, mean)
golub$AML.mean <- apply(golub[, samples.AML], 1, mean)
golub$logFC <- log2(golub$AML.mean/golub$ALL.mean)
calc.p.value <- function(v, x, y) {
if ( sd(v[c(x,y)]) == 0 ) {
return(1)
} else {
return(t.test(v[x], v[y])$p.value)
}
}
golub$p.value <- apply(golub[,c(samples.AML,samples.ALL)], 1, calc.p.value, samples.AML, samples.ALL)
golub$q.value <- p.adjust( golub$p.value , method="fdr" )
signif.golub.df <- subset(golub, q.value<= 0.05 & abs(logFC) >= log2(1.5))
signif.golub.df[, c("Golubtemplate", "AML.mean", "ALL.mean", "logFC", "q.value")]
> > golub <- read.delim("http://www.bioinf.ucd.ie/people/ian/Golub.txt", sep="\t")
> head(golub)
Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5 ALL.6 ALL.7 ALL.8 ALL.9
1 AFFX-BioB-5_at 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 AFFX-BioB-M_at 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 AFFX-BioB-3_at 6.643856 8.033423 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 7.409391 8.285402
4 AFFX-BioC-5_at 8.357552 8.204571 8.108524 7.507795 8.055282 6.643856 7.149747 7.977280 7.442943 7.149747
5 AFFX-BioC-3_at 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.832890
6 AFFX-BioDn-5_at 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
ALL.10 ALL.11 ALL.12 ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18 ALL.19 AML AML.1
1 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
4 7.467606 7.011227 7.189825 7.721099 7.491853 6.643856 7.888743 7.247928 8.076816 7.426265 6.643856 8.271463
5 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
6 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
AML.2 AML.3 AML.4 AML.5 AML.6 AML.7 AML.8 AML.9 AML.10 AML.11 AML.12 AML.13
1 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 6.643856 7.727920 6.672425 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
4 7.409391 8.294621 8.317413 7.475733 8.409391 6.643856 8.629357 7.139551 6.643856 7.189825 8.751544 6.643856
5 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
6 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
ALL.20 ALL.21 ALL.22 ALL.23 ALL.24 ALL.25 ALL.26 ALL.27 ALL.28 ALL.29 ALL.30 ALL.31
1 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 6.643856 6.643856 6.643856 8.049849 6.643856 7.748193 7.894818 6.643856 6.727920 6.643856 6.643856 6.643856
4 6.643856 8.144658 8.271463 6.643856 7.392317 6.643856 6.643856 6.643856 8.066089 7.774787 6.643856 6.643856
5 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
6 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
ALL.32 ALL.33 ALL.34 ALL.35 ALL.36 ALL.37 ALL.38 ALL.39 ALL.40 ALL.41 ALL.42 ALL.43
1 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
4 7.434628 7.924813 7.219169 8.005625 8.233620 6.643856 7.693487 7.768184 7.044394 6.643856 6.643856 7.994353
5 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
6 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
ALL.44 ALL.45 ALL.46 AML.14 AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21 AML.22
1 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
2 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
3 6.643856 6.643856 6.643856 6.643856 6.643856 7.087463 6.954196 6.643856 6.643856 6.643856 6.643856 6.643856
4 8.303781 6.643856 6.643856 6.658211 7.044394 8.312883 8.344296 8.614710 7.912889 7.592457 8.285402 7.845490
5 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
6 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
AML.23 AML.24
1 6.643856 6.643856
2 6.643856 6.643856
3 6.643856 6.643856
4 8.366322 8.396605
5 6.643856 6.643856
6 6.643856 6.643856
> colnames(golub)
[1] "Golubtemplate" "ALL" "ALL.1" "ALL.2" "ALL.3" "ALL.4"
[7] "ALL.5" "ALL.6" "ALL.7" "ALL.8" "ALL.9" "ALL.10"
[13] "ALL.11" "ALL.12" "ALL.13" "ALL.14" "ALL.15" "ALL.16"
[19] "ALL.17" "ALL.18" "ALL.19" "AML" "AML.1" "AML.2"
[25] "AML.3" "AML.4" "AML.5" "AML.6" "AML.7" "AML.8"
[31] "AML.9" "AML.10" "AML.11" "AML.12" "AML.13" "ALL.20"
[37] "ALL.21" "ALL.22" "ALL.23" "ALL.24" "ALL.25" "ALL.26"
[43] "ALL.27" "ALL.28" "ALL.29" "ALL.30" "ALL.31" "ALL.32"
[49] "ALL.33" "ALL.34" "ALL.35" "ALL.36" "ALL.37" "ALL.38"
[55] "ALL.39" "ALL.40" "ALL.41" "ALL.42" "ALL.43" "ALL.44"
[61] "ALL.45" "ALL.46" "AML.14" "AML.15" "AML.16" "AML.17"
[67] "AML.18" "AML.19" "AML.20" "AML.21" "AML.22" "AML.23"
[73] "AML.24"
> library(limma)
> samples.AML <- colnames(golub)[ grep("^AML", colnames(golub)) ]
> samples.ALL <- colnames(golub)[ grep("^ALL", colnames(golub)) ]
> samples.AML <- colnames(golub)[ grep("^AML\.", colnames(golub)) ]
Error: '\.' is an unrecognized escape in character string (<input>:1:44)
> samples.ALL <- colnames(golub)[ grep("^ALL\.", colnames(golub)) ]
Error: '\.' is an unrecognized escape in character string (<input>:1:44)
> samples.AML <- colnames(golub)[ grep("^AML\\.", colnames(golub)) ]
> samples.ALL <- colnames(golub)[ grep("^ALL\.", colnames(golub)) ]
Error: '\.' is an unrecognized escape in character string (<input>:1:44)
> samples.ALL <- colnames(golub)[ grep("^ALL\\.", colnames(golub)) ]
> golub$ALL.mean <- apply(golub[, samples.ALL], 1, mean)
> golub$AML.mean <- apply(golub[, samples.AML], 1, mean)
> golub$logFC <- log2(golub$AML.mean/golub$ALL.mean)
> golub$ALL.mean <- apply(golub[, samples.ALL], 1, mean)
> golub$AML.mean <- apply(golub[, samples.AML], 1, mean)
> golub$logFC <- log2(golub$AML.mean/golub$ALL.mean)
> calc.p.value <- function(v, x, y) {
+ if ( sd(v[c(x,y)]) == 0 ) {
+ return(1)
+ } else {
+ return(t.test(v[x], v[y])$p.value)
+ }
+ }
> golub$p.value <- apply(golub[,c(samples.AML,samples.ALL)], 1, calc.p.value, samples.AML, samples.ALL)
> golub$q.value <- p.adjust( golub$p.value , method="fdr" )
> View(golub)
> subset(golub, q.value<= 0.05)
Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5 ALL.6 ALL.7
48 AFFX-HUMTFRR/M11507_5_at 7.219169 7.768184 8.957102 7.988685 7.761551 8.515700 7.312883 9.403012
50 AFFX-HUMTFRR/M11507_3_at 7.839204 6.643856 9.746514 7.754888 7.794416 7.108524 6.643856 9.634811
65 AB000449_at 7.219169 6.643856 8.144658 8.154818 9.702173 7.569856 8.463524 8.539159
88 AB002559_at 10.245553 10.076816 10.040290 9.396605 9.467606 10.299208 10.037547 10.770664
92 AB003698_at 6.845490 7.451211 7.562242 7.312883 9.461479 8.124121 8.219169 7.636625
98 AC000064_cds1_at 7.820179 9.303781 7.982994 7.965784 8.651052 9.041659 7.098032 8.159871
119 AF000231_at 6.643856 6.643856 6.643856 7.022368 7.707359 6.643856 7.199672 6.643856
131 AF002700_at 7.807355 8.607330 8.791163 8.361944 8.154818 9.483816 6.643856 8.900867
134 AF005043_at 6.643856 6.643856 6.643856 6.977280 7.761551 6.643856 6.643856 6.643856
136 AF005775_at 9.299208 10.779719 8.169925 8.479780 8.344296 9.544964 8.011227 9.211888
149 AF009426_at 6.643856 6.643856 6.643856 6.643856 7.266787 6.643856 6.930737 7.303781
155 AF015913_at 9.419960 8.577429 9.144658 9.052568 10.483816 9.850187 9.881114 8.370687
ALL.8 ALL.9 ALL.10 ALL.11 ALL.12 ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18
48 6.882643 7.129283 6.643856 7.266787 7.483816 7.851749 7.924813 7.622052 8.060696 7.247928 8.257388
50 6.643856 6.643856 7.276124 7.011227 8.252665 7.768184 7.189825 7.864186 8.787903 8.033423 7.988685
65 7.960002 8.103288 6.845490 6.643856 6.643856 6.643856 8.438792 8.577429 6.845490 6.643856 8.005625
88 10.366322 9.511753 9.688250 10.066089 9.596190 9.430453 10.127994 9.758223 9.139551 9.903882 10.427313
92 7.851749 7.209453 6.643856 6.643856 6.643856 6.643856 8.066089 7.247928 6.643856 7.266787 8.294621
98 7.179909 6.988685 7.523562 6.965784 7.507795 7.219169 8.832890 8.471675 8.321928 8.087463 7.629357
119 6.643856 6.643856 6.643856 6.643856 6.643856 7.339850 8.049849 7.592457 7.209453 6.643856 6.643856
131 8.842350 8.614710 7.546894 8.164907 8.647458 8.596190 9.505812 7.199672 7.515700 8.438792 8.640245
134 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.930737 6.700440 6.643856 6.643856 6.643856
136 7.033423 8.658211 8.149747 7.366322 9.370687 9.438792 8.546894 8.665336 9.507795 8.826548 8.693487
149 7.434628 6.643856 6.643856 6.643856 6.643856 6.643856 7.357552 6.643856 6.643856 6.643856 6.643856
155 9.405141 9.507795 9.893302 8.577429 8.991522 8.562242 10.716819 11.046442 7.826548 9.831307 8.870365
ALL.19 AML AML.1 AML.2 AML.3 AML.4 AML.5 AML.6 AML.7 AML.8
48 6.700440 7.118941 8.693487 8.876517 8.507795 7.483816 7.118941 7.813781 7.714246 8.294621
50 6.643856 6.643856 9.926296 9.592457 9.479780 7.754888 6.643856 8.939579 7.614710 8.252665
65 7.607330 6.643856 6.643856 6.643856 6.643856 7.491853 6.643856 6.643856 7.651052 6.643856
88 10.406205 9.743151 10.870365 10.637531 10.952013 11.381543 10.139551 11.814182 11.072133 10.646559
92 7.679480 6.643856 7.400879 6.643856 6.643856 8.224002 6.965784 6.832890 7.055282 6.643856
98 7.515700 7.139551 9.413628 9.182394 9.187352 9.098032 8.357552 8.581201 8.417853 9.579316
119 7.199672 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
131 7.607330 7.330917 9.014020 9.255029 9.152285 8.447083 8.238405 8.603626 8.897845 9.121534
134 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
136 8.700440 8.455327 9.890264 10.393390 9.794416 9.537218 8.848623 7.948367 8.413628 9.688250
149 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
155 8.299208 7.339850 9.027906 7.971544 8.836050 9.670656 9.060696 8.592457 8.038919 8.463524
AML.9 AML.10 AML.11 AML.12 AML.13 ALL.20 ALL.21 ALL.22 ALL.23 ALL.24 ALL.25
48 7.748193 6.643856 8.238405 8.731319 8.417853 8.044394 9.002815 9.250298 7.515700 8.184875 7.936638
50 7.247928 7.434628 7.741467 7.894818 6.643856 8.219169 8.154818 7.546894 8.022368 8.164907 7.794416
65 6.643856 7.139551 8.016808 7.392317 6.643856 6.643856 7.400879 8.487840 8.098032 8.280771 7.857981
88 9.586840 9.802516 11.295195 11.503329 11.068106 9.567956 10.693487 10.142107 10.052568 9.682995 9.994353
92 6.794416 7.118941 7.754888 7.607330 6.643856 7.011227 7.507795 8.724514 7.199672 7.924813 7.087463
98 8.531381 7.882643 8.807355 10.060696 8.388017 8.082149 9.216746 8.836050 8.693487 8.632995 7.179909
119 6.643856 7.000000 7.209453 6.643856 6.643856 7.400879 6.643856 6.643856 6.643856 7.614710 6.643856
131 8.396605 7.936638 9.836050 9.060696 7.864186 8.543032 8.906891 9.741467 8.982994 7.734710 8.693487
134 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 7.348728 7.118941 6.741467 6.643856
136 9.157347 8.845490 9.396605 9.972980 10.250298 8.816984 9.337622 9.717676 9.052568 9.005625 8.927778
149 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.906891 6.643856
155 8.071462 8.977280 9.440869 8.686501 7.599913 8.894818 9.255029 9.731319 9.432542 10.348728 8.434628
ALL.26 ALL.27 ALL.28 ALL.29 ALL.30 ALL.31 ALL.32 ALL.33 ALL.34 ALL.35 ALL.36
48 7.417853 8.927778 8.535275 6.643856 8.044394 7.888743 7.523562 7.312883 7.434628 7.066089 8.936638
50 7.139551 6.845490 8.335390 8.321928 7.707359 7.607330 7.366322 8.700440 8.144658 7.257388 8.154818
65 7.033423 6.643856 8.290019 7.672425 8.262095 7.303781 9.923327 8.299208 8.960002 6.965784 8.515700
88 10.060696 10.506803 10.101976 10.184875 9.370687 9.556506 9.507795 9.842350 9.881114 9.330917 11.400346
92 6.781360 6.643856 8.554589 7.686501 8.348728 6.643856 8.873444 8.543032 7.774787 7.108524 9.022368
98 8.962896 9.049849 8.262095 8.252665 7.636625 8.463524 8.209453 9.236014 8.592457 8.224002 8.479780
119 7.189825 6.643856 7.044394 6.643856 6.672425 6.643856 8.194757 6.643856 6.714246 8.000000 6.643856
131 9.592457 9.707359 8.810572 7.000000 7.247928 8.503826 7.651052 8.113742 7.129283 7.076816 9.974415
134 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 7.238405 6.643856 6.643856 6.643856 6.781360
136 8.654636 9.390169 8.942515 8.721099 8.455327 8.459432 8.839204 8.994353 8.727920 9.036174 9.147205
149 6.643856 7.400879 6.643856 6.643856 6.643856 6.643856 6.942515 6.643856 8.066089 7.417853 6.643856
155 9.350939 9.103288 9.689998 8.882643 8.778077 7.900867 10.605480 9.461479 9.879583 9.129283 10.415742
ALL.37 ALL.38 ALL.39 ALL.40 ALL.41 ALL.42 ALL.43 ALL.44 ALL.45 ALL.46 AML.14
48 8.247928 7.348728 8.199672 6.643856 8.154818 8.409391 7.960002 8.348728 7.426265 7.238405 8.417853
50 8.894818 7.266787 8.044394 7.523562 7.357552 7.665336 6.643856 9.252665 7.303781 7.651052 8.519636
65 7.366322 8.154818 9.667112 9.060696 6.643856 7.066089 8.569856 7.700440 8.154818 7.622052 6.643856
88 9.839204 9.403012 10.288866 7.108524 10.604553 9.550747 9.668885 9.381543 9.266787 10.213104 10.814582
92 7.475733 7.011227 8.710806 8.455327 6.643856 7.491853 8.700440 6.643856 6.882643 7.228819 6.918863
98 8.174926 7.832890 8.820179 8.087463 8.321928 8.413628 9.808964 7.977280 7.643856 9.771489 9.219169
119 6.643856 7.339850 9.459432 7.665336 6.643856 6.643856 7.936638 6.643856 7.228819 6.643856 6.643856
131 8.285402 8.134426 8.703904 7.577429 8.438792 8.829723 8.499846 8.483816 8.447083 9.927778 9.030667
134 6.643856 6.643856 7.459432 7.614710 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
136 10.278449 8.209453 9.306062 10.016808 8.842350 8.873444 9.442943 8.357552 8.679480 9.627534 9.552669
149 6.643856 6.643856 8.060696 7.686501 6.643856 6.643856 7.076816 6.643856 6.643856 6.643856 6.643856
155 8.527477 9.189825 11.264443 6.643856 7.686501 9.068778 10.234817 8.396605 8.588715 9.066089 8.727920
AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21 AML.22 AML.23 AML.24
48 8.754888 8.071462 9.656425 10.650154 9.453271 9.126704 8.675957 8.409391 8.636625 9.552669
50 8.370687 8.189825 10.741467 11.565578 9.497852 9.665336 8.764872 8.413628 8.330917 10.525521
65 6.781360 6.643856 7.055282 6.794416 7.475733 7.499846 6.643856 6.906891 6.643856 7.994353
88 10.000000 11.590587 10.946906 10.082149 10.090112 9.826548 11.272630 10.675075 10.635718 10.369597
92 7.768184 7.434628 7.247928 6.643856 6.714246 6.643856 6.930737 6.643856 6.643856 6.643856
98 9.596190 8.487840 9.002815 9.243174 8.826548 8.326429 9.066089 9.142107 9.019591 9.219169
119 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 7.033423 6.643856
131 8.731319 9.105909 8.596190 8.554589 9.483816 8.038919 9.415742 9.113742 9.686501 10.062046
134 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
136 9.224002 10.286558 10.176173 8.829723 8.977280 9.592457 9.611025 10.288866 9.317413 9.539159
149 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
155 8.228819 8.686501 8.076816 7.982994 8.731319 8.434628 8.459432 8.888743 8.945444 8.679480
ALL.mean AML.mean logFC p.value q.value
48 7.835278 8.487009 0.11527185 2.762802e-03 2.777999e-02
50 7.750385 8.656344 0.15948980 2.981440e-03 2.911601e-02
65 7.843955 6.996906 -0.16486430 1.500405e-06 8.291771e-05
88 9.863342 10.711407 0.11899959 1.845546e-06 9.892405e-05
92 7.593135 7.023488 -0.11250803 2.449189e-04 4.523386e-03
98 8.245878 8.943189 0.11711628 1.104133e-05 3.987625e-04
119 7.040835 6.698494 -0.07190969 5.449118e-04 8.109972e-03
131 8.413089 8.901810 0.08146316 5.147477e-03 4.252186e-02
134 6.765392 6.643856 -0.02615261 4.831756e-03 4.081231e-02
136 8.899127 9.480508 0.09130059 7.214263e-04 9.947869e-03
149 6.849978 6.643856 -0.04407846 7.498680e-04 1.024726e-02
155 9.235872 8.594995 -0.10375149 3.155276e-04 5.355706e-03
[ reached 'max' / getOption("max.print") -- omitted 900 rows ]
> subset(golub, q.value<= 0.01 & abs(logFC) >= 2)
[1] Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5
[8] ALL.6 ALL.7 ALL.8 ALL.9 ALL.10 ALL.11 ALL.12
[15] ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18 ALL.19
[22] AML AML.1 AML.2 AML.3 AML.4 AML.5 AML.6
[29] AML.7 AML.8 AML.9 AML.10 AML.11 AML.12 AML.13
[36] ALL.20 ALL.21 ALL.22 ALL.23 ALL.24 ALL.25 ALL.26
[43] ALL.27 ALL.28 ALL.29 ALL.30 ALL.31 ALL.32 ALL.33
[50] ALL.34 ALL.35 ALL.36 ALL.37 ALL.38 ALL.39 ALL.40
[57] ALL.41 ALL.42 ALL.43 ALL.44 ALL.45 ALL.46 AML.14
[64] AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21
[71] AML.22 AML.23 AML.24 ALL.mean AML.mean logFC p.value
[78] q.value
<0 rows> (or 0-length row.names)
> subset(golub, q.value<= 0.01 & abs(logFC) > 1)
[1] Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5
[8] ALL.6 ALL.7 ALL.8 ALL.9 ALL.10 ALL.11 ALL.12
[15] ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18 ALL.19
[22] AML AML.1 AML.2 AML.3 AML.4 AML.5 AML.6
[29] AML.7 AML.8 AML.9 AML.10 AML.11 AML.12 AML.13
[36] ALL.20 ALL.21 ALL.22 ALL.23 ALL.24 ALL.25 ALL.26
[43] ALL.27 ALL.28 ALL.29 ALL.30 ALL.31 ALL.32 ALL.33
[50] ALL.34 ALL.35 ALL.36 ALL.37 ALL.38 ALL.39 ALL.40
[57] ALL.41 ALL.42 ALL.43 ALL.44 ALL.45 ALL.46 AML.14
[64] AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21
[71] AML.22 AML.23 AML.24 ALL.mean AML.mean logFC p.value
[78] q.value
<0 rows> (or 0-length row.names)
> subset(golub, q.value<= 0.05 & abs(logFC) > 1)
[1] Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5
[8] ALL.6 ALL.7 ALL.8 ALL.9 ALL.10 ALL.11 ALL.12
[15] ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18 ALL.19
[22] AML AML.1 AML.2 AML.3 AML.4 AML.5 AML.6
[29] AML.7 AML.8 AML.9 AML.10 AML.11 AML.12 AML.13
[36] ALL.20 ALL.21 ALL.22 ALL.23 ALL.24 ALL.25 ALL.26
[43] ALL.27 ALL.28 ALL.29 ALL.30 ALL.31 ALL.32 ALL.33
[50] ALL.34 ALL.35 ALL.36 ALL.37 ALL.38 ALL.39 ALL.40
[57] ALL.41 ALL.42 ALL.43 ALL.44 ALL.45 ALL.46 AML.14
[64] AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21
[71] AML.22 AML.23 AML.24 ALL.mean AML.mean logFC p.value
[78] q.value
<0 rows> (or 0-length row.names)
> subset(golub, q.value<= 0.05 & abs(logFC) >= log2(1.5))
Golubtemplate ALL ALL.1 ALL.2 ALL.3 ALL.4 ALL.5 ALL.6 ALL.7 ALL.8
758 D88270_at 10.970106 10.816184 11.602699 10.909893 13.310187 12.914572 11.922956 11.986909 11.943614
1779 M19507_at 8.936638 6.643856 9.949827 9.463524 7.636625 10.336507 6.643856 10.244364 6.643856
1882 M27891_at 6.643856 6.643856 7.965784 6.741467 6.643856 8.622052 6.643856 6.643856 8.971544
2288 M84526_at 6.643856 6.643856 6.643856 7.665336 6.643856 6.643856 6.643856 6.643856 6.643856
4680 X82240_rna1_at 11.168672 6.643856 6.857981 7.727920 13.965784 12.826548 10.424166 12.735768 12.816184
ALL.9 ALL.10 ALL.11 ALL.12 ALL.13 ALL.14 ALL.15 ALL.16 ALL.17 ALL.18
758 12.763835 11.772315 11.689124 8.276124 10.999295 12.942698 12.667112 7.219169 10.516685 11.827740
1779 6.643856 6.643856 9.840778 6.643856 8.280771 6.643856 9.601771 6.643856 11.147841 8.209453
1882 6.643856 6.643856 7.900867 6.643856 8.823367 7.348728 6.643856 11.029287 8.487840 9.111136
2288 6.643856 6.643856 6.643856 6.643856 6.643856 6.845490 6.643856 7.219169 6.643856 6.643856
4680 12.749031 13.384379 11.942148 8.546894 12.319954 11.027906 13.815583 6.643856 11.583083 13.854576
ALL.19 AML AML.1 AML.2 AML.3 AML.4 AML.5 AML.6 AML.7 AML.8
758 13.002288 6.643856 6.643856 6.643856 6.643856 7.139551 6.643856 6.643856 6.643856 6.643856
1779 11.388556 9.544964 13.136671 13.312174 13.157663 13.965784 13.285691 13.796648 13.965784 13.965784
1882 7.491853 13.965784 13.965784 11.422590 13.512494 13.390303 11.629357 13.792892 11.264443 6.942515
2288 6.643856 12.816384 13.965784 13.266054 12.687157 12.092096 13.304922 12.527233 11.859923 11.056638
4680 13.707790 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 7.599913 6.643856 6.643856
AML.9 AML.10 AML.11 AML.12 AML.13 ALL.20 ALL.21 ALL.22 ALL.23 ALL.24
758 6.643856 7.748193 7.584963 7.700440 6.643856 10.940314 7.629357 9.177420 12.252961 11.724087
1779 13.965784 7.832890 12.525276 11.935533 11.347621 6.643856 6.643856 6.643856 6.643856 9.873444
1882 13.636511 11.038233 13.965784 11.685187 13.965784 8.243174 10.407268 7.988685 6.643856 6.643856
2288 12.188280 6.643856 12.952013 11.401413 12.357827 6.643856 6.643856 6.643856 6.643856 8.179909
4680 6.643856 6.643856 6.643856 6.643856 7.936638 13.824859 6.643856 6.643856 12.346237 12.549303
ALL.25 ALL.26 ALL.27 ALL.28 ALL.29 ALL.30 ALL.31 ALL.32 ALL.33 ALL.34 ALL.35
758 6.643856 10.557464 10.285402 6.643856 8.511753 6.643856 9.475733 12.874405 6.643856 12.397140 11.895197
1779 6.643856 6.643856 6.643856 6.643856 8.614710 6.643856 10.503826 6.643856 6.643856 8.379378 11.798472
1882 10.189825 6.643856 6.643856 8.392317 8.960002 6.643856 6.741467 6.643856 9.679480 6.643856 6.643856
2288 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
4680 6.643856 6.643856 13.241089 6.643856 6.643856 6.643856 8.118941 12.700006 6.643856 13.854479 9.663558
ALL.36 ALL.37 ALL.38 ALL.39 ALL.40 ALL.41 ALL.42 ALL.43 ALL.44 ALL.45
758 6.643856 6.643856 6.643856 10.109831 13.206557 12.079485 6.643856 12.004923 12.321646 12.891784
1779 10.492855 7.523562 8.900867 8.566054 6.643856 6.643856 6.643856 6.643856 7.076816 7.614710
1882 8.550747 6.643856 6.643856 7.679480 6.643856 6.643856 8.731319 6.643856 6.643856 6.643856
2288 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
4680 6.643856 11.606868 11.706496 13.965784 10.049849 11.680799 6.643856 11.771902 11.234817 13.369734
ALL.46 AML.14 AML.15 AML.16 AML.17 AML.18 AML.19 AML.20 AML.21 AML.22
758 11.884171 6.832890 6.643856 6.643856 7.022368 6.643856 7.392317 6.643856 6.643856 7.409391
1779 6.643856 12.764872 13.946815 10.969387 9.970106 12.477252 8.744834 11.105909 8.463524 13.451855
1882 6.643856 10.585901 10.532356 13.965784 13.698705 10.697836 10.067434 12.765286 13.829227 12.821375
2288 6.643856 12.237807 6.643856 12.158925 12.722594 10.615630 6.643856 12.681458 12.847840 12.238703
4680 12.890454 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856 6.643856
AML.23 AML.24 ALL.mean AML.mean logFC p.value q.value
758 7.189825 7.546894 10.531606 6.940867 -0.6015376 5.651072e-14 6.714415e-11
1779 7.707359 12.978710 7.932549 12.032247 0.6010496 2.716821e-10 8.420963e-08
1882 13.590938 11.835261 7.554570 12.275083 0.7003111 1.065609e-13 1.085247e-10
2288 7.139551 10.965784 6.716345 11.383300 0.7611706 4.819829e-10 1.227163e-07
4680 6.643856 6.643856 10.534503 6.737558 -0.6448246 1.024434e-11 4.564495e-09
> dim (subset(golub, q.value<= 0.05 & abs(logFC) >= log2(1.5)))
[1] 5 78
> signif.golub.df <- subset(golub, q.value<= 0.05 & abs(logFC) >= log2(1.5))
> signif.golub.df[, c("Golubtemplate")]
[1] "D88270_at" "M19507_at" "M27891_at" "M84526_at" "X82240_rna1_at"
> signif.golub.df[, c("Golubtemplate", "AML.mean", "ALL.mean", "LogFC", "q.value")]
Error in `[.data.frame`(signif.golub.df, , c("Golubtemplate", "AML.mean", :
undefined columns selected
> signif.golub.df[, c("Golubtemplate", "AML.mean", "ALL.mean", "logFC", "q.value")]
Golubtemplate AML.mean ALL.mean logFC q.value
758 D88270_at 6.940867 10.531606 -0.6015376 6.714415e-11
1779 M19507_at 12.032247 7.932549 0.6010496 8.420963e-08
1882 M27891_at 12.275083 7.554570 0.7003111 1.085247e-10
2288 M84526_at 11.383300 6.716345 0.7611706 1.227163e-07
4680 X82240_rna1_at 6.737558 10.534503 -0.6448246 4.564495e-09
merge() é bom pra melhorar a anotação dos genes nos datas frames já existentes.
svg é vetorial, melhor tipo de extensão para gerar uma imagem.
inkscape é ótimo pra trabalhar com imagens svg no linux.
ggsave também funciona para baixar o gráfico no R.
escala Z score padroniza os resultados quando eles são muito próximos. Maiores que a média é uma cor e menores que a média terão outro valor.
heatdiff<-heatmap.2(
as.matrix(signif.golub.df[,golub.samps]), ##golub.samples são os dados puros
na.rm=TRUE,
*scale="row"*, #aqui é que muda a escala para se tornar z score.
col=greenred,
breaks=100,
key=TRUE,
symkey=FALSE,
symbreaks=TRUE,
density.info="none",
trace="none",
Rowv=T,
Colv=T,
cexRow=0.45,
cexCol=1,
keysize=1,
margins = c(10,10),
dendrogram=c("both"),
main = "HeatMap"
)
metodos de ligação em distância euclediana.
o dendodragrama é uma representação da ligação entre os genes em grupos pela distância. Existem diversas maneiras programar isso e escolher. Em geral se usa a distância euclediana.
aheatmap - anotação do heatmap, tem o pheatmap que é legal.
regexpr - pesquisar mais sobre
aqui dá pra alterar deixar o # de threads fixos, mas aí tem que alterar os parâmetros pra frente na sequencia do input depois.
arquivos .bad no prinseq - são as reads que não passam nos critério definidos antes. No nosso acabou que não gerou, mas pode na vida real.
aqui indica que tem que utilizar os arquivos pre processados que tem _1 e _2 para indicar o foward e reverse. E se o padrão que eu estou usando de dados é diferente eu tenho que mudar aqui.
o -s nessa parte testa para ver se tem ou não. Caso não tenha, ele cria.
**caminhos = caminhos no grafo**, para a construção de um transcrito.
--min_kmer_cov 3 = parâmetro que determina do threshold, o mínimo para considerar o crescimento do transcrito.
rf r1 é reverse do senso e r2 é foward pro senso. do RNA, do transcrito.
./rnaseq-denovo.sh output/processed/cleaned ./output/ 8 20 (8threads e 20G de ram para rodar)
in silico??
não é necessário usar todos os conjuntos de dados, é uma normalização mínima suficiente para que sejam feita as contas. - distribuição de kmers, se n~ão está dentro do esperado é deletado. Isso é, esses dados a mais não vai ser usado na montagem, ele não apaga, só não usa.
dentro do raw
makeblastdb -title
Indice para o blast com base no transcritoma que a gente fez.
"RefTrans" \
-dbtype nucl \
-out RefTrans \
-in transcriptoma.fa \
> makebleastdb.log.out.txt \
2> makeblastdb.log.err.txt
fora do raw
blastn -max_hsps 1 \
-max_target_seqs 1 \
-num_threads 8 \
-query ./output/trinity_assembled/Trinity.fasta \
-task blastn \
-db ./raw/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
hsps - match quase que exato od transcrito com a base de dados. (procurar, não ficou super claro)
-max_target_seqs 1 1 hit - ou seja que seja específico.
-outfmt "6 qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore" \ - formato tabular com as seguintes informações..
O script do detonator tá certo, só a versão dele compilada que está problemática.
crbb - é o melhor hit de uma é igual a da outra, tipo os homólogos do exercício de blast da Claudia.
TransDecoder.LongOrfs -t ./output/trinity_assembled/Trinity.fasta.gene_trans_map -S -- complete_orfs_only --output_dir ./output/
PWM (position weight matrix ou PSSM - position specific scoring matrix.
mergeR.R - une coisas, tabelas no R.
o exepmlo do script trinity abundance funciona só pro trinity, porque tem a sequência certa das colunas.
jeito certo do input
o id tem que ficar o que sai o nome sim pode mudar.
RSEM é um normalizador da estimativa de abundância.
o kalisto gera arquivos tsv. POr algum motivo as versões do Salmon não estão batendo, no trinity não vai direito no novo e o outro não vai bem no velho.
cpm - counts per million
# Para fazer o heatmap:
```
data <- read.delim(file = "/data/bioaat/bmoreno/pvulgaris/output/star_stringtie_cuffdiff/gene_exp.diff",
header = TRUE,
sep= "\t",
dec =".",
stringsAsFactors = FALSE)
```
Sign in with Wallet
Connect another wallet