• Ejercitándonos con datos genómicos (read.table, rownames, colnames, str)
• Ordenando un tabla completa (sort y order)
• Tamaño y composición de los genomas (Operaciones)
• Reutilizando el código

url <- "http://datos.langebio.cinvestav.mx/~cei/cursos/BP_2018/data/ensembl_info.tab"

rownames(ensTab)

##  [1] "Anole_lizard"    "C_elegans"       "Cat"            
##  [4] "Chicken"         "Chimpanzee"      "Cow"            
##  [7] "Dog"             "Dolphin"         "Duck"           
## [10] "Elephant"        "Fruitfly"        "Fugu"           
## [13] "Gorilla"         "Horse"           "Human"          
## [16] "Microbat"        "Mouse"           "Panda"          
## [19] "Pig"             "Platypus"        "Rat"            
## [22] "Stickleback"     "Tasmanian_devil" "Xenopus"        
## [25] "Zebrafish"

colnames(ensTab)

##  [1] "coding_genes"            "ncrna_genes"            
##  [3] "pseudo_genes"            "coding_gene_avg_length" 
##  [5] "ncrna_gene_avg_length"   "pseudo_gene_avg_length" 
##  [7] "coding_trans_avg_length" "ncrna_trans_avg_length" 
##  [9] "pseudo_trans_avg_length" "genome_length"

str(ensTab)

## 'data.frame':    25 obs. of  10 variables:
##  $ coding_genes           : int  18595 20362 19493 18346 18759 19994 19856 16550 15634 20033 ...
##  $ ncrna_genes            : int  7168 922 1855 6490 8681 3825 11898 3790 567 2644 ...
##  $ pseudo_genes           : int  157 1658 542 43 572 797 950 925 249 568 ...
##  $ coding_gene_avg_length : num  35333 3088 36567 24355 53192 ...
##  $ ncrna_gene_avg_length  : num  10319 334 103 12638 126 ...
##  $ pseudo_gene_avg_length : num  878 1390 792 1013 902 ...
##  $ coding_trans_avg_length: num  2572 1400 2170 2073 2541 ...
##  $ ncrna_trans_avg_length : num  548 256 103 859 126 ...
##  $ pseudo_trans_avg_length: num  869 856 743 994 883 ...
##  $ genome_length          : num  1.80e+09 1.00e+08 2.46e+09 1.23e+09 3.31e+09 ...

ensTab["Human",]

Podemos especificar tanto el nombre del renglón como el nombre de la columna que necesitan:

ensTab["Human","genome_length"]

## [1] 3096649726

Quizá sea más fácil de visualizar si lo representamos como el tamaño pero en millones de bases (normalmente hablamos de megabases o Mb):

round(ensTab["Human","genome_length"]/1e6,1)

## [1] 3096.6

• Les interesa explorar más de cerca el tamaño de los genomas. Generen un nuevo objeto de nombre genomes que tenga solamente los valores de los tamaños de los genomas, expresados en Megabases, y manteniendo el nombre de cada genoma.

• Con la función sort ordenen el resultado para ver fácilmente tanto los valores más pequeños como los más grandes, y sobreescriban el objeto original para preservarlo ordenado.

• Conviertan todos los tamaños de genomas a porcentajes del genoma humano. Es decir, el genoma humano ahora va a medir 100, uno que mida la mitad será 50 y uno que mida tan solo la décima parte medirá 10.

genomes <- round(ensTab[, 'genome_length']/1e6, 1)
names(genomes) <- rownames(ensTab)
genomes <- sort(genomes)
genomes <- genomes/genomes['Human'] * 100

La función sort es muy útil para ordenar vectores (una sola dimensión).

En el caso de tablas existe una función diferente llamada order. Veamos un ejemplo muy sencillo:

x = c("hongo","bacteria","planta","virus")
y = c(100, 10, 1000, 1)
z = data.frame(x,y)
z[order(z$y), ] 

##          x    y
## 4    virus    1
## 2 bacteria   10
## 1    hongo  100
## 3   planta 1000

• Usen la función order para reordenar la tabla ensTab de acuerdo a la longitud de los genomas (y guárdenlo así).

• Si los genomas varían tanto de tamaño, pero el número de genes permanece más o menos constante, ¿qué es lo que ocupa el resto del espacio en los genomas?

Debido a que tenemos una tabla con 25 genomas, la tarea no se vuelve tan trivial. Podriamos elaborar una funcion que indique las instrucciones elaboradas para el genoma humano como a continuación

feature_size <- function(x) {
  gene_size <- x['coding_genes'] * x[ 'coding_gene_avg_length'] 
  trans_size <- x['coding_genes'] * x[ 'coding_trans_avg_length']
  
  intr_size <- gene_size - trans_size
  
  genome_pct <- abs(gene_size - x['genome_length']) / x['genome_length']
  trans_pct <- trans_size / x['genome_length']
  intr_size_pct <- intr_size / x['genome_length']

  
  out <- data.frame(genome_pct, trans_pct, intr_size_pct)
  out <- out * 100
  return(out)
}

Finalmente, aplicamos esta funcion a cada uno de los genomas de la tabla ensTab

apply_features <- apply(ensTab, 1, feature_size)
features <- do.call(rbind, apply_features)
names(features) <- c('Genomic', 'Exonic', 'Intronic')

## 'data.frame':    25 obs. of  3 variables:
##  $ Genomic : num  63.5 37.3 71 63.7 69.9 ...
##  $ Exonic  : num  2.66 28.42 1.72 3.09 1.44 ...
##  $ Intronic: num  33.9 34.3 27.3 33.2 28.7 ...

## Warning: package 'ggplot2' was built under R version 3.5.2

Avanzado: Hacer el gráfico

library(reshape2)

features$Specie <- rownames(features)

mfeat <- melt(features,id.vars = 'Specie', 
              value.name = 'size',variable.name = 'feature', 
              factorsAsStrings = FALSE)

library(ggplot2)
library(RColorBrewer)

getPalette = colorRampPalette(brewer.pal(3, "Paired"))
colourCount <- 3

ggplot(mfeat, aes(x = Specie,y = size, fill = feature)) + 
  geom_bar(stat = "identity", position = "stack", color = "black") + 
  coord_flip() + theme_bw() + labs(y = 'Genome Fraction') +
  scale_fill_manual(values = getPalette(colourCount))