Fig_14

pdf("fig_14.pdf", width = 5, height = 5)
par(bty = "n", las = 0, mar = c(4, 4, 2, 1) + 0.2)

# función suma armónicos
sum_arm <- function(n) {
  a <- numeric(n)
  for (i in 1:n) {
    for (j in 1:i) {
      a[j] <- a[j] + 1/i
    }
  }
  a/n
}

# la secuencia de un archivo txt
sec <- readChar("secuencia.txt", file.info("secuencia.txt")$size)
longitud <- nchar(sec)

# posiciones de las metioninas
mets <- unlist(gregexpr("M", sec))
# distancias entre ellas
n <- length(mets) - 1
distancias <- numeric(n)
for (i in 1:n) {
  distancias[i] <- mets[i + 1] - mets[i]
}
# ordenar de mayor a menor distancia
distancias <- sort(distancias, decreasing = TRUE)

# prueba estadística (coeficiente de variación)
# distribución de probabilidad
repetir <- 10000
cvS <- numeric(repetir)
trozosS <- numeric(n)
for (j in 1:repetir){
  p_azar <- longitud*runif(n + 1)    # elegir los puntos de corte al azar
  p_azar <- sort(p_azar)               # ordenar los cortes
  for (i in 1:n){
    trozosS[i] <- p_azar[i + 1] - p_azar[i]
  }
  trozosS <- sort(trozosS, decreasing = TRUE)
  cvS[j] <- sd(trozosS)/mean(trozosS)
}

# dibujamos el histograma
hist(cvS, freq = FALSE, main = "", border = "lightgrey", xlab = expression("Coeficiente de variación"), ylab = "Densidad de probabilidad")

# vemos dónde quedan los límites de confianza clásicos
limites <- quantile(cvS, probs = c(0.05, 0.95))
abline(v = limites, col = "red", lty = 2)
axis(3, at = limites, labels = round(limites, 2), lwd = 0, cex.axis = 0.8, col.axis = "red")

# calculamos cv en la proteína y lo representamos
cvP <- sd(distancias)/mean(distancias)
abline(v = cvP, col = "darkgreen", lwd = 2)
axis(3, at = cvP, labels = round(cvP, 2), lwd = 0, cex.axis = 0.8, col.axis = "darkgreen")

dev.off()