1 Independencia y Homogeneidad

En este ejemplo generamos muestras sintéticas de cadenas de longitud 1000 de las bases “a”,“c”,“g” y “t” y calculamos la base de máxima frecuencia.

Primero fijamos la semilla de generación aleatoria por [reproducibilidad}(https://es.wikipedia.org/wiki/Reproducibilidad_y_repetibilidad)

set.seed(31415)

La generación se produce en tres escenarios/tipos distintos “A”, “B” y “C”. Cada escenario aparece con una probabilidad 0.4, 0.3 y 0.3 respectivamente. En cada escenario las probabilidades de aparación de las cuatro bases es distinta:

Generaremos 10000 muestras de estos tres escenarios en función de su probabilidad.

2 Generación de las muestras.

Generamos los tipos en función de sus proporciones.

tipo=sample(x=c("A","B","C"),size=10000,replace=TRUE,prob=c(0.4,0.3,0.3))
head(tipo)
## [1] "B" "A" "A" "C" "A" "A"
str(tipo)
##  chr [1:10000] "B" "A" "A" "C" "A" "A" "A" "C" "B" "B" "A" "B" "B" "A" ...
table(tipo)
## tipo
##    A    B    C 
## 4026 2997 2977

Ahora calculamos para cada uno la base que predomina (si hay empate escogemos la primera que aparece):

bases=c("a","c","g","t")
tipo.probs=data.frame(
  p.a=c(0.25,0.26,0.24),
  p.c=c(0.25,0.26,0.25),
  p.g=c(0.25,0.24,0.26),
  p.t=c(0.25,0.24,0.25)
)
tipo.probs
##    p.a  p.c  p.g  p.t
## 1 0.25 0.25 0.25 0.25
## 2 0.26 0.26 0.24 0.24
## 3 0.24 0.25 0.26 0.25
row.names(tipo.probs)=c("A","B","C")

tipo.probs["A",]
##    p.a  p.c  p.g  p.t
## A 0.25 0.25 0.25 0.25
tipo.probs["B",]
##    p.a  p.c  p.g  p.t
## B 0.26 0.26 0.24 0.24
tipo.probs["C",]
##    p.a  p.c  p.g  p.t
## C 0.24 0.25 0.26 0.25
base.predominante= function(tipo,tipo.probs)
  {aux=which.max(table(sample(bases,size=1000,replace=TRUE,prob=tipo.probs[tipo,])))
   names(aux)[1]
}

# Ejemplos generación muestra
base.predominante("A",tipo.probs)
## [1] "t"
base.predominante("A",tipo.probs)
## [1] "t"
base.predominante("B",tipo.probs)
## [1] "g"
base.predominante("B",tipo.probs)
## [1] "c"
base.predominante("C",tipo.probs)
## [1] "g"
base.predominante("C",tipo.probs)
## [1] "t"

Ahora generamos la muestra para cada una de las observaciones de la muestra tipos que ya hemos generado.

muestra=data.frame(tipo)
muestra$max.frec=unlist(lapply(muestra$tipo,FUN=base.predominante,tipo.probs))

head(muestra)
##   tipo max.frec
## 1    B        g
## 2    A        a
## 3    A        a
## 4    C        c
## 5    A        t
## 6    A        a
str(muestra)
## 'data.frame':    10000 obs. of  2 variables:
##  $ tipo    : Factor w/ 3 levels "A","B","C": 2 1 1 3 1 1 1 3 2 2 ...
##  $ max.frec: chr  "g" "a" "a" "c" ...
##Guardar datos
#Pon tu path
setwd("./")

write.table(muestra,file="MuestraTotalBases.txt",sep="\t",row.names=FALSE)

write.table(muestra[muestra$tipo=="A",],file="muestraTipoAbases.txt",sep="\t",row.names=FALSE)
write.table(muestra[muestra$tipo=="B",],file="muestraTipoBbases.txt",sep="\t",row.names=FALSE)
write.table(muestra[muestra$tipo=="C",],file="muestraTipoCbases.txt",sep="\t",row.names=FALSE)

Podeis bajar los y el código datos aquí..