## Bioinformática II: Análises de transcriptômica 21\10 e 25/10
diretorio a ser usado: cd /state/partition1
porta: hammer 22
via linux: ssh -p 22 hammer.fcav.unesp.br
df -vh: pra saber quanto está disponivel para vc usar, e quanto está sendo usado.
Profundidade de Cobertura: quantas vezes eu estou mostrando aquela mesma base.
Programas diferentes pode ter resultados diferentes, pois a heuristica de cada um é diferente.
Importante: quantidade de ciclos, e profundidade. Eucariotos: no caso de humanos pelo menos 80 milhoes.
Resolucação de extremidades 3' (Sítios de poli(A)): Conseguimos identificar também a orientação ad ORF.
Importante: sequenciamento strand: fita específica (ss)
Quantificação e qualidade
rRNA 28S e 18S (Razao de 2:1) e procariotos: 23S/16S
- Formato de dados:
- FASTA .fasta, .fa .fna
- FASTAQ .fastq ou .fq
- SRA Sequence read archive: .sra
- SRF Sequence read format: .srf
- SFF Sequence flowgram format .sff
- PHD PHreD: .phd ou .phd.1
- SCF Sequence Chromatogram format: .scf
Pacote emboss: Alteração de formatos
sequence.qual: igual ao fasta mais no lugar das bases, há números separados por espaçõs; Qualidade.
Qualidade: Qphred
```
mkdir ./refs
cd ./refs
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.fna.gz
wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/025/095/GCF_003025095.1_Triha_v1.0/GCF_003025095.1_Triha_v1.0_genomic.gff.gz
```
Utilize gunzip para descompactar
`gunzip *.gz`
Limpeza dos cabeçalhos FASTA
Criar arquivo com o script (cleanfasta.sh) abaixo:
`nano cleanfasta.sh`
Ctrl+c & Ctrl+v
```
#!/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_003025095.1_Triha_v1.0_genomic.fna > genome.fa`
Ajuste do arquivo GFF:
Executar o script fixNCBIgff.sh:
`fixNCBIgff.sh GCF_003025095.1_Triha_v1.0_genomic.gff genome.gff`
Script:
```
#!/bin/bash
infile=$1
outfile=$2
if [ ! ${infile} ]; then
echo "Missing input file"
exit
else
if [ ! -e ${infile} ]; then
echo "Wrong input file (${infile})"
exit
fi
fi
infbn=`basename ${infile} .gff`;
if [ ! ${outfile} ]; then
outfile="${infbn}_fixed.gff"
fi
```
```
grep -v -P '\t(cDNA_match|match|repeat_region|region|D_loop|sequence_feature|binding_site|mobile_genetic_element|origin_of_replication|STS|sequence_alteration|Shine_Dalgarno_sequence|ribosome_entry_site|recombination_feature|stem_loop|RNase_P_RNA|telomere|centromere|telomerase_RNA|RNase_MRP_RNA|riboswitch|minus_35_signal|minus_10_signal|protein_binding_site)\t' ${infile} |
grep -v -P '^(##sequence-region|##species|#!genome-build|!processor|#!gff-spec-version)' | \
fixGFFgenestruct.pl | \
fixGFFdupID.pl | \
sortGFF.pl > ${outfile}
```
Importante: gffread - extrai o transcriptoma ou proteoma
gffread -h consulto o help
-g (arquivo da seq fasta)
gff read (nome do arquivo gff)
gffread genome.gff -g genoma.fa -w
-w: gerar a seq, reunir todos os exons por transcrito
se quero a tradução: com o -x eu pego o cds.
gffread genome.gff -g genoma.fa -x
se quero a seq de proteina -y
gffread genome.gff -g genoma.fa -y proteina.fa
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`
para indexar o genoma:
`bowtie2-build -f genome.fa genome_index`
```
bowtie2 --help
```
Returning from initFromVector
Wrote 44026413 bytes to primary EBWT file: genome_index.rev.1.bt2
Wrote 29870612 bytes to secondary EBWT file: genome_index.rev.2.bt2
Re-opening _in1 and _in2 as input streams
Returning from Ebwt constructor
Headers:
len: 119482427
bwtLen: 119482428
sz: 29870607
bwtSz: 29870607
lineRate: 6
offRate: 4
offMask: 0xfffffff0
ftabChars: 10
eftabLen: 20
eftabSz: 80
ftabLen: 1048577
ftabSz: 4194308
offsLen: 7467652
offsSz: 29870608
lineSz: 64
sideSz: 64
sideBwtSz: 48
sideBwtLen: 192
numSides: 622305
numLines: 622305
ebwtTotLen: 39827520
ebwtTotSz: 39827520
color: 0
reverse: 1
```
bowtie2 -x Ref/genome_index -c "GGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTACATT
TGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTATGTTT
GGACATTTATTGTCATTCTTACTCCTTTGTG" > bowtie2-test.txt> bowtie2-test-statistics.txt
```
ou
```
bowtie2 -x Ref/genome_index -c "GGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCTTCATTTGTTATATTGGATACAAGCTTTGCTACGATCTACATT
TGGGAATGTGAGTCTCTTATTGTAACCTTAGGGTTGGTTTATCTCAAGAATCTTATTAATTGTTTGGACTGTTTATGTTT
GGACATTTATTGTCATTCTTACTCCTTTGTG" -S bowtie2-test.txt
```
indexação: contruir um indice para facilitar buscas no genoma.
bowtie2-inspects : recupera a seq fasta a partir do indice, caso vc perca por algum motivo a seq fasta original.
Bowtie:
- end to end: não eh global, somente as extremidades das reads devem ser alinhadas, e nao ambas(reads e ref)
- local: as extremidades das reads podem não alinhar com a referencia.
time antes do comanda mostra o tempo de execução.
quanto menor o tamanho da seed maior a sensibilidade, pois mais regioes vc encontra no genoma.
-N = numero de mismatches que aceita
L= comprimento da seed
`bowtie2 -N 1 -L 12 -i S,1,0.25 -D 30 -R 5 -x ./genome_index -1 /home/jcardoso/sim/raw/SAMPLEA1_R1.fastq -2 /home/jcardoso/sim/raw/SAMPLEA1_R2.fastq -S bowtie2-algA1.txt`
`time bowtie2 -x /state/partition1/jcardoso/Ref/genome -q -1 /home/jcardoso/sim/raw/SAMPLEA1_R1.fastq -2 /home/jacardoso/sim/raw/SAMPLEA1_R2.fastq -S bowtie2-algA1.sam`
25001 reads; of these:
25001 (100.00%) were paired; of these:
12301 (49.20%) aligned concordantly 0 times
12700 (50.80%) aligned concordantly exactly 1 time
0 (0.00%) aligned concordantly >1 times
----
12301 pairs aligned concordantly 0 times; of these:
2127 (17.29%) aligned discordantly 1 time
----
10174 pairs aligned 0 times concordantly or discordantly; of these:
20348 mates make up the pairs; of these:
12318 (60.54%) aligned 0 times
8030 (39.46%) aligned exactly 1 time
0 (0.00%) aligned >1 times
75.36% overall alignment rate
real 0m13.745s
user 0m11.800s
sys 0m1.182s
```
samtools view -b bowtie2-algA1.sam > bowtie2-algA1.bam
samtools sort bowtie2-algA1.bam -o bowtie2-algA1-sorted.bam
```
`samtools index bowtie2-algA1-sorted.bam`
MAPQ = Score de alinhamento = quanto maior melhor.
Importante: checar qual o metodo para calcular p, pois algumas são fixos.
Tophat>>> Star
Tophat> detecção abinitio de sitios de splicing:nao precisa passar nenhuma coordenada de sitios de junção, ele é responsável pela identificação desses sitios.
```
STAR --runThreadN 7 \
--runMode genomeGenerate \
--genomeFastaFiles ./genome.fa \
--genomeDir ./star_index \
--sjdbGTFfile ./genome.gtf \
--sjdbOverhang 149
```
STAR --runThreadN 7 \
--genomeDir ./state/partition1/jcardoso/Ref/star_index \
--readFilesIn ./home/jcardoso/sim/raw/SAMPLEA1_R1.fastq \
./home/jcardoso/sim/raw/SAMPLEA1_R2.fastq \
--sjdbGTFfile ./state/partition1/jcardoso/Ref/genome.gtf \
--outFileNamePrefix ./state/partition1/jcardoso/Ref/star_out/align
STAR --runThreadN 7 --genomeDir ./state/partition1/jcardoso/Ref/star_index --readFilesIn ./home/jcardoso/sim/raw/SAMPLEA1_R1.fastq ./home/jcardoso/sim/raw/SAMPLEA1_R2.fastq --sjdbGTFfile ./state/partition1/jcardoso/Ref/genome.gtf --outFileNamePrefix ./state/partition1/jcardoso/Ref/star_out/align
Sign in with Wallet
Connect another wallet