2019-09-25
BiostringsBioconductor comprende 1741 paquetes interoperables (Septiembre de 2019) aportados por una gran y diversa comunidad de científicos. Los paquetes cubren una gama de aplicaciones bioinformáticas y estadísticas.
1. Huber, W., et. al. (2015). Orchestrating high-throughput genomic analysis with Bioconductor. Nature methods, 12(2), 115.
library(Biobase)
A partir de la versión 3.5.0 de R, la instalación de Bioconductor cambió. Si tienen una versión anterior, lo siguiente debe bastar:
source("http://bioconductor.org/biocLite.R")
biocLite("Biostrings")
Para las versiones más nuevas de R (de la 3.5.0 en adelante), pueden instalar Bioconductor de la siguiente manera:
if (!requireNamespace("BiocManager"))
install.packages("BiocManager")
BiocManager::install("Biostrings") # Para instalar un paquete en particular
library(Biobase) library(Biostrings)
dnaSeq = DNAString("TTCAGATCTAGTTCGTGTGTGACTGATGATCTGTCACACGTTTTTCTGATCTTCTGACTAGTCGAT")
dnaSeq
## 66-letter "DNAString" instance ## seq: TTCAGATCTAGTTCGTGTGTGACTGATGATCTGTCACACGTTTTTCTGATCTTCTGACTAGTCGAT
class(dnaSeq)
## [1] "DNAString" ## attr(,"package") ## [1] "Biostrings"
Biostrings define una nueva clase de objeto XString, el cual puede ser de tipo BString (cadenas genéricas), DNAString, RNAString o AAString. Veamos como trabajar con un DNAString.
Podemos accesar a parte de la secuencia utilizando [ ]:
dnaSeq[4]
## 1-letter "DNAString" instance ## seq: A
dnaSeq[4:10]
## 7-letter "DNAString" instance ## seq: AGATCTA
length(dnaSeq) reverseComplement(dnaSeq) alphabetFrequency(dnaSeq)
También podemos leer archivos FASTA directamente. Para este ejercicio pueden descargar un archivo que contiene las secuencias de algunos genes de la bacteria E. coli.
fasta_url <- 'http://datos.langebio.cinvestav.mx/~cei/cursos/BP_2018/data/ecoliORFs.fa' dnaSeqs = readDNAStringSet(fasta_url) dnaSeqs
## A DNAStringSet instance of length 11 ## width seq names ## [1] 66 ATGAAACGCATTAGCACCAC...CAGGTAACGGTGCGGGCTGA thrL ## [2] 2463 ATGCGAGTGTTGAAGTTCGG...CATGGAAGTTAGGAGTCTGA thrA ## [3] 933 ATGGTTAAAGTTTATGCCCC...CACGAGTACTGGAAAACTAA thrB ## [4] 1287 ATGAAACTCTACAATCTGAA...TGATGATGAATCATCAGTAA thrC ## [5] 504 ATGCCGGGCAACAGCCCGCA...ACATAAAACACTATCAATAA insB ## ... ... ... ## [7] 606 ATGGCAGAGAAATTTATCAA...AACCTGCGTTTATGAATTAA leuD ## [8] 1401 ATGGCTAAGACGTTATACGA...ACATTCGCAACATTAAATAA leuC ## [9] 1092 ATGTCGAAGAATTACCATAT...ATGTAGCAGAAGGGGTGTAA leuB ## [10] 1572 ATGAGCCAGCAAGTCATTAT...ACAACAAGGAAACCGTGTGA leuA ## [11] 87 ATGACTCACATCGTTCGCTT...TGAGCGGCATCCAGCATTAA leuL
length(dnaSeqs) names(dnaSeqs) width(dnaSeqs) subseq(dnaSeqs,11,30) dinucleotideFrequency(dnaSeqs) reverseComplement(dnaSeqs) translate(dnaSeqs)
subseq(dnaSeqs,11,30)
## A DNAStringSet instance of length 11 ## width seq names ## [1] 20 TTAGCACCACCATTACCACC thrL ## [2] 20 TGAAGTTCGGCGGTACATCA thrA ## [3] 20 TTTATGCCCCGGCTTCCAGT thrB ## [4] 20 ACAATCTGAAAGATCACAAC thrC ## [5] 20 ACAGCCCGCATTATGGGCGT insB ## ... ... ... ## [7] 20 AATTTATCAAACACACAGGC leuD ## [8] 20 CGTTATACGAAAAATTGTTC leuC ## [9] 20 ATTACCATATTGCCGTATTG leuB ## [10] 20 AAGTCATTATTTTCGATACC leuA ## [11] 20 TCGTTCGCTTTATCGGTCTA leuL
dinucleotideFrequency(dnaSeqs)
## AA AC AG AT CA CC CG CT GA GC GG GT TA TC TG TT ## [1,] 3 10 2 5 11 7 3 1 2 4 4 2 4 1 3 3 ## [2,] 176 113 102 161 130 143 199 143 160 220 170 142 86 139 221 157 ## [3,] 52 40 52 50 52 50 80 49 59 94 85 56 31 47 77 58 ## [4,] 103 66 65 70 74 64 105 73 93 119 79 73 34 67 115 86 ## [5,] 36 31 21 33 34 18 45 29 29 50 43 27 22 27 40 18 ## [6,] 21 20 9 14 25 16 18 21 7 25 19 15 12 19 19 15 ## [7,] 50 37 29 33 36 34 45 35 44 55 37 26 19 24 51 50 ## [8,] 131 86 54 74 90 101 133 62 80 133 110 72 44 66 98 66 ## [9,] 69 60 50 73 75 73 100 50 70 119 67 52 38 46 91 58 ## [10,] 138 91 78 103 111 97 138 65 100 136 97 89 61 87 109 71 ## [11,] 3 7 4 6 6 1 8 8 4 7 3 4 7 8 3 7
reverseComplement(dnaSeqs)
## A DNAStringSet instance of length 11 ## width seq names ## [1] 66 TCAGCCCGCACCGTTACCTG...GTGGTGCTAATGCGTTTCAT thrL ## [2] 2463 TCAGACTCCTAACTTCCATG...CCGAACTTCAACACTCGCAT thrA ## [3] 933 TTAGTTTTCCAGTACTCGTG...GGGGCATAAACTTTAACCAT thrB ## [4] 1287 TTACTGATGATTCATCATCA...TTCAGATTGTAGAGTTTCAT thrC ## [5] 504 TTATTGATAGTGTTTTATGT...TGCGGGCTGTTGCCCGGCAT insB ## ... ... ... ## [7] 606 TTAATTCATAAACGCAGGTT...TTGATAAATTTCTCTGCCAT leuD ## [8] 1401 TTATTTAATGTTGCGAATGT...TCGTATAACGTCTTAGCCAT leuC ## [9] 1092 TTACACCCCTTCTGCTACAT...ATATGGTAATTCTTCGACAT leuB ## [10] 1572 TCACACGGTTTCCTTGTTGT...ATAATGACTTGCTGGCTCAT leuA ## [11] 87 TTAATGCTGGATGCCGCTCA...AAGCGAACGATGTGAGTCAT leuL
translate(dnaSeqs)
## A AAStringSet instance of length 11 ## width seq names ## [1] 22 MKRISTTITTTITITTGNGAG* thrL ## [2] 821 MRVLKFGGTSVANAERFLRV...TAAGVFADLLRTLSWKLGV* thrA ## [3] 311 MVKVYAPASSANMSVGFDVL...EGFVHICRLDTAGARVLEN* thrB ## [4] 429 MKLYNLKDHNEQVSFAQAVT...SHNLPADFAALRKLMMNHQ* thrC ## [5] 168 MPGNSPHYGRWPQHDFTSLK...SVELHDKVIGHYLNIKHYQ* insB ## ... ... ... ## [7] 202 MAEKFIKHTGLVVPLDAANV...LQHDDAIAAYEAKQPAFMN* leuD ## [8] 467 MAKTLYEKLFDAHVVYEAEN...AMAAAAAVTGHFADIRNIK* leuC ## [9] 364 MSKNYHIAVLPGDGIGPEVM...AVSTDEMGDIIARYVAEGV* leuB ## [10] 524 MSQQVIIFDTTLRDGEQALQ...VEKELQRKAQHNENNKETV* leuA ## [11] 29 MTHIVRFIGLLLLNASSLRGRRVSGIQH* leuL
Y al final siempre podemos accesar a la secuencia como una cadena de texto normal:
as.character(dnaSeqs[[1]][1:20])
## [1] "ATGAAACGCATTAGCACCAC"
Todos los paquetes de Bioconductor tienen manuales o vignettes que los acompañan. Para accesar la lista de manuales disponibles, usen la función
openVignette('Biostrings')
¿Cuántas secuencias miden más de 1000?
¿Cuál es el %GC de la secuencia más corta? (tip: ?alphabetFrequency o ?letterFrequency)
Obtén el primer codón de cada secuencia. ¿Para qué aminoácido codifican? (tip: prueba la función table para agrupar o tabular tus resultados)
Ahora obtén el último codón de cada secuencia. Intenta hacerlo de más de una forma. ¿Cómo puedes demostrar que tu resultado está bien?
barplot)Nota: todos estos ejercicios se pueden solucionar exclusivamente con las funciones vistas dentro de las prácticas.
sum(width(dnaSeqs) > 1000)
## [1] 5
genes <- names(dnaSeqs) lengths <- width(dnaSeqs) min_seq <- genes[width(dnaSeqs) == min(lengths)] gc <- sum(alphabetFrequency(dnaSeqs[min_seq])[,c(2,3)]) / min(lengths) gc
## [1] 0.5151515
fcodon <- subseq(dnaSeqs,1,3) translate <- translate(fcodon) table(fcodon)
## fcodon ## ATG GTG ## 10 1
table(translate)
## translate ## M V ## 10 1
lcodon <- subseq(dnaSeqs, lengths - 2, lengths) translate <- translate(lcodon) table(lcodon)
## lcodon ## TAA TGA ## 8 3
table(translate)
## translate ## * ## 11
freq <- letterFrequency(dnaSeqs, letters="TGCA", OR=0)
rownames(freq) <- genes
total <- rowSums(freq)
freq_pct <- freq / total
barplot(t(freq_pct), col=topo.colors(4),
legend = TRUE,
args.legend = list(ncol = 1, bty = "n", x = "right"))
.bioc_packages <- c("Biostrings", "DECIPHER", "phangorn", 'ggtree')
.inst <- .bioc_packages %in% installed.packages()
if(any(!.inst)) {
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(.bioc_packages[!.inst], ask = F)
}
# Load packages into session, and print package version
sapply(c(.bioc_packages), require, character.only = TRUE)
## Warning: package 'DECIPHER' was built under R version 3.5.1
## Warning: package 'phangorn' was built under R version 3.5.2
## Warning: package 'ape' was built under R version 3.5.2
## Warning: package 'ggtree' was built under R version 3.5.2
## Biostrings DECIPHER phangorn ggtree ## TRUE TRUE TRUE TRUE
fasta_url <- 'https://raw.githubusercontent.com/RJEGR/July_2018_bioinfo/master/seqs.fasta' seqs = readDNAStringSet(fasta_url) seqs
## A DNAStringSet instance of length 19 ## width seq names ## [1] 957 CACAAAGACATTGGGACACT...TCGCCACACTCCACGGAAGC CYTC007-12|Homo s... ## [2] 1350 CACAAAGATATTGGAACCCT...CTATAGGTTCATTTATCTCT CYTC1709-12|Loxod... ## [3] 690 ATTGGTGGTTTTGGTAATTG...TGGATTCACGTGCTTATTTT GBMIN138350-18|Ca... ## [4] 639 GGTGCTTGAGCAGGGATGGT...TCCTTTACCAACACTTATTC ACLB006-06|Xenopu... ## [5] 657 ACACTATATCTACTATTCGG...TCCTATACCAACACTTATTC ABRMM054-06|Goril... ## ... ... ... ## [15] 658 TACTTTATACTTACTATTTG...TCCTATATCAACACCTATTC GBMA3857-12|Canis... ## [16] 639 GGAGTATGATCCGGAATAGT...TTTTATACCAACACTTATTC ACMC007-04|Aedes ... ## [17] 957 CACAAAGATATTGGAACACT...TCGCTACACTTCACGGAAGC CYTC1028-12|Pan p... ## [18] 639 GGCGCATGAGCTGGAGTCCT...TCCTATATCAACACCTATTC ACLB019-06|Pongo ... ## [19] 671 ATTGGTACCTTGTATTTAAT...CATTTGTTTTGATTTTTTGG GBMNA2144-19|Stro...
GTR+G+I (Generalized time-reversible with Gamma rate variation) neighbor-joining tree
genes <- names(seqs) strsplit(genes[1], "\\|")[[1]][2]
## [1] "Homo sapiens"
species <- sapply(strsplit(genes, "\\|"), `[`, 2) species
## [1] "Homo sapiens" "Loxodonta cyclotis" ## [3] "Caenorhabditis elegans" "Xenopus laevis" ## [5] "Gorilla gorilla" "Equus ferus caballus" ## [7] "Gasterosteus aculeatus" "Felis catus" ## [9] "Drosophila melanogaster" "Sus scrofa" ## [11] "Sus scrofa domesticus" "Danio rerio" ## [13] "Panthera tigris" "Canis lupus" ## [15] "Canis lupus familiaris" "Aedes aegypti" ## [17] "Pan paniscus" "Pongo pygmaeus" ## [19] "Strongyloides fuelleborni"
names(seqs) <- species
alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA, verbose = FALSE)
phangAlign <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phangAlign)
treeNJ <- NJ(dm)
fit = pml(treeNJ, data=phangAlign)
fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
rearrangement = "stochastic", control = pml.control(trace = 0))
tree <- fitGTR$tree p <- ggtree(tree, ladderize=FALSE) + geom_tiplab() + xlim(0,5) p
p + geom_taxalink('Gasterosteus aculeatus', 'Danio rerio',
color='red', linetype = 'dashed')
p + geom_taxalink('Gasterosteus aculeatus', 'Danio rerio',
color='red', linetype = 'dashed') +
geom_hilight(node=33, fill="darkgreen", alpha=0.5, extend = 1) +
geom_hilight(node=23, fill="steelblue", alpha=0.5, extend = 1) +
geom_cladelabel(30, "Hominidae", barsize=2, offset=1,
angle=90, offset.text=0.1,
hjust=0.5, fontsize=3)
# View nodo position
p + geom_text2(aes(subset=!isTip, label=node))