AnalysisSeq-2
Cardenas-Rivera, A.
2023-04-20
Rho - continuación
library(ggplot2)
library(seqinr)
library(msa)## Loading required package: Biostrings## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## Loading required package: S4Vectors## Loading required package: stats4##
## Attaching package: 'S4Vectors'## The following objects are masked from 'package:base':
##
## expand.grid, I, unname## Loading required package: IRanges## Loading required package: XVector## Loading required package: GenomeInfoDb##
## Attaching package: 'Biostrings'## The following object is masked from 'package:seqinr':
##
## translate## The following object is masked from 'package:base':
##
## strsplitseqinr::choosebank("genbank")
sars2 <- seqinr::query("sars2", "AC=MN908947")
covid19 <- unlist(seqinr::getSequence(sars2))
rho1 <- seqinr::rho(covid19) # guardamos el resultado de rho en un objeto
rho1 <- rho1 - 1 # le restamos 1
rho1 <- as.data.frame(rho1) # convertimos el objeto en data frame para que este en dos columnas
ggplot(data = rho1, aes(fill=Var1, y=Freq, x=Var1)) +
geom_bar(position="stack", stat="identity") +
labs(title = "FRECUENCIA RELATIVA DE DINUCLEOTIDOS",
subtitle = "SARS-CoV-2 genoma",
x = "Dinucleotido",
y = "Frecuencia relativa",
caption = Sys.Date()) +
theme(legend.position = "none")
Como podemos ver, es muy evidente que el dinucleotido “cg” esta subrepresentado, y eso significa que la frecuencia en la que ocurre este dinucleotido es mucho menos que la mitad de lo que se espera. Esto sugiere alguna presión evolutiva que evita la formación de este dinucleótido. Por otro lado, el dinucleótido “tg” está sobrerepresentado.
Ahora vamos a analizar la presencia de los codones. Como sabemos, el
código genético se rige por una codificación de tres bases nitrogenadas
= un aminåcido. Para ello vamos a modificar el parámetro
wordsize=3, para que haga el análisis en trinucleótidos y
no en dinucleótidos.
rho2 <- seqinr::rho(covid19, wordsize=3)
rho2 <- rho2 - 1 # le restamos 1
rho2 <- as.data.frame(rho2) # convertimos el objeto en data frame para que este en dos columnas
nrow(rho2)## [1] 64Como la cantidad de filas es grande y lo que nos intersa saber es que
trinucleotidos están sobre- e subrrepresentados, entonces vamos
seleccionar solo los seis más representados y los seis más
subrepresentados usando las funciones head() y
tail() respectivamente.
rho2 <- rho2[order(rho2$Freq),] # primero lo ordenamos de menor a mayor
rho2_ext <- rbind(head(rho2), tail(rho2))
nrow(rho2_ext)## [1] 12Ahora si podemos usar ggplot para graficar y ver si encontramos algún dato interesante
ggplot(data = rho2_ext, aes(fill=Var1, y=Freq, x=Var1)) +
geom_bar(position="stack", stat="identity") +
labs(title = "FRECUENCIA RELATIVA DE DINUCLEOTIDOS",
subtitle = "SARS-CoV-2 genoma",
x = "Dinucleotido",
y = "Frecuencia relativa") +
theme(legend.position = "none")
De forma interesante el codon “CGA” que es el menos representado en el genoma del SARS-CoV-2 también es el menos representado en el genoma del ser Humano, y se ha reportado que este codon cumple un rol particular en el proceso de transcripción del Virus de la leucemia murina. Por otro lado, el trinucleotido mas representado en el genoma es el “ACA”; curiosamente un estudio del 2011 mostro que la expresión de una endonucleasa espercifica de la secuencia ACA inhibe la infeccion de VIH en un modelo experimental in vitro. Este podría ser un tema interesante para estudiar como una alternativa de tratamiento en SARS-CoV-2 (tal vez para la siguiente pandemia).
GC content
Como vimos en la gráfica anterior, el dinucleótido menos representado
por mucho es el “GC”. Curiosamente, los sitios ricos en GC son reguladores -Cis
positivos de la transcripción y están asociados a regiones de inicio
de transcripción, promotores y enhancers. Para saber rapidamente el
contenido (frecuencia relativa) de GC en una secuencia podemos usar la
función GC() del paquete seqinr.
seqinr::GC(covid19)## [1] 0.3797278Para analizar la presencia de GC a lo largo de la secuencia, podemos definir una ventana e ir haciendo el calculo sub-secuencia por subsecuencia, corriendo uno a la vez. Para eso vamos a tener que definir la amplitud de la ventana y hacer uso de la loop for para ir haciendo las iteraciones
win.size <- 100 # definimos el tamaño de la ventana
pos.end <- length(covid19) - win.size # definimos la última posición posible para la ventana
res <- vector() # creamos el vector resultante
for (i in 1:pos.end) { # iniciamos el loop for
j <- i + win.size # definimos la posicion final de la ventana
p1 <- seqinr::GC(covid19[i:j]) # calculamos la concentracion de GC en la subsecuencia especifica
res <- c(res, p1) # agregamos el valor al vector resultante
}
plot(res, type = "l") # generamos un grafico de línea
Con esta forma de analizar tenemos dos limitaciones. La primera es que si desemos escoger diferentes amplitudes de ventana, vamos a tener que repetir las lineas de codigo y la segunda es que necesitamos probar diferentes dimensiones, porque hasta este punto no sabemos si las regiones regulatorias son de 20, 50, 100, 200 o 300 nucleotidos de longitud; lo único que sabemos es que la frecuencia relativa de GC en todo el genoma es de 0.3797 ~ 38%.
Para resolver el primer problema definiremos una función
personalizada llamada GCcontet, que tendrá dos parámetros
bd1 que será la secuencia a analizar y win.len
que será la amplitud de la ventana. Usaremos la misma estructura que
usamos en el ejercicio anterior pero ahora lo metermos dentro de una
función.
GCcontet <- function(db1, win.len){ # inicializamos la función con los dos parámetros
pos.end <- length(db1) - win.len # definimos la última posición posible
res <- vector() # creamos un vector vacío
for (i in 1:pos.end) { # iniciamos el loop
j <- i + win.len # definimos la ultima posicion de la ventana
p1 <- seqinr::GC(db1[i:j]) # calculamos el contenido de GC
res <- c(res, p1) # Agregamos el resultado a la cola
}
return(res) # OUTPUT: devolvemos el resultado
}Con esta función ya podemos calcular el contenido de GC a lo largo de la secuencia con una amplitud de ventana variable
ven20 <- GCcontet(db1 = covid19, win.len = 20)
plot(ven20,type = "l", col = "green")
ven10000 <- GCcontet(db1 = covid19, win.len = 10000)
lines(ven10000, col = "red")
El problema que se genero ahora es que no podemos graficar las dos lineas de forma adecuada porque los vectores son de longitudes diferentes y las “posiciones” de cada dato realmente no corresponden a la posición en la secuencia. Por otro lado, si es que quisieramos usar el paquete ggplot para graficarlo, sería necesario convertirlo a data.frame. Entonces vamos a cambiar un poco la función para que el OUTPUT sea un data.frame y que una columna corresponda a la posición (el centro de la ventana) y la otra columan sea la frecuencia. Así podremos graficarlo usando cordenadas “X” y “Y”.
GCcontet <- function(db1, win.len){
pos.end <- length(db1) - win.len # ultima posicion posible de la ventana
res <- vector() # vector para la frecuencia
pos <- vector() # vector para la posicion
for (i in 1:pos.end) { # iniciamos el loop for
j <- i + win.len # definimos la ultima posicion dentro del loop
p1 <- seqinr::GC(db1[i:j]) # calculo de la frecuencia de CG
res <- c(res, p1) # agregamos el valor a la cola
p2 <- mean(c(i,j)) # calculamos la mitad de la ventana
pos <- c(pos,p2) # agregamos la posicion a la cola
}
res <- data.frame(pos = pos, freq = res) # creamos el data.frame
return(res) # devolvemos el OUTPUT
}Ahora con esta función ya deberíamos poder graficarlo sin problemas
ven20 <- GCcontet(db1 = covid19, win.len = 20)
ven10000 <- GCcontet(db1 = covid19, win.len = 10000)
plot(x=ven20$pos, y=ven20$freq,type="l", col="green")
lines(x=ven10000$pos, y=ven10000$freq, col="red")
Ahora aprovecharemos la versatilidad de las listas para generar el analisis con diferentes amplitudes de listas y poderlos guardar en un mismo objeto pese a que tienen diferentes dimensiones y son data.frames
mainL <- list() # creamos una lista vacia
ventanas <- c(10,20,50,100,200,300,400,500,1000,2000) # definimos el tamaño de las ventanas
for (i in 1:length(ventanas)) { # iniciamos el loop
mainL[[i]] <- GCcontet(db1 = covid19, # en cada slot guardamos data data.frame
win.len = ventanas[i])
}Verifiquemos que son de diferente dimension
summary(mainL)## Length Class Mode
## [1,] 2 data.frame list
## [2,] 2 data.frame list
## [3,] 2 data.frame list
## [4,] 2 data.frame list
## [5,] 2 data.frame list
## [6,] 2 data.frame list
## [7,] 2 data.frame list
## [8,] 2 data.frame list
## [9,] 2 data.frame list
## [10,] 2 data.frame listsapply(mainL, nrow)## [1] 29893 29883 29853 29803 29703 29603 29503 29403 28903 27903Antes de ver las graficas de lineas de los resultados, veamos si la distribucion “normal” de los resultados. Para ello usaremos histogramas, pero nos convendria hacerlo uno encima del otro (overplot). Para ello, nuevamente haremos uso de una lista y un loop para llenar en cada slot de la lista el resultado del histograma
mainHist <- list()
for(i in 1:length(mainL)){
tmp1 <- mainL[[i]]
mainHist[[i]] <- hist(tmp1[,2])
}
Antes de hacer los histogramas, vamos a ver como podemos definir los colores de nuestras graficas. Esto va a ser de suma utilidad cuando publiquemos nuestros articulos, porque podremos: definir los colores que querramos usar, utilizar los mismo colores para cada grafica, y utilizar un codigo de colores como identificador de cada grupo de datos. Para ello, aunque pueden definir los colores que Uds gusten, yo les recomiendo usen paletas de colores ya diseñadas por personas que se dedican a esto. Como la que pueden encontrar en las paletas de adobe
main.col <- c(rgb(196/256,51/256,2/256,0.2),
rgb(237/256,170/256,37/256,0.2),
rgb(183/256,191/256,153/256,0.2),
rgb(10/256,115/256,115/256,0.2),
rgb(1/256,2/256,33/256,0.2))Ahora ya podemos hacer nuestros histogramas, usando nuevamente un loop
plot( mainHist[[1]], col=main.col[1],ylim=c(0,12), freq=F) # first histogram
for(i in 2:5){
plot( mainHist[[i]], col=main.col[i], freq=F, add=T) # second
}
Ahora necesitamos crear un data frame con 3 columnas, una con los valores de “X” (localización), otra con los valores de “Y” (valor de CG) y una tercera columna para el grupo de datos al que pertenece cada dato
res <- data.frame(pos=vector(),freq=vector(),group=vector()) # creamos el data.frame vacio
for(i in 1:length(mainL)){
tmp1 <- cbind(mainL[[i]],
group=rep(paste("window-",
ventanas[i],
sep = ""
),
nrow(mainL[[i]])
)
)
res <- rbind(res,tmp1)
}Ahora que ya tenemos el data.frame en el formato que necesitamos,
podemos usar la funcion geom_lineplot() del paquete
ggplot2
ggplot(data = res, aes(x= pos, y=freq, color = group)) +
geom_line() +
geom_hline(yintercept = seqinr::GC(covid19), color = "red", size = 1.5) +
geom_hline(yintercept = 2*seqinr::GC(covid19), color = "red", size = 1, linetype="dashed")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Como podemos ver, con la ventana de 10 bases tenemos picos cuya
frecuencia es superior al doble de la frecuencia esperada para toda la
secuencia. Para encontrar esas posiciones, podemo usar la funcion
order para poner en orden creciente el data.frame de la
ventana de 10 y luego seleccionar las ultimas filas con la funcion
tail()
ven10 <- mainL[[1]] # extraemos los datos correspondeintes a la ventana 10
ven10 <- ven10[order(ven10$freq),] # ordenamos el data.frame en orden creciente segun la freq
peaks <- tail(ven10,16) # extraemos los ultimos 16 valores
peaks # imprimimos en pantalla el resultado## pos freq
## 28390 28395 0.8181818
## 28391 28396 0.8181818
## 28394 28399 0.8181818
## 28800 28805 0.8181818
## 29194 29199 0.8181818
## 29196 29201 0.8181818
## 29732 29737 0.8181818
## 29733 29738 0.8181818
## 29734 29739 0.8181818
## 29735 29740 0.8181818
## 29736 29741 0.8181818
## 29737 29742 0.8181818
## 15955 15960 0.9090909
## 23603 23608 0.9090909
## 23606 23611 0.9090909
## 29195 29200 0.9090909Para ver las secuencias correspondientes, usamos el valor de la columna “posicion” que sabemos define posición de la mitad de la ventana. Lo que sería interesante es alinear estas secuencias y ver si hay una región concenso que puede sugerir alguna secuencia regulatoria o un blanco terapeutico.
seq <- vector()
tam <- 5
for(i in 1:nrow(peaks)){
p1 <- seqinr::c2s(covid19[(peaks$pos[i]-tam):(peaks$pos[i]+tam)])
seq <- rbind(seq,p1)
}
# generamos el alineamiento
msa(unlist(as.list(seq)), type="dna")## use default substitution matrix## CLUSTAL 2.1
##
## Call:
## msa(unlist(as.list(seq)), type = "dna")
##
## MsaDNAMultipleAlignment with 16 rows and 18 columns
## aln
## [1] ---CCTCGGCGGGC----
## [2] ------CGGCGGGCACG-
## [3] ---GCAGAGGCGGC----
## [4] -----GGGGCCGGCTG--
## [5] --CCGAGGCCACG-----
## [6] ---CGAGGCCACGC----
## [7] ----GAGGCCACGCG---
## [8] ------GGCCACGCGGA-
## [9] -------GCCACGCGGAG
## [10] -----AGGCCACGCGG--
## [11] ACGTCGGCCCC-------
## [12] -CGTCGGCCCCA------
## [13] ------CCCCCAGCGCT-
## [14] -----GCCCCCAGCGC--
## [15] ----TGCCCCCAGCG---
## [16] ----CGGCCCCAAGG---
## Con ----??GGCCC?GCG---Interesantemente tenemos una secuencia “concenso” de 5 bases nitrogenadas, que como ya vimos anteriormente es altamente improbable que sucedan de forma aleatoria. La pregunta es si está relacionado con alguna zona “regulatoria cis”, como secuencia promotora o un enhancer. El articulo de Chen et al (1999) muestra que esta secuencia llamada “Elementos de respuesta a metales de alta afinidad”y que se une a los dedos de Zinc del factor de transcripción MTF-1. Ademas esta secuencia tambien se encuentra presente en la region promotora del factor de necrosis tumoral alfa y su presencia se encuentra disminuida en los pacientes con enfermedad de Takayasu (Naqiang et al 2011).
Si ampliamos un poco mas la ventana de las secuecias a 10 nucleotidos por lado, tenemos tambien un resultado interesante
seq <- vector()
tam <- 10
for(i in 1:nrow(peaks)){
p1 <- seqinr::c2s(covid19[(peaks$pos[i]-tam):(peaks$pos[i]+tam)])
seq <- rbind(seq,p1)
}
# generamos el alineamiento
msa(unlist(as.list(seq)), type="dna")## use default substitution matrix## CLUSTAL 2.1
##
## Call:
## msa(unlist(as.list(seq)), type = "dna")
##
## MsaDNAMultipleAlignment with 16 rows and 29 columns
## aln
## [1] --AAACAACGTCGGCCCCAAGGT------
## [2] ---AACAACGTCGGCCCCAAGGTT-----
## [3] ------AACGTCGGCCCCAAGGTTTAC--
## [4] ------CAATTTGCCCCCAGCGCTTCA--
## [5] -------AATTTGCCCCCAGCGCTTCAG-
## [6] --------ATTTGCCCCCAGCGCTTCAGC
## [7] -------ACCGAGGCCACGCGGAGTACG-
## [8] --------CCGAGGCCACGCGGAGTACGA
## [9] ------CACCGAGGCCACGCGGAGTAC--
## [10] -----TCACCGAGGCCACGCGGAGTA---
## [11] ---TTTCACCGAGGCCACGCGGAG-----
## [12] ----TTCACCGAGGCCACGCGGAGT----
## [13] ATTCTCCTCGGCGGGCACGTA--------
## [14] ---CTCCTCGGCGGGCACGTAGTG-----
## [15] ---AGGGAGCAGAGGC-GGCAGTCA----
## [16] ------TCCTAGGGGCCGGCTGTTTTG--
## Con ------?AC?G?GGCCACGCGG??T?---Esta secuencia concenso “GGCCACGCGG” de 10 nt, mostró en una simulación in silico que puede unirse a los farmacos como Amonafie o Azonafide (Bear & Remers 1996). Y estos fármacos estan en fase III de experimentación como tratamiento contra Leucemia mieloide aguda (Allen & Lundberg 2011). Lo que abre la posiblidad con un blanco terapeutico en contra del COVID.
ANALISIS DE PROTEINAS
El paquete seqinr nos ofrece una función muy útil
para hace una analisis descriptivo de las secuencias de aminoacidos de
una proteína, como por ejemplo: el punto isoelectrico, el peso
molecular, la frecuencia absoluta de los aminoacidos, asi com la
proporcion de aminoacidos basico, acidos, apolares, pequeños, entre
otros. Esta función se llama AAstat(), y el único parámetro
que necesita es meter la secuencia como una cadena de caracteres. Para
este ejemplo vamos a trabajar con tres secuencias:
- Amiloide beta (uniprot=P05067) relacionado con la enfermedad de Alzheimer (AD)
- Actina (uniprot= P60709)relacionado con el citoesqueleto y los microfilamentos
- TREM2 (uniprot=Q9NZC2) un receptor de membrana expresado en microglias en el cerebro y relacionado con AD
Primero para llamar las secuencias usaremos la función
query() del paquete seqinr. Luego
seleccionaremos el base de datos de Uniprot y llamaremos a cada
secuencia usando su ID.
library(seqinr)
seqinr::choosebank(infobank = T)## bank status
## 1 genbank on
## 2 embl on
## 3 emblwgs on
## 4 genbankseqinr on
## 5 swissprot on
## 6 ensembl on
## 7 hogenom7dna on
## 8 hogenom7 on
## 9 hogenom on
## 10 hogenomdna on
## 11 hovergendna on
## 12 hovergen on
## 13 hogenom5 on
## 14 hogenom5dna on
## 15 hogenom4 on
## 16 hogenom4dna on
## 17 homolens on
## 18 homolensdna on
## 19 hobacnucl on
## 20 hobacprot on
## 21 phever2 on
## 22 phever2dna on
## 23 refseq on
## 24 refseq16s on
## 25 greviews on
## 26 bacterial off
## 27 archaeal on
## 28 protozoan on
## 29 ensprotists on
## 30 ensfungi on
## 31 ensmetazoa on
## 32 ensplants on
## 33 ensemblbacteria on
## 34 mito on
## 35 polymorphix on
## 36 emglib on
## 37 refseqViruses on
## 38 ribodb on
## 39 taxodb on
## info
## 1 GenBank Release 246 (15 October 2021) Last Updated: Nov 19, 2021
## 2 EMBL Nucleotide Archive Release 143 (March 2020) Last Updated: Nov 21, 2021
## 3 EMBL Whole Genome Shotgun sequences (July 2018)
## 4 GenBank Release 231 (15 April 2019) Last Updated: Jun 8, 2019
## 5 UniProt Knowledgebase Release 2021_03 of 09-Jun-2021, Last Updated: Aug 6, 2021
## 6 Ensembl 85 - (10/03/16) Last Updated: Oct 3, 2016
## 7 HOGENOM - genomic data - Release 07 (Nov 3,2015) Last Updated: Apr 19, 2017
## 8 HOGENOM - protein data - Release 07 (Nov 3,2015) Last Updated: Jun 3, 2019
## 9 HOGENOM - protein data - Release 06 (Oct 30,2011) Last Updated: May 10, 2012
## 10 HOGENOM - genomic data - Release 06 (Oct 30,2011) Last Updated: Nov 14, 2011
## 11 HOVERGEN - genomic data - Release 49 (Dec 22 2009) Last Updated: Dec 22, 2009
## 12 HOVERGEN - protein data - Release 49 (Dec 22 2009) Last Updated: Dec 22, 2009
## 13 HOGENOM5 (AA)
## 14 HOGENOM5 (DNA)
## 15 HOGENOM4 (AA)
## 16 HOGENOM4 (DNA)
## 17 HOMOLENS 5 - Homologous genes from Ensembl(60)\t Last Updated: Feb 17, 2011
## 18 HOMOLENS 5 - Homologous genes from Ensembl(60)\t Last Updated: Feb 17, 2011
## 19 HOBACGEN - genomic data - Release 10 (February 12 2002)
## 20 HOBACGEN - protein data - Release 10 (February 12 2002)
## 21 PhEVER - protein data - Release 2 (June 1 2010) Last Updated: Jul 22, 2010
## 22 PhEVER - genomic data - Release 2 (June 1 2010) Last Updated: Jul 22, 2010
## 23 Refseq RNA Sequences. Last Updated: Jan 27, 2018
## 24 Refseq RNA 16S 23S Sequences. Last Updated: Jan 29, 2018
## 25 Genome Review from EBI Last Updated: Jan 24, 2013
## 26 NCBI Bacterial Genomes
## 27 Archaeal Genomes from NCBI Last Updated: Aug 25, 2015
## 28 Protozoan Genomes from NCBI Last Updated: Feb 23, 2011
## 29 Ensembl protists 86 - (10/05/16) Last Updated: Oct 5, 2016
## 30 Ensembl fungi 86 - (10/09/16) Last Updated: Oct 9, 2016
## 31 Ensembl protists 86 - (10/05/16) Last Updated: Oct 5, 2016
## 32 Ensembl protists 86 - (10/16/16) Last Updated: Oct 16, 2016
## 33 Ensembl Bacterial Genomes 21 - (02/23/14) Last Updated: Feb 27, 2014
## 34 Mitochondrial sequences - Release 41 (May 19, 2010) Last Updated: Jul 9, 2010
## 35 POLYBASE - Release 1 (June 20, 2003)
## 36 EMGLib Release 5 (December 9, 2003)
## 37 RNA sequences - numrel1 (daterel1) Last Updated: May 10, 2012
## 38 RiboDB
## 39 taxonomic databaseseqinr::choosebank("swissprot")
app <- seqinr::query("app", "AC=P05067")
actin <- seqinr::query("actin", "AC=P60709")
trem2 <- seqinr::query("trem2", "AC=Q9NZC2")Luego extraeremos las secuencias usando la función
getSequence() del mismo paquete (recuerda que las tienes
que sacar de la lista)
app.seq <- unlist(seqinr::getSequence(app))
actin.seq <- unlist(seqinr::getSequence(actin))
trem2.seq <- unlist(seqinr::getSequence(trem2))Finalmente, usaremos la función AAstat() para obtener la
estadistica descriptiva de cada una de ellas
seqinr::AAstat(app.seq)
## $Compo
##
## * A C D E F G H I K L M N P Q R S T V W Y
## 0 63 18 50 92 21 38 25 24 41 56 24 31 35 36 37 35 50 65 9 20
##
## $Prop
## $Prop$Tiny
## [1] 0.2649351
##
## $Prop$Small
## [1] 0.5
##
## $Prop$Aliphatic
## [1] 0.1883117
##
## $Prop$Aromatic
## [1] 0.0974026
##
## $Prop$Non.polar
## [1] 0.4844156
##
## $Prop$Polar
## [1] 0.5155844
##
## $Prop$Charged
## [1] 0.3181818
##
## $Prop$Basic
## [1] 0.1337662
##
## $Prop$Acidic
## [1] 0.1844156
##
##
## $Pi
## [1] 4.730445seqinr::AAstat(actin.seq)
## $Compo
##
## * A C D E F G H I K L M N P Q R S T V W Y
## 0 29 6 23 26 13 28 9 28 19 27 17 9 19 12 18 25 26 22 4 15
##
## $Prop
## $Prop$Tiny
## [1] 0.304
##
## $Prop$Small
## [1] 0.4986667
##
## $Prop$Aliphatic
## [1] 0.2053333
##
## $Prop$Aromatic
## [1] 0.1093333
##
## $Prop$Non.polar
## [1] 0.5546667
##
## $Prop$Polar
## [1] 0.4453333
##
## $Prop$Charged
## [1] 0.2533333
##
## $Prop$Basic
## [1] 0.1226667
##
## $Prop$Acidic
## [1] 0.1306667
##
##
## $Pi
## [1] 5.290482seqinr::AAstat(trem2.seq)
## $Compo
##
## * A C D E F G H I K L M N P Q R S T V W Y
## 0 16 6 13 12 7 19 11 9 6 36 2 4 14 11 12 17 15 10 7 3
##
## $Prop
## $Prop$Tiny
## [1] 0.3173913
##
## $Prop$Small
## [1] 0.4956522
##
## $Prop$Aliphatic
## [1] 0.2391304
##
## $Prop$Aromatic
## [1] 0.1217391
##
## $Prop$Non.polar
## [1] 0.5608696
##
## $Prop$Polar
## [1] 0.4391304
##
## $Prop$Charged
## [1] 0.2347826
##
## $Prop$Basic
## [1] 0.126087
##
## $Prop$Acidic
## [1] 0.1086957
##
##
## $Pi
## [1] 5.838939