Generación de muestras sintéticas para pruebas de independencia y homogeneidad

En este ejemplo generamos una muestra sintética formada por 10000 cadenas de longitud 1000 de las bases “a”, “c”, “g” y “t” y calculamos la base de máxima frecuencia de cada cadena.

Vamos a producir cadenas de tres tipos distintos: “A”, “B” y “C”. Estos tipos van a caracterizarse por sus probabilidades de aparición de las cuatro bases “a”, “c”, “g” y “t”:

Vamos a imponer que los tipos de cadena “A”, “B” y “C” se den con probabilidades 0.4, 0.3 y 0.3 respectivamente. Así pues, generaremos nuestra muestra de cadenas al azar teniendo en cuenta estas probabilidades.

Generación de las muestras

Generamos al azar los tipos de cadenas 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))

Demos un vistazo al vector de tipos de cadenas que hemos generado:

str(tipo)

##  chr [1:10000] "B" "A" "A" "C" "A" "A" "A" "C" "B" "B" "A" "B" "B" "A" ...

head(tipo)

## [1] "B" "A" "A" "C" "A" "A"

table(tipo)

## tipo
##    A    B    C 
## 4026 2997 2977

Las frecuencias de los tres tipos han quedado muy cerca de las esperadas.

Ahora especificamos las probabilidades de las bases que definirán las cadenas de cata tipo:

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)
)
row.names(tipo.probs)=c("A","B","C")
tipo.probs

##    p.a  p.c  p.g  p.t
## A 0.25 0.25 0.25 0.25
## B 0.26 0.26 0.24 0.24
## C 0.24 0.25 0.26 0.25

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

La función base.predominante siguiente genera una cadena de longitud 1000 del tipo entrado como argumento y calcula su base de frecuencia máxima (si hay empate en el máximo, escoge una al azar):

base.predominante= function(tipo,tipo.probs)
  {TT=table(sample(bases,size=1000,replace=TRUE,prob=tipo.probs[tipo,])) 
  aux=names(which(TT==max(TT)))
  sample(aux,1)
}

Veamos algunos ejemplos de aplicación de esta función:

# 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 de pares (tipo de cadena, base predominante), aplicando la función base.predominante a cada uno de los tipos de cadenas que forman la muestra tipo que hemos generado al principo. Organizamos los resultados en un data frame llamado muestra.

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

Demos un vistazo a la tabla de datos que hemos generado:

head(muestra)

##   tipo max.frec
## 1    B        g
## 2    A        a
## 3    A        c
## 4    C        t
## 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" "c" "t" ...

Finalmente, creamos los ficheros descargables en el directorio actual:

El fichero “MuestraTotalBases.txt” contiene la tabla de datos completa.
Los ficheros “MuestraTipoAbases.txt”, “MuestraTipoBbases.txt” y “MuestraTipoCbases.txt” contienen, respectivamente, las filas de la muestra global que corresponden a cadenas de tipo “A”, de tipo “B” y de tipo “C”.

En las filas de estos cuatro ficheros, las dos variables están separadas por una coma.

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

write.table(muestra[muestra$tipo=="A",],file="MuestraTipoAbases.txt",sep=",",row.names=FALSE)
write.table(muestra[muestra$tipo=="B",],file="MuestraTipoBbases.txt",sep=",",row.names=FALSE)
write.table(muestra[muestra$tipo=="C",],file="MuestraTipoCbases.txt",sep=",",row.names=FALSE)