Pré-processamento de
Sequências de Transcritos
Daniel Guariz Pinheiro, PhD.
Laboratório de Genética Molecular e Bioinformática
Departamento de Genética
Faculdade de Medicina de Ribeirão Preto
Universidade de São Paulo
Planejamento
• Preparação da árvore de diretórios
• Obtenção dos dados
– dataset1 (Roche 454 GS FLX)
– dataset2 (Illumina Genome Analyzer)
• Pré-processamento
– Mapeamento
– Montagem “de novo”
Introdução
FERRAMENTAS PARA O PRÉPROCESSAMENTO DE SEQUÊNCIAS DE
TRANSCRITOS
Preparação (aula1)
•
Organização da estrutura de diretórios para a aula 1
raw
dataset1
processed
aligned
bwa
raw
processed
dataset2
bwa
aligned
bowtie
hg19
classes
ref
refGene
contaminants
bwa
hg19
bowtie
indexes
scripts
contaminants
blast
Preparação (aula2)
•
Organização da estrutura de diretórios para a aula 2
raw
processed
dataset1
aligned
bwa
assembled
newbler
raw
processed
bwa
dataset2
aligned
bowtie
assembled
classes
velvet
hg19
ref
refGene
contaminants
bwa
hg19
bowtie
indexes
scripts
contaminants
blast
Comandos úteis
#Criar diretório (mkdir)
mkdir –p /work/CBAB/nomedoaluno/
#Informar diretório atual (pwd)
pwd
#Trocar de diretório (cd)
cd /work/CBAB/nomedoaluno/
#Listar todo o conteúdo do diretório (ls)
ls /work/CBAB/nomenoaluno/*
#Determinar o tipo do arquivo (file)
file undeterminedfiletype.unk
#Criar um atalho (ln)
ln -s /source/file.txt /destiny/linktofile.txt
#Descompactar arquivos no formato .gz
gunzip file.gz
#Descompactar arquivos no formato .tar.gz ou .tgz
tar -zxvf file.tar.gz
#Descompactar arquivos no formato .tar
tar -xvf file.tar
#Descompactar arquivos no formato .tar.bz2 ou .tar.bz
tar -jxvf file.tar.bz2
#Descompactar arquivos no formato .bz2 ou .bz
bunzip2 file.bz2
#Descompactar arquivos no formato .zip
unzip file.zip
#Atribuir permissão de execução (chmod)
chmod a+x script.sh # p/ todos os usuários
#Imprimir as n linhas de um arquivo (head)
head -10 file1.txt # primeiras 10 linhas
#Imprimir todas linhas de arquivo(s) (cat)
cat file1.txt file2.txt
Repositórios públicos
•
•
•
SRA (NCBI Sequence Read Archive): http://www.ncbi.nlm.nih.gov/sra
ENA (EBI European Nucleotide Archive): http://www.ebi.ac.uk/ena/
DRA (DDBJ Sequence Read Archive): http://trace.ddbj.nig.ac.jp/dra/index_e.shtml
[http://trace.ddbj.nig.ac.jp/dra/documentation_e.shtml]
dataset1
• Pool de 2 amostras de culturas de melanócitos de epiderme humana
normal (454 GS FLX)
– SRA
• http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=download_reads
– GEO
• http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM566260
– SRR063336.sra
• NCBI SRA Toolkit (sff-dump)
–
http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=show&f=software&m=software&s=software
sff-dump -A SRR063336.sra
– SRR063336.sff
– Aspera Download (command-line)
• ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.putty -Q -l100m
[email protected]:/sra/srainstant/reads//ByRun/sra/SRR/SRR063/SRR063336/SRR063336.sra .
dataset2
• Amostra de uma linhagem celular de linfoblastos humanos
obtidas de paciente com câncer de maama (Illumina GA)
– http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=download_reads
– ERR022690.sra
• NCBI SRA Toolkit (fastq-dump)
–
http://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?cmd=show&f=software&m=software&s=software
fastq-dump -A ERR022690.sra
– ERR022690.sra_1.fastq
– ERR022690.sra_2.fastq
– Aspera Download (command-line)
• ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.putty -Q -l100m
[email protected]:/sra/srainstant/reads//ByRun/sra/SRR/SRR022/SRR022690/SRR022690.sra .
Manipulação arquivos fastq
• FASTX-Toolkit
– http://hannonlab.cshl.edu/fastx_toolkit/
– Ferramentas executadas via linha de comando para
manipulação das sequências no formato FASTQ. Por
exemplo, produzem estatísticas de qualidade das leituras,
realizam uma poda de qualidade ou de adaptadores, etc.
• Apropriado para dados de Illumina;
• Pode ser utilizado para dados de Roche 454 convertidos para o
formato FASTQ;
– Galaxy (Giardine et al., 2005)
• http://g2.bx.psu.edu/
Checagem de qualidade
fastx_quality_stats -Q 33 -i ERR022690.sra_1.fastq -o ERR022690.sra_1.fastq_qual_stat
fastx_quality_stats -Q 33 -i ERR022690.sra_2.fastq -o ERR022690.sra_2.fastq_qual_stat
The output TEXT file will have the following fields (one row per column):
column = column number (1 to 36 for a 36-cycles read solexa file)
count
= number of bases found in this column.
min
= Lowest quality score value found in this column.
max
= Highest quality score value found in this column.
sum
= Sum of quality score values for this column.
mean
= Mean quality score value for this column.
Q1
= 1st quartile quality score.
med
= Median quality score.
Q3
= 3rd quartile quality score.
IQR
= Inter-Quartile range (Q3-Q1).
lW
= 'Left-Whisker' value (for boxplotting).
rW
= 'Right-Whisker' value (for boxplotting).
A_Count = Count of 'A' nucleotides found in this column.
C_Count = Count of 'C' nucleotides found in this column.
G_Count = Count of 'G' nucleotides found in this column.
T_Count = Count of 'T' nucleotides found in this column.
N_Count = Count of 'N' nucleotides found in this column.
max-count= max. number of bases (in all cycles)
fastq_quality_boxplot_graph.sh -i ERR022690.sra_1.fastq_qual_stat -o ERR022690.sra_1.fastq_qual_stat.png
fastq_quality_boxplot_graph.sh -i ERR022690.sra_2.fastq_qual_stat -o ERR022690.sra_2.fastq_qual_stat.png
FastQC
Ferramenta para análise e controle de qualidade
• http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/
fastqc seqfile1 seqfile2 .. seqfileN
fastqc
[-o output dir] [--(no)extract] [-f fastq|bam|sam]
[-c contaminant file] seqfile1 .. seqfileN
Estatísticas
•
assemblystats
– http://community.g2.bx.psu.edu/tool/
– Métricas em arquivos fasta
•
•
•
•
•
•
•
Min read length
Max read length
Mean read length
Standard deviation of read length
Median read length
N50 read length
PRINSEQ
– http://prinseq.sourceforge.net
– Métricas em arquivos fasta, qual e fastq
– Filtros
•
•
•
•
•
qualidade
poly(A)
Conteúdo de GC
duplicações
FASTQC
– http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/
– Necessita conversão para fastq
•
sff2fastq
–
https://github.com/indraniel/sff2fastq
Colapsar sequências
Objetivo: Evitar duplicação de leituras ocorridas na etapa de PCR. Adequado para detecção de
mutações, mas não em expressão gênica diferencial, o melhor é utilizar um modelo que considere a
ocorrência dessas duplicações.
•
FASTX-Toolkit
fastx_collapser -Q 33 -i ERR022690.sra.merged.fastq -o ERR022690.sra.merged.collapsed.fasta
•
PRINSEQ (requer MUITA memória!)
–
•
http://prinseq.sourceforge.net/
in-house perl scripts (requer preferencialmente “nsort” - http://www.ordinal.com caso contrário demora muito
tempo!). Entrada é FASTQ.
–
–
http://lgmb.fmrp.usp.br/redmine-1.0/projects/bitutils/repository
dedup.sh
Obs.: Nas três abordagens acima, no caso de leituras paired-end é necessário concatenar os pares, colapsar as
sequências e depois separá-las novamente.
•
Goby (http://campagnelab.org/software/goby/)
–
–
–
goby 3g
goby 3g
goby 3g
fasta-to-compact
tally-reads
-i
compact-to-fasta
--paired-end --quality-encoding Sanger -d
-x
PE_1.txt
PE.compact-reads
-o
myfilter
-t fastq
-f
myfilter-keep.filter
-i PE.compact-reads
-o PE.compact-reads
-o PE_p1.txt
-p PE_pair.txt
Podas de qualidade (1) /
Cauda poli-A/T
•
prinseq-lite.pl
A filtragem de cauda poli-A/T pode reduzir o número de falsos positivos nos
alinhamentos.
A eliminação de regiões de baixa qualidade pode reduzir o número de falsos
negativos.
Principais argumentos:
-derep
: opção 1 (idêntico)
-min_len
: tamanho mínimo
-out_format
: opção 2 (FASTA e QUAL)
-trim_tail_right
: tamanho mínimo
-trim_tail_left : tamanho mínimo
-trim_qual_step
: passo para o deslize da janela
-trim_qual_window
: tamanho da janela
-out_good
: arquivo de saída
Podas de qualidade (2)
• fastq_quality_trimmer -Q 33 -t 31 -i ERR022690.sra_2.fastq o ERR022690.sra_2_trim31.fastq
–
–
–
–
–
-Q: quality score (33 Phred/64 Illumina)
-t : quality threshold
-l : minimum length
-i : input
-o: output
• Trim.pl (by Nik Joshi)
– http://wiki.bioinformatics.ucdavis.edu/index.php/Trim.pl
– Ideal para leituras paired-ends
• perl Trim.pl --type 2 --qual-threshold 30 --length-threshold 31 -qual-type 0 --pair1 dataset2/input/ERR022690.sra_1.fastq --pair2
dataset2/input/ERR022690.sra_2.fastq --outpair1
dataset2/input/ERR022690.sra_1_trim20.fastq --outpair2
dataset2/input/ERR022690.sra_2_trim20.fastq --single
dataset2/input/ERR022690.sra_trim.fastq
Manipulação de arquivos SFF
• Arquivos .sff (standard flowgram format)
– Converter sff para fasta/qual
sffinfo -seq INPUTREADS.sff > READS.fasta
sffinfo -qual INPUTREADS.sff > READS.qual
– Pirograma
sffinfo -flow INPUTREADS.sff > READS.flow
– Somente ids
sffinfo -a INPUTREADS.sff > ACCS.txt
– Gerar outro sff (lista)
sfffile -i ACCS.txt -o OUTREADS_ACCS.sff INPUTREADS.sff
– Gerar outro sff (aleatório 10k reads)
sfffile –pickr 10k –o OUTREADS_10k.sff
INPUTREADS.sff
Homopolímeros
6
5
4
A
C
T
G
3
2
1
0
1
T
2
C
3
A
4
G
5
A
6
7
?c
GG
key sequence (TCAG) – Calibragem do sinal
• Detecção entre sinais 1 e 2 = 100%.
• Detecção entre sinais 5 e 6 = 20%.
8
9
- AAAAA ?a
Protocolos especiais 454
• Paired-Ends
– Orientação fwd-fwd
• Multiplex
– Adaptadores MID
• Scripts úteis
– http://www.genome.ou.edu/informatics.html
Biblioteca Padrão x Biblioteca MID
Montagem Transcriptoma
com leituras 454
• Desafios extras
– Cauda poly(A)
– Genes ribossomais
• Bancos de dados pré-montados (est2assembly)
• http://www.ncbi.nlm.nih.gov/
– Genoma mitocondrial
• Bancos de dados pré-montados (est2assembly)
• http://www.ncbi.nlm.nih.gov/
– Elementos repetitivos (e.g. elementos transponíveis)
• http://www.girinst.org/repbase/
– Adaptadores/Primers
• http://vbc.med.monash.edu/cgi-bin/wiki.pl/454_HowTo
seqclean
http://sourceforge.net/projects/seqclean/
seqclean <seqfile> [-v <vecdbs>] [-s <screendbs>] [-r <reportfile>]
[-o <outfasta>] [-n slicesize] [-c {<num_CPUs>|<PVM_nodefile>}]
[-l <minlen>] [-N] [-A] [-L] [-x <min_pid>] [-y <min_vechitlen>]
[-m <e-mail>]
Parameters
<seqfile>: sequence file to be analyzed (multi-FASTA)
-c use the specified number of CPUs on local machine(default 1)
-n number of sequences taken at once in each search slice (default 2000)
-v comma delimited list of sequence files to use for end-trimming of <seqfile> sequences
(usually vector sequences)
-l during cleaning, consider invalid the sequences sorter than <minlen> (default 100)
-s comma delimited list of sequence files to use for screening <seqfile> sequences for
contamination
(mito/ribo or different species contamination)
-r write the cleaning report into file <reportfile> (default: <seqfile>.cln)
-o output the "cleaned" sequences to file <outfasta> (default: <seqfile>.clean)
-x minimum percent identity for an alignemnt with a contaminant (default 96)
-y minimum length of a terminal vector hit to be considered(>11, default 11)
-N disable trimming of ends rich in Ns (undetermined bases)
-M disable trashing of low quality sequences
-A disable trimming of polyA/T tails
-L disable low-complexity screening (dust)
-I do not rebuild the cdb index file
-m send e-mail notifications to <e-mail>
Reproduzir a poda no arquivo .qual
cln2qual <cln_report> <qual_file>
RepeatMasker
• Mascarar elementos repetitivos espalhados no
genoma
– e.g., transposons, retrotransposons, ISs
• http://www.repeatmasker.org/
RepeatMasker -qq -no_is -nolow <FASTA FILE> -lib <REPEATLIBFILE> -species
<SPECIES>
-qq
: mais rápido porém menos sensível
-no_is
: não mascarar inserção de sequências de bactérias (IS)
-nolow : não mascarar sequências de baixa complexidade
-pa
: número de processadores
-lib
: arquivo com as sequências dos elementos repetitivos no formato FASTA
Filtra reads com mais de 70% de mascaramento.
prinseq-lite.pl -fasta input.fasta.masked -qual input.fasta.qual \
-ns_max_p 70 -out_good input.fasta.masked.cleaned ;
Clipping adapters/primers/barcodes
(FASTQ - Illumina)
• Lista Sequências contaminates (Illumina)
– http://bioinfo-core.org/index.php/Sequence_Contaminant_List
cat ERR022690.sra_2.fastq |\
fastx_clipper -Q 33 -l 31 -v -a ACACTCTTTCCCTACACGACGCTCTTCCGATCT |\
fastx_clipper -Q 33 -l 31 -v -a CGGTCTCGGCATTCCTACTGAACCGCTCTTCCGATCT |\
fastx_clipper -Q 33 -l 31 -v -a ATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATC |\
fastx_clipper -Q 33 -l 31 -v -a CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATC |\
fastx_artifacts_filter -Q 33 -v |\
fastq_quality_filter -Q 33 -q 20 -p 50 -v -o ERR022690.sra_2_cleaned.fastq;
fastx_artifacts_filter – remoção de sequências baixa complexidade;
fastx_quality_filter – remoção das sequências que não possuem –p % de bases com
qualidade maior ou igual a –q ;
Reparar a ordem das leituras
http://lgmb.fmrp.usp.br/redmine-1.0/projects/bitutils/repository
pair-ends.pl
[-h/--help] [
-g1 PEExp_1.fastq -g2 PEExp_2.fastq
-i1 PEExp_1_cleaned.fastq -i2 PEExp_2_cleaned.fastq
-o1 PEExp_1_cleaned_paired.fastq -o2 PEExp_2_cleaned_paired.fastq
–os EExp_cleaned_single.fastq
-h
--help
Help
-g1
-g2
--guidefile1
--guidefile2
Guide file 1 (Original fastq p1 file - pre-filtering)
Guide file 2 (Original fastq p2 file - pre-filtering)
-i1
-i2
--inputfile1
--inputfile2
Input file 1 (Filtered fastq p1 file - post-filtering)
Input file 2 (Filtered fastq p2 file - post-filtering)
-o1
-o2
-os
--outputfile1
--outputfile2
--outputfiles
Output file 1
Output file 2
Output file s
Preparação da Entrada para o
Newbler
• Converter o arquivo pré-processado no formato FASTQ para FASTA e QUAL
prinseq-lite.pl -fastq input.fastq -out_format 2 \
> -out_good /tmp/input
Arquivos gerados:
input.fasta
input.qual
Preparação da Entrada (paired-ends)
para o Velvet
• Une os arquivos em pares P1 (forward) e P2 (reverse)
shuffleSequences_fasta.pl P1.fasta P2.fasta input.fasta
shuffleSequences_fastq.pl P1.fastq P2.fastq input.fastq
@2009-05-05:4:1:2:762#0/1
CCGATTTTCCGGAAAAAGGCTAAAACTACAAAGNNN
+2009-05-05:4:1:2:762#0/1
ababa`a`aaaaababaaaa_aaba`W`aabaYDDD
@2009-05-05:4:1:2:762#0/2
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+2009-05-05:4:1:2:762#0/2
DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
@2009-05-05:4:1:2:1736#0/1
ATGCGCATGGCCACCCCGCTGCTGATGCAGGCGNNN
+2009-05-05:4:1:2:1736#0/1
`]`a^aaT`\a`a`a`a^`\][KO\RX`[MM\PDDD
@2009-05-05:4:1:2:1736#0/2
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
+2009-05-05:4:1:2:1736#0/2
DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD
P1
P2
Referências
• Schmieder R, Edwards R. Quality control and preprocessing
of metagenomic datasets. Bioinformatics. 2011 Mar
15;27(6):863-4. Epub 2011 Jan 28. PubMed
PMID:21278185; PubMed Central PMCID: PMC3051327;
• Tarailo-Graovac M, Chen N. Using RepeatMasker to
identify repetitive elements in genomic sequences. Curr
Protoc Bioinformatics. 2009 Mar;Chapter 4:Unit 4.10.
PubMed PMID: 19274634;
• http://sourceforge.net/projects/seqclean/
• http://seqanswers.com/
• http://trace.ddbj.nig.ac.jp/dra/documentation_e.shtml
• http://prinseq.sourceforge.net/manual.html
• http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/
• http://hannonlab.cshl.edu/fastx_toolkit/index.html
Daniel Guariz Pinheiro
[email protected]
Tutorial
• http://lgmb.fmrp.usp.br/~daniel/downloads/TutorialCBAB.docx
• os arquivos ecoli_... rDNA... mito... (têm que ser
descomprimidos de est2assembly_dataC.tar)
• Arquivos já foram baixados !!! Não façam download dos dados
aqui no curso...
– /home/labinfo/DanielGP
• Os arquivos neste diretório devem ser copiados para as respectivas pastas
dentro da estrutura de diretórios organizada (classes/).
• Alternativa: criar links simbólicos (somente usuários avançados);
• Diretório de trabalho
– /cbab/labinfo/
• ALMOÇO: 13:00 RETORNO: 14:00 hrs