Fig_13

pdf("fig_13.pdf", width = 5, height = 5)
par(bty = "n", las = 0, mar = c(4, 4, 1, 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)

# 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)
longitud <- sum(distancias)

# primero calculamos 'd' observada en la proteína
d <- sum(abs(distancias - longitud * sum_arm(n)))

# ahora construimos la distribución de probabilidad
repetir <- 10000
dS <- numeric(repetir)
distanciaa <- numeric(n)
for (j in 1:repetir) {
  p_azar <- nchar(sec) * runif(n + 1)  # elegir los puntos de corte al azar
  p_azar <- sort(p_azar)  # ordenar los cortes
  for (i in 1:n) {
    distanciaa[i] <- p_azar[i + 1] - p_azar[i]
  }
  distanciaa <- sort(distanciaa, decreasing = TRUE)
  longitua <- sum(distanciaa)
  dS[j] <- sum(abs(distanciaa - longitua * sum_arm(n)))
}

# calculamos donde cae el valor observado 'd'
probd <- sum(dS <= d)/repetir

# dibujamos el histograma
hist(dS, freq = FALSE, main = "", border = "lightgrey", xlab = expression("Valor " ~ italic(d)), ylab = "Densidad de probabilidad")
# vemos dónde quedan los límites de confianza clásicos
limites <- quantile(dS, probs = c(0.05, 0.95))
abline(v = limites, col = "red", lty = 2)
# y dibujamos el valor observado 'd' en la proteína
abline(v = d, col = "darkgreen", lwd = 2)

dev.off()