El propósito de estos ejercicios es obtener la anotación del genoma con el que estamos trabajando, y utilizar esta anotación para catalogar nuestro mapeo de la secuenciación masiva. El resultado de esta práctica será obtener una figura mostrando la frecuencia con la que hemos secuenciado fragmentos de distintos tipos de elementos genómicos (regiones codificantes rRNA, snoRNA, etc). Haremos todo esto dentro del ambiente de R/Bioconductor.
Exactamente el mismo principio puede aplicarse a proyectos sin un genoma completo, siempre y cuando se cuente con algún tipo de archivo describiendo la anotación de distintas regiones de las secuencias de referencia (contigs, ESTs, cDNAs, etc).
Durante esta práctica vamos a utilizar:
| R y BioConductor | Instalados |
| Paquetes de BioC: IRanges, biomaRt, ShortRead, GenomicRanges | Instalados |
En esta práctica manipularemos muchos datos que se nos presentan en forma de rangos (coordenadas de genes, posiciones de secuencias mapeadas, etc). Una de las maneras más flexibles de hacer esto en R es aprovechando el paquete IRanges.
library(Biobase) library(IRanges)
Podemos generar un objecto de clase IRanges si sabemos los inicios y finales de cada rango:
x = IRanges(start=c(11,35,40), end=c(20,50,63)) x
IRanges of length 3
start end width
[1] 11 20 10
[2] 35 50 16
[3] 40 63 24
Para extraer la información básica de estos rangos usamos:
start(x) # Todos los inicios end(x) # Todos los finales width(x) # Ancho de cada rango range(x) # Rango total de todos
Pero existen muchas más funciones de gran utilidad:
coverage(x) # La suma de la cobertura en cada posición
'integer' Rle of length 63 with 6 runs Lengths: 10 10 14 5 11 13 Values : 0 1 0 1 2 1
reduce(x) # Une los rangos que estan encimados
IRanges of length 2
start end width
[1] 11 20 10
[2] 35 63 29
y para darles una idea de lo que pretendemos:
exons = reduce(x) reads = IRanges(start=c(1,21,30,50,80), width=20) reads
IRanges of length 5
start end width
[1] 1 20 20
[2] 21 40 20
[3] 30 49 20
[4] 50 69 20
[5] 80 99 20
countOverlaps(exons, reads)
[1] 1 3
Para mayores detalles, les recomiendo el vignette de IRanges (openVignette(package="IRanges")).
Lo primero que necesitamos es información sobre los elementos que nos interesa analizar de un genoma. Por ahora trabajaremos únicamente con genes conocidos del genoma humano. Para cada uno de estas regiones anotadas del genoma, necesitamos saber sus coordenadas y la cadena en la que se encuentran. También nos conviene saber el tipo biológico de cada transcrito, si representa una proteina, o un RNA no codificante. Vamos a extraer toda esta información de Ensembl, utilizando biomaRt, que es la inferase en R a BioMart.
Por cuestiones de tiempo, veremos un ejemplo muy básico de como usar biomaRt y procederemos a bajar una serie de datos obtenidos previamente.
library(biomaRt)
mart = useMart("ensembl", dataset="hsapiens_gene_ensembl")
# Los atributos que queremos obtener
attribs = c("ensembl_transcript_id", "chromosome_name", "strand", "exon_chrom_start", "exon_chrom_end")
# Obtener la información de la bases de datos
transInfo = getBM(attributes=attribs, filters="ensembl_transcript_id", values="ENST00000337138", mart=mart)
transInfo
ensembl_transcript_id chromosome_name strand exon_chrom_start exon_chrom_end 1 ENST00000337138 14 1 62162258 62162557 2 ENST00000337138 14 1 62187100 62187290 3 ENST00000337138 14 1 62188227 62188372 4 ENST00000337138 14 1 62188457 62188541 5 ENST00000337138 14 1 62193424 62193536 ...
Si quieren ver la lista completa de los atributos que podemos extraer, usen la función listAttributes(mart). Con listFilters(mart) obtenemos una tabla con todos los filtros válidos. Si quieren aprender más, vienen varios ejercicios en el vignette que se instala junto con biomaRt (openVignette(package="biomaRt")).
Para continuar con la práctica, preparé una tabla usando biomaRt, pidiendo una lista más grande de atributos, pero solamente para el cromosoma 14 mediante los parámetros filters="chromosome_name" y values="14".
Bajen entonces la tabla con la anotación del cromosoma 14.
Leamos la tabla que bajaron en R y veamos que contiene:
coords = read.table("human_chr14_annotation.tab.gz", header=TRUE, sep="\t", check.names=FALSE)
head(coords)
ensembl_gene_id ensembl_transcript_id chromosome_name strand exon_chrom_start exon_chrom_end transcript_biotype 1 ENSG00000249549 ENST00000493038 chr14 - 19940995 19941024 lincRNA 2 ENSG00000249549 ENST00000493038 chr14 - 19858667 19861377 lincRNA 3 ENSG00000251363 ENST00000515218 chr14 + 41423916 41424157 lincRNA 4 ENSG00000251363 ENST00000515218 chr14 + 41431994 41432047 lincRNA 5 ENSG00000251363 ENST00000515218 chr14 + 41443545 41443595 lincRNA 6 ENSG00000251363 ENST00000515218 chr14 + 41444546 41445068 lincRNA
Hay una clase de objeto muy util para guardar información acerca de coordenadas genómicas, llamado GRanges. Para usarlo, necesitamos asegurarnos de tener cargado el paquete GenomicRanges.
library(GenomicRanges) exonGRs = GRanges( seqnames = Rle(factor(coords$chromosome_name)), ranges = IRanges(start=coords$exon_chrom_start, end=coords$exon_chrom_end), strand = Rle(factor(coords$strand)), biotype = as.character(coords$transcript_biotype), id = as.character(coords$ensembl_transcript_id) ) exonGRs
GRanges with 24051 ranges and 2 elementMetadata values
seqnames ranges strand | biotype id
<Rle> <IRanges> <Rle> | <character> <character>
[1] chr14 [19940995, 19941024] - | lincRNA ENST00000493038
[2] chr14 [19858667, 19861377] - | lincRNA ENST00000493038
[3] chr14 [41423916, 41424157] + | lincRNA ENST00000515218
... ... ... ... ... ... ...
[24049] chr14 [ 21145178, 21145259] - | tRNA_pseudogene ENST00000461087
[24050] chr14 [ 24027031, 24027100] - | tRNA_pseudogene ENST00000474722
[24051] chr14 [ 74055529, 74055601] - | tRNA_pseudogene ENST00000489070
seqlengths
chr14
NA
Podemos obtener acceso a los componentes de este nuevo objecto:
seqnames(exonGRs) ranges(exonGRs) strand(exonGRs) values(exonGRs)
Lo último corresponde a una tabla donde podemos tener todas las columnas que querramos con información extra.
Podemos manipular los datos sencillamente:
table(values(exonGRs)$biotype)
IG_C_gene IG_C_pseudogene IG_D_gene
57 8 30
IG_J_gene IG_J_pseudogene IG_V_gene
71 3 153
IG_V_pseudogene Mt_tRNA_pseudogene TR_V_gene
110 10 36
TR_V_pseudogene antisense lincRNA
4 108 397
miRNA misc_RNA non_coding
92 46 122
nonsense_mediated_decay polymorphic_pseudogene processed_pseudogene
1132 2 153
processed_transcript protein_coding pseudogene
1600 18815 17
rRNA rRNA_pseudogene retained_intron
10 3 782
retrotransposed scRNA_pseudogene snRNA
4 27 65
snRNA_pseudogene snoRNA snoRNA_pseudogene
4 97 4
tRNA_pseudogene transcribed_processed_pseudogene transcribed_unprocessed_pseudogene
4 5 60
unitary_pseudogene unprocessed_pseudogene
6 14
Una advertencia. Aunque no lo vamos a tomar en cuenta aquí, muchos de los rangos de nuestra tabla podrían estar contenidos unos dentro de otros. Como lo que queremos es contar ocurrencias de nuestras secuencias cortas dentro de estos rangos, podríamos terminar contando algunas de las secuencias más de una vez.
Ejercicios:
Si cerraron su sesión de R después de la práctica anterior, vamos a tener que leer una vez más los resultados del alineamiento con Bowtie. Si no los tienen, pueden bajarlas aquí: clean_reads.bwt.zip.
library(ShortRead)
aln = readAligned("clean_reads.bwt", type="Bowtie")
aln
class: AlignedRead length: 393848 reads; width: 51 cycles chromosome: chr14 chr14 ... chr14 chr14 position: 74190113 53512945 ... 105477651 23536730 strand: + + ... - + alignQuality: NumericQuality alignData varLabels: similar mismatch
Lo primero que tenemos que hacer es convertir este objeto en uno que podamos comparar con nuestras coordenadas de exons:
alnGRs = as(aln, "GRanges") alnGRs
GRanges with 393848 ranges and 3 elementMetadata values
seqnames ranges strand | id similar mismatch
<Rle> <IRanges> <Rle> | <BStringSet> <integer> <character>
[1] chr14 [ 74190113, 74190163] + | GAII06_0003:3:1:1040:19544#0/1 0
[2] chr14 [ 53512945, 53512995] + | GAII06_0003:3:1:1148:2517#0/1 0
[3] chr14 [ 54458302, 54458352] + | GAII06_0003:3:1:1180:16552#0/1 0
... ... ... ... ... ... ... ...
[393846] chr14 [ 38149602, 38149652] - | GAII06_0003:3:120:19796:3235#0/1 0
[393847] chr14 [105477651, 105477701] - | GAII06_0003:3:120:19801:16968#0/1 0
[393848] chr14 [ 23536730, 23536780] + | GAII06_0003:3:120:19805:4487#0/1 0
seqlengths
chr14
NA
Con la función countOverlaps podemos encontrar cuantas secuencias alineadas traslapan con cada uno de nuestros exones. De una vez guardamos el resultado directamente como anotación de los exones:
values(exonGRs)$counts = countOverlaps(exonGRs, alnGRs) values(exonGRs)
DataFrame with 24051 rows and 3 columns
biotype id counts
<character> <character> <integer>
1 lincRNA ENST00000493038 0
2 lincRNA ENST00000493038 4
3 lincRNA ENST00000515218 1
... ... ... ...
24049 tRNA_pseudogene ENST00000461087 0
24050 tRNA_pseudogene ENST00000474722 0
24051 tRNA_pseudogene ENST00000489070 0
Finalmente, podemos agregar estos valores de la manera que más nos convenga. En este caso podría ser por tipo de transcrito, o por transcrito individual.
typeCounts = aggregate(values(exonGRs)$counts, by=list("biotype"=values(exonGRs)$biotype), sum)
head(typeCounts)
biotype x
...
11 antisense 347
12 lincRNA 502
13 miRNA 87
14 misc_RNA 167
15 non_coding 5
16 nonsense_mediated_decay 12548
17 polymorphic_pseudogene 0
18 processed_pseudogene 2627
19 processed_transcript 12472
20 protein_coding 329940
21 pseudogene 469
22 rRNA 32
23 rRNA_pseudogene 1
24 retained_intron 16184
25 retrotransposed 345
26 scRNA_pseudogene 5
27 snRNA 11
28 snRNA_pseudogene 1
29 snoRNA 27
...
transCounts = aggregate(values(exonGRs)$counts, by=list("id"=values(exonGRs)$id), sum)
head(transCounts)
id x
1 ENST00000005198 267
2 ENST00000013070 382
3 ENST00000013125 235
4 ENST00000031146 602
5 ENST00000045347 29
6 ENST00000163416 327
Las cuentas por transcrito las podríamos usar directamente en un proceso de análisis de expresión diferencial, por ejemplo. La distribución de secuencias cortas en los distinto biotipos podríamos querer visualizarlo en una figura:
minCount = 1000
typeCountsHigh = typeCounts[typeCounts$x > minCount,]
typeCountsHigh = typeCountsHigh[order(typeCountsHigh$x),]
typeCountsHigh = rbind(data.frame("biotype"="other",
"x"=sum(typeCounts$x[typeCounts$x <= minCount])),
typeCountsHigh)
pie(typeCountsHigh$x, labels=typeCountsHigh$biotype,
main="Number of aligned reads per biotype")
