Download Análisis de Datos en R: Aplicaciones Genómicas

Análisis de Datos en R: Aplicaciones Genómicas
Introducción a RNASeq
UNSAM
1,2 y 3 de junio de 2015
UNSAM
Dı́a 2
2/6/2015
1 / 30
UNSAM
Dı́a 2
2/6/2015
2 / 30
RNAseq
UNSAM
Dı́a 2
2/6/2015
3 / 30
Workflow
UNSAM
Dı́a 2
2/6/2015
4 / 30
Diseño Experimental
Cuáles son mis preguntas?
Anotar el transcriptoma?
Si tengo referencia:
genoma bien anotado?
transcritpoma bien anotado?
Quiero cuantificar la expresión?
Nivel Gen?
Nivel Isoforma?
Nivel Exon?
Cuánto puedo secuenciar? ($$$)
solo mRNA?
todo el RNA? (Puedo ver otros RNAs: RNA pequeños, lncRNAs, etc...)
UNSAM
Dı́a 2
2/6/2015
5 / 30
Genoma Referencia
Si existe genoma de referencia, consiste en un archivo fasta que
contene las secuencias.
>chr1
ATCGTAGCTGATGCTGACTGATGCTGACTGACTGACTGACTGACTGCA
CTGACTGACTGACTGACTGCAATCGTAGCTGATGCTGACTGATGCTGA
CTGACTGACTGACTGCAATCGTAGCTGATGCTGACTGATGCTGACTGA
Los genomas son ser publicados usando versiones
En general, son los cromosomas o scaffolds si aúnn no ha sido
cerrados
También podemos tener como referencia los CDS o regiones
codificantes
En caso del transcriptoma, normalmente están en archivos de formato
GTF/GFF
UNSAM
Dı́a 2
2/6/2015
6 / 30
Genoma Referencia
ATENCION!!!
Controlar cuidadosamente los nombres y las posiciones del genoma y
transcriptoma
MUCHAS VECES son diferentes!!! Ej: no coinciden los sı́mbolos, las
versiones genoma-tr, los nombres de los cromosomas, posiciones, etc, etc
UNSAM
Dı́a 2
2/6/2015
7 / 30
Genoma Referencia: Dónde buscarlos?
Ensembl: http://www.ensembl.org
BioMart: http://www.biomart.org
UCSC: https://genome.ucsc.edu
Phytozome: http://phytozome.jgi.doe.gov
Gramene: http://www.gramene.org
Los sitios propios de los organismos en estudio:
TAIR (Arabidopsis)
Drosophila http://flybase.org/
Bioconductor tiene varios genomas y transcriptomas en paquetes, con
formatos fácilmente interconvertibles a objetos para trabajar con R
*.org : nombres de genes y vı́as metabólicas
TxDb.* cordenadas genómicas a diferentes niveles
BSGenome.*: secuencias genómicas
biomaRt, rtracklayer: para interaccionar con las DBs mencionadas
arriba
UNSAM
Dı́a 2
2/6/2015
8 / 30
Diseño Experimental
UNSAM
Dı́a 2
2/6/2015
9 / 30
Cobertura
Es el número promedio de lecturas que alinean o cubren sobre una
referencia.
Coverage (X)= total número de lecturas * longitud de la lectura longitud
estimada target
ejemplo:
Para un transcriptoma 50Mb, si se produce 1M de lecturas de 100 pb, la
cobertura teórica (coverage) será: 1000000 * 100 \5000000 = 20X
Dependerá de la aplicación en la que estemos trabajando
Ensamblado de novo ≥ 50 X
Detección de mutaciones, SNPs, rearreglos, 10-30X
ChipSeq se recomienda aprox 100x.
RNAseq se recomienda 10-30X.
calcular la cobertura teórica es complicado porque la expresión de lso
genes no es uniforme.
puede haber variaciones de hasta 5 ordenes de magnitud
UNSAM
Dı́a 2
2/6/2015
10 / 30
Diseño Experimental: Cobertura
https://genohub.com/next-generation-sequencing-guide/
UNSAM
Dı́a 2
2/6/2015
11 / 30
Cobertura
UNSAM
Dı́a 2
2/6/2015
12 / 30
Réplicas
UNSAM
Dı́a 2
2/6/2015
13 / 30
Diseño Experimental: Réplicas
UNSAM
Dı́a 2
2/6/2015
14 / 30
Diseño Experimental:
UNSAM
Dı́a 2
2/6/2015
15 / 30
Mapeo
UNSAM
Dı́a 2
2/6/2015
16 / 30
Mapeo
UNSAM
Dı́a 2
2/6/2015
17 / 30
Mapeo
TopHat: discovering splice junctions with RNA-Seq.
Cole Trapnell,1, Lior Pachter,2 and Steven L. Salzberg1 (2009) Bioinformatics
UNSAM
Dı́a 2
2/6/2015
18 / 30
Output files
SAM (texto plano)
BAM (binarios)
BED
WIG
Contienen información sobre las posiciones de la referencia donde fueron
alineadas las lecturas
Algunos son intercorvertibles con herramientas como: samtools, bedtools
UNSAM
Dı́a 2
2/6/2015
19 / 30
Output files
En nuestro caso utilizaremos archivos BAM (8) que se encuentran en un
paquete de R que está en Bioconductor. Cada archivo BAM fue obtenido:
Alineando las lecturas (paired-end) al genoma hg19 Homo Sapiens
con TopHat2
Extrayendo sólo las lecturas que alinearon al cromosoma 14
Los datos del experimento original están en ArrayExpress, con el número
de acceso: E-MTAB-1147
El objetivo es analizar el impacto sobre la remodelación del
transcriptoma ante el silenciamiento de una proteı́na del spliceosoma
(hnRNP C)
Fue elegido porque está disponible como paquete de bioconductor
Son 4 réplicas por condición
Condiciones Control (CT) , KnockDown (KD)
Para cargar el paquete:
source("http://bioconductor.org/biocLite.R")
biocLite("RNAseqData.HNRNPC.bam.chr14")
UNSAM
Dı́a 2
2/6/2015
20 / 30
Tablas de conteos
Luego del alineamiento hay que contar cuántas lecturas han caı́do en cada
región de interés y dependerá nuevamente de nuestra pregunta biológica y
de la referencia contra la que hayamos decidido alinear.
Hay 2 modelos de trabajo mas relevantes:
1
métodos que se limitan a contar las lecturas que alinean en regiones
especı́ficas
2
métodos que intentan reconstruir las isoformas
UNSAM
Dı́a 2
2/6/2015
21 / 30
Tablas de conteos
Los archivos de alineamiento tienen información sobre las posiciones
en el genoma donde fueron alineadas las lecturas.
Si superponemos cada una de estas cordenadas con las cordenads
genómicas que nos interesan (ej. genes), podemos saber cuántas
lecturas cayeron en cada una de los genes.
Repitiendo este proceso para cada uno de los genes, tendramos una
tabla con todos los genes y la cantidad de lecturas para cada uno:
UNSAM
Dı́a 2
2/6/2015
22 / 30
HTSEQ
Es una herramienta muy popular y veloz para generar las tablas de counts.
UNSAM
Dı́a 2
2/6/2015
23 / 30
Bioconductor
En nuestro caso, usaremos tablas generadas con los paquetes:
GenomicRanges,GenomicFeatures,GenomicAlignments
Con los siguientes comandos:
bamfiles <- RNAseqData.HNRNPC.bam.chr14_BAMFILES
targets <- data.frame(row.names=names(RNAseqData.HNRNPC.bam.chr14_BAMFILES),
bam=matrix(unlist(RNAseqData.HNRNPC.bam.chr14_BAMFILES), ncol = 1, byrow = TRUE),
condition=c(rep("CT",4),rep("KD",4)))
files=as.character(targets$bam)
reads<-lapply(files,function(x) { aln <- readGAlignments(x)}); names(reads)<-rownames(targets)
c14<-loadDb("chr14.sqlite")
feature <-exonsBy(c14, by="gen")
hits<-lapply(reads, function(x){countOverlaps(feature, x, ignore.strand = TRUE) })
hits.ul<-do.call(cbind.data.frame, hits)
genes.hits<-write.table(file="genes.hits.tab")
Análogo al modo de union de HTSEQ
No se han tenido en cuenta que los reads son paired-end
El archivo gene.hits.tab está en la página de práctica
UNSAM
Dı́a 2
2/6/2015
24 / 30
Tablas de conteos: Modo reconstrucción
Cufflinks / Cudffdiff: Sin descubrimiento
UNSAM
Dı́a 2
2/6/2015
25 / 30
Tablas de conteos: Modo reconstrucción
Cufflinks / Cudffdiff: Con descubrimiento
UNSAM
Dı́a 2
2/6/2015
26 / 30
Tablas de conteos: Tuxedo
UNSAM
Dı́a 2
2/6/2015
27 / 30
Tablas de conteos: Tuxedo
RPKM/FPKM = 109 ∗ C /(N ∗ L)
UNSAM
Dı́a 2
2/6/2015
28 / 30
Tablas de conteos: Tuxedo
Para generar la carpeta output con la que trabajaremos con el paquete
cummeRbund, Cuffdiff fue ejecutado en modo ”No descubrimiento” con
los siguientes comandos:
cuffdiff chr14.gtf [archivos.bam] -L CT,KD -o output
RPKM/FPKM
= 109 ∗ C /(N ∗ L)
C: Número de lecturas que caen en los exones
N: total de lecturas de la biblioteca
L: Longitud del gen
UNSAM
Dı́a 2
2/6/2015
29 / 30
Cufflinks
UNSAM
Dı́a 2
2/6/2015
30 / 30