Pregunta:
dividir el genoma en ventanas que no se superponen usando R
Anna1364
2017-11-24 11:25:08 UTC
view on stackexchange narkive permalink

Tengo casi 10 millones de SNP ubicados en 10 cromosomas. Quiero dividir el genoma en ventanas no superpuestas de 15, 20 y 30 kb. Aquí está parte de mi tabla SNP: 

  head (sap_ids) snp_id chr pos Chr01__15043 1 15043 Chr01__15079 1 15079 Chr01__15139 1 15139 Chr01__15165 1 15165 ... ... ... Chr17__214708424 17 214708424 Chr17__214708451 17 214708451 Chr17__214708484 17 214708484 Chr17__214708508 17 214708508 Chr17__214708574 17 214708574

He estado usando este código, pero me da una salida incorrecta porque cada cromosoma comienza en la posición 1 

  win_size<-c (15000,30000,50000,100000) res<- cbind (snp_ids, data.frame (lapply (setNames (win_size, paste ("ventana", win_size, sep = "_")), función (x) as.numeric (techo (snp_ids $ pos / x)))))

Por ejemplo, como puede ver en el siguiente ejemplo en la ventana 4, obtengo SNP de los cromosomas 1, 3 y 5: 

  snp_id chrome poistion windowChr01__58332 1 58332 4Chr01__58335 1 58335 4Chr01__58341 1 58341 4Chr01__58450 1 58450 4Chr01__58471 1 58471 4Chr01__58530 1 58530 4Chr01__58542 1 58542 4Chr01__58641 1 58641 4Chr03__45457 3 45457 4Chr03__45604 3 45604 4Chr04__56873 4 56873 4Chr04__57387 4 57387 4Chr04__57399 4 57399 4Chr04__57528 4 57528 4Chr04__58419 4 58419 4Chr04__59670 4 59670 4Chr04__59704 4 59 704 4

Supongo que la solución sería ejecutar el código anterior recorriendo cada cromosoma. Probé algo como esto: 

  for (i in 1:10) {res <- cbind (gatos, data.frame (lapply (setNames (win_size, paste ("ventana", win_size, sep = "_")), función (w)
paste (gatos $ chr [gatos $ chr == i], techo (gatos $ pos / w), sep = "_"))))}

Pero aún no funciona correctamente ( la columna 4 debe ser 1_1, para el chr 1)! 

  Chr01__912 1 912 2_1Chr01__944 1 944 2_1Chr01__1107 1 1107 2_1Chr01__1118 1 1118 2_1

¿Alguien puede ayudarme para averiguar cómo puedo modificar el código para crear estas ventanas que no se superponen correctamente.

¿Por qué tiene 4 campos en `sap_ids` para chr17 y 3 para chr1? ¿Por qué cambia eso? ¿Sus posiciones aumentan con el número de cromosomas (son todas las posiciones en cr2 mayores que aquellas en cr1) o cuenta cada cromosoma por separado (por lo que cada cromosoma comienza desde la posición 1)? Si es lo último, tiene sentido que su enfoque falle.
No obtengo la pieza con ventanas de 15, 20 y 30 kb. ¿Quiere que las ventanas sean de varios tamaños dependiendo de la cantidad de SNP en el interior o desea tener un tamaño de ventana que corresponda en promedio a ~ 200 SNP?
@terdon, Acabo de editar un poco mi publicación. ver mis ediciones. sí, en mi caso cuento cada cromosoma por separado, lo que significa que cada cromosoma comienza en la posición 1. ¿Alguna idea de cómo hacerlo en R? @ Kamil S Jaron, ¡solo quiero hacer mi análisis con diferentes tamaños de ventana! Eso es.
@Anna Para evitar eliminar las respuestas, volví a la versión anterior. Si su pregunta no tiene una respuesta completa, [edítela] o comente explicando por qué no. Siempre puede anular la selección de la respuesta aceptada.
¡Deje de cambiar radicalmente su pregunta! Si necesita más ayuda, haga una nueva pregunta que explique qué más necesita. Pero dado que esto ha sido respondido, cambiarlo hace que las respuestas sean obsoletas.
One responder:
Sebastian Müller
2017-11-25 20:47:56 UTC
view on stackexchange narkive permalink

La pregunta es un poco confusa para mí, en esencia, entiendo que desea ventanas que no se superpongan entre los cromosomas.

Una forma de lograrlo es usar GenomicRanges :: tileGenome función, que necesita las longitudes de los cromosomas como entrada, por ejemplo, 

  biblioteca (Biostrings) mygenome <- readDNAStringSet (list.files (mypath, "mygenome.fa $", full = TRUE)) chrSizes <- width (mygenome) names (chrSizes) <- names (mygenome) print (chrSizes) # Chr1 Chr2 Chr3 Chr4 Chr5 # 18585056 19698289 23459830 26975502 30427671

Esto puede ser introducido en tileGenomes : 

  biblioteca (GenomicRanges) bins <- tileGenome (chrSizes, tilewidth = 5e5, cut.last.tile.in.chrom = T)

Este crea un objeto GRange que le brinda un repertorio completo de opciones de posprocesamiento, es decir, superpuesto con lecturas de SNP o NGS (utilizo este último para gráficos de cobertura amplia del genoma): 

  print (bins) # Objeto GRanges con 240 rangos y 0 columnas de metadatos: # rangos de seqnames strand # <Rle> <IRanges> <Rle> # [1] Chr1 [1, 500000] * # [2] Chr1] [500001 * 3 1000000] Chr1 [1500001, 2000000] * # [5] Chr1 [2000001, 2500000] * # ... ... ... ... # [236] Chr5 [28000001, 28500000] * # [237] Chr5 [28500001, 29000000] * # [238] Chr5 [29000001, 29500000] * # [239] Chr5 [29500001, 30000000] * # [240] Chr5 [30000001, 30427671] * # ------- # seqinfo: 5 secuencias de un genoma no especificado
En relación con esto, existe el comando `tile ()` de GenomicRanges que producirá un objeto GRangesList de la entrada. No producirá exactamente lo que OP deseaba, pero puede ser una alternativa más conveniente.
Muchas gracias @Sebastian Müller. si. En realidad, este es un enfoque mejor que escribir código para eso.
¡Me alegro de que sea útil! ¿Quizás podría hacer su pregunta de SNP en una pregunta diferente y sacarla aquí para que la pregunta real esté más en línea con el título (que es casi autoexplicativo por sí mismo)?
También para completar, también existe la función `GenomicRanges :: slideWindows` que toma un objeto GRanges diseñado para dar como resultado ventanas superpuestas, pero también se puede hacer que produzca ventanas no superpuestas estableciendo` width` = `step`


ⓘ
Esta pregunta y respuesta fue traducida automáticamente del idioma inglés.El contenido original está disponible en stackexchange, a quien agradecemos la licencia cc by-sa 3.0 bajo la que se distribuye.
Loading...