Evaluación estadística de las bases de datos de frecuencias

Este Notebook está diseñado para ejecutar mispitools y evaluar el poder estadístico de una base de datos de frecuencias alélicas. Primero, instalamos las herramientas necesarias:

#install.packages("mispitools")
#install.packages("pedtools")
rmarkdown::render("~/Escritorio/mispitools4refs.Rmd")

Ahora cargamos las bibliotecas.

library(mispitools)
library(pedtools)
library(forrel)

Ahora configuramos nuestros sistemas de prueba. En este caso, vamos a analizar los siguientes parentescos:

• Paternidad.
• Tío.
• Hermanos.
• Medio hermanos.
• Abuelos.

paternity <-linearPed(1)
uncle <- avuncularPed()
sibling <- nuclearPed(2)
halfSib <- halfSibPed()
grandparent <-linearPed(2)

par(mfrow = c(2, 3))
plot(paternity, main = "Paternity", hatched = 1,  fill = c(`3` = "2"))
plot(uncle, main = "Avuncular", hatched = 3,  fill = c(`6` = "2"))

plot(sibling, main = "Sibling", hatched = 3,  fill = c(`4` = "2"))
plot(halfSib, main = "Half Sibling", hatched = 4,  fill = c(`5` = "2"))

plot(grandparent, main = "Grandparent", hatched = 1,  fill = c(`5` = "2"))

Tenga en cuenta que los individuos representados con líneas discontinuas serán los genotipados, y los individuos en rojo serán las Personas de Interés.

Trabajando con la base de datos

Primero, debemos proporcionar una base de datos de frecuencias alélicas. Por defecto, mispitools cuenta con bases de datos de todo el mundo. Además, si es necesario agregar una nueva base de datos, debemos tener en cuenta el formato. Veamos la base de datos de frecuencias alélicas de Argentina:

Argentina

La primera columna (“Allele”) es obligatoria y debe tener ese nombre. Contiene todos los alelos observados (en todos los marcadores). Luego, cada columna representa un marcador, y los valores de frecuencia deben sumar 1 en todos ellos. La siguiente función reorganiza la base de datos a un formato más amigable (desde el punto de vista computacional). Dicho esto, se puede agregar una nueva base de datos usando, por ejemplo, read_csv, teniendo en cuenta estos detalles de formato. Después de leerla, deberíamos poder visualizarla de la misma manera que vemos la base de datos argentina.

Arg <- getfreqs(Argentina)
head(Arg, 1) # showing one marker

$D8S1179
    2.2     3.2       4     4.2       5       6       7       8     8.3       9     9.1     9.2 
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00898 0.00000 0.00762 0.00000 0.00000 
    9.3      10    10.2    10.3      11    11.2    11.3      12    12.1    12.2    12.3      13 
0.00000 0.06872 0.00000 0.00000 0.07247 0.00000 0.00000 0.14353 0.00000 0.00000 0.00000 0.30446 
   13.2    13.3      14    14.2    14.3      15    15.2    15.3      16    16.2    16.3    16.4 
0.00000 0.00000 0.22236 0.00000 0.00000 0.13732 0.00000 0.00000 0.03047 0.00000 0.00000 0.00000 
     17    17.1    17.3      18    18.1    18.2    18.3      19    19.1    19.2    19.3      20 
0.00337 0.00000 0.00000 0.00065 0.00000 0.00000 0.00000 0.00005 0.00000 0.00000 0.00000 0.00000 
   20.1    20.2    20.3      21    21.2    21.3      22    22.2    22.3      23    23.2    23.3 
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
     24    24.2    24.3      25    25.2    25.3      26    26.2      27    27.2      28    28.2 
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
     29    29.2      30    30.2      31    31.2      32    32.1    32.2      33    33.1    33.2 
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
   33.3      34    34.2      35    35.1    35.2      36    36.2      37 
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 

paternity <-setMarkers(paternity, locusAttributes = Arg[1:22]) # Database Arg, after processing, with 22 markers
uncle <- setMarkers(uncle, locusAttributes = Arg[1:22]) 
sibling <- setMarkers(sibling, locusAttributes = Arg[1:22]) 
halfSib <- setMarkers(halfSib, locusAttributes = Arg[1:22]) 
grandparent <-setMarkers(grandparent, locusAttributes = Arg[1:22]) 

Ahora podemos simular nuestros genotipos. Comenzamos simulando los individuos de referencia (representados con líneas discontinuas en el árbol genealógico).

numsims = 50 # Here put the number of pairs to be simulated
paternity <-profileSim(paternity, N = numsims, ids = 1)
uncle <- profileSim(uncle, N = numsims, ids = 3)
sibling <- profileSim(sibling, N = numsims, ids = 3)
halfSib <- profileSim(halfSib, N = numsims, ids = 4)
grandparent <-profileSim(grandparent, N = numsims, ids = 1)

Bucle Importante - Permite simular LRs a partir de pares de individuos.

Para cada individuo de referencia, simulamos una Persona de Interés (ten en cuenta que podemos modificar el código fácilmente para simular varias Personas de Interés por cada individuo de referencia; en este caso, estamos simulando 2 por referencia, por lo que con 50 referencias obtenemos 100 LRs).

NOTA IMPORTANTE: Después de ejecutar, si deseas volver a realizar la simulación, necesitarás re-ejecutar desde el bloque de gráfico del árbol genealógico (esto se debe a que es necesario limpiar el objeto ped).

Y ahora podemos extraer los resultados utilizando algunas funciones personalizadas:

# Función para aplicar simLR2dataframe a cada simulación dentro de la lista y combinar los resultados
combine_sim_results <- function(sim_list) {
  # Crear una lista vacía para almacenar los dataframes de cada simulación
  df_list <- list()
  
  # Iterar sobre cada simulación dentro de la lista sim_list
  for (i in 1:length(sim_list)) {
    # Aplicar la función simLR2dataframe a cada simulación
    df <- simLR2dataframe(sim_list[[i]])
    df_list[[i]] <- df  # Guardar el dataframe en la lista
  }
  
  # Combinar todos los dataframes en uno solo
  combined_df <- do.call(rbind, df_list)
  
  return(combined_df)
}

# Crear dataframes combinados para cada relación de parentesco
df_patsim_combined <- combine_sim_results(patsim)
df_unclesim_combined <- combine_sim_results(unclesim)
df_sibsim_combined <- combine_sim_results(sibsim)
df_halfsim_combined <- combine_sim_results(halfsim)
df_grandsim_combined <- combine_sim_results(grandsim)

Ahora puedes visualizar y analizar fácilmente los datos obtenidos :)

Los Related son los LRs obtenidos bajo H1, y los Unrelated son los LRs obtenidos bajo H2.

head(df_patsim_combined)

Contamos con varias métricas para evaluar el rendimiento. Por favor, consulta ITpaper y Decisions para obtener más información.

Conclusiones

En este análisis, hemos evaluado el poder estadístico de una base de datos de frecuencias alélicas utilizando diferentes relaciones de parentesco. Los resultados obtenidos nos permiten comprender mejor cómo los marcadores genéticos seleccionados contribuyen a la discriminación entre hipótesis de parentesco y no parentesco.

Este enfoque puede ser extendido y adaptado para incorporar nuevas bases de datos, diferentes marcadores genéticos y otras relaciones de parentesco, proporcionando una herramienta flexible para la genética forense y estudios de parentesco.