Hola a todos.

Tengo una pregunta sobre un proyecto que estamos analizando, y que os resumo:


Hemos amplificado el promotor de un gen para ver su metilacion. Para eso hemos
hecho bisulfito y amplificacion con primers adaptados a los cambios de C a T.
Las PCRs han dado regulin, porque daban una banda de tamaño inferior muy
prominente en muchas muestras (pero no hubo manera de mejorarlo, eran primers
que venian de un paper y nos teniamos que ceñir a eso).

Para eliminar los productos de tamaño inferior lo filtre con Prinseq y
quedaron de cientos a miles de lecturas por muestra (hasta aqui todo OK, como
es un solo amplicon nos parece suficiente cobertura).


A continuacion hice un mapeo con Bowtie, usando un amplicon modificado (con
los cambios C a T a excepcion de 5 posibles islas CpG) aplicando las
condiciones Very Sensitive para que aceptara los cambios introducidos por el
bisulfito. Los % de mapeo son como del 50% que para una amplicon tan malo, lo
vemos bien tambien.


Para cuantificar la modificacion C a T de esas 5 islas CpG he usado VarScan2 y
ahi es donde me han surgido dudas.

Hice primero 3 muestras (nº 160 166 y 167) y una vez visto que todo corria
bien, ya hice todas (son casi 300).


Ahora vienen mis dudas /preguntas:


VarScan hace algun tipo de normalizacion de lecturas? (o se puede modificar si
lo hace?) Es que no me ha dado el mismo resultado cuando he analizado esas 3
muestras ellas solas que en el conjunto de las 300. No solo el % de
conversion, sino sobre todo la cobertura total y por alelo en cada posicion.


VarScan (SNPs) es una buena opcion para hacer esto?


En el servidor, dentro de la carpeta GemaDiaz_001-293 estan los datos del
piloto y de los mapeos y llamadas de variantes. Tambien del amplicon
modificado en formato fasta

Os adjunto las tablas de resumen de los resultados y en todo caso, como
entiendo que es un poco complicado, hablamos el lunes y os cuento mejor.


Muchas gracias como siempre por la ayuda


Ricardo




Attachments

Hola Ricardo,

Varscan imagino que la diferencias que reportas es porque estas haciendo es un
pileup donde consideras todas las muestras para hacer una llamada global
considerando todas las muestras como una poblacion.

Esto se hace mucho en variantes somaticas para cazar aquellas variantes de
baja frecuencia. Esa es la razon por la cual no te da el mismo resultado
haciendo 3 que 300, porque haciendo 300 habra algunas

variantes que no pasen el fitro y otras que si. No normaliza sino que aplica
criterios de frecuencias y alineamiento. Yo haria con todas las muestras
porque sera mas fiable que si lo haces con unas pocas.


No obstante, no estoy seguro de que varscan te valga para lo que planteas pero
tampoco veo que no se pueda hacer. De hecho la informacion que quieres la
tienes en el vcf.

Si quieres cuantificar la format variable SPD de format te da las cuentas raw
y si quieres las frecuencias tienes las variables FREQ y si quieres las
cuentas raw que soportan la variante tienes la variable ADF para las cuentas
forward y ADR para las cuentas reversas.

Esto logicamente tienes que interpretarlo por porcentages, por eso es quiza ir
directamente a la variable FREQ.


Ya me dices si eso aclara las dudas o si necesitas que lo hablemos por
telefono

Carlos




Gracias, Carlos.

En realidad no es que vea unas variantes u otras en los dos analisis (entre
hacerlo con unas pocas muestras de prueba o con todas). Es un amplicon pequeño
y en los dos casos estan todas las variantes de las posiciones que queria ver,
pero el problema es que veo cambios muy grandes en la cobertura (tanto ADF
como ADR como el total) entre ambas comparaciones y por consiguiente tambien
en la frecuencia, que es mi dato final.

Se que de partida tengo mucha diferencia de cobertura entre las muestras y sin
embargo eso no se refleja apenas en el vcf (por eso pienso que esta haciendo
algun tipo de normalizacion).

Con SPD, ADF, ADR FREQ etc. te refieres a los encabezados del vcf ¿no? Mi
pregunta es si puedo llegar a hacer el vcf por otra via que sea mas adecuada
en un caso como el mio. He aplicado el mpileup (pero no necesito hacer una
poblacion de mis muestras) asi que eso lo podria cambiar. En principio he
escogido la opcion SNP en lugar de indel aunque quiza sea mejor la combinada?
(no estoy seguro) o quiza sea mejor usar otra herramienta de GATK.

Hoy no tengo mucho acceso a los datos, mañana por la tarde puedo probar alguna
opcion mas y hablamos.

Gracias as always por la ayuda!
ricardo






Hola de nuevo Ricardo,

La cuestion es la misma, cuando haces una llamada de variantes, aplicas una
serie de criterios, filtrado por calidad etc, realineamientos, otras
variantes llamadas te pego aqui el manual de varscan2


|

VarScan - Variant Detection in Massively Parallel Sequencing Data - User’s
Manual

varscan.sourceforge.net




El resultado es general pero si que el vcf reporta las cuentas raws lo que
pasa es que after filters. Por eso, puedes tener tantas diferencias desde una
muestra y otra.


Por eso te decia que no estoy seguro de que una aplicacion para hacer llamada
de variantes sea la mejor opcion aqui, (Varscan o gatk) aunque te podrian
valer. No obstante si lo que quieres

ver es simplemente cuantificar los cambios de C a T, quiza lo ideal seria
hacer un script especifico que considerara el amplicon (no es muy grande la
region, verdad?) dado que siendo una secuencia unica se podria contar las
posiciones que deberian tener C y ver cuales hayan sido sustituidas por T
directamente en los fastq. Guillermo te atreverias tu a hacer un script en
python por ejemplo haciendo algo asi?


Yo creo que eso seria lo mas apropiado, simplemente estarias cuantificando
directamente desde los fastq pero en esencia eso es lo que quieres, ver si el
bisulfito te ha funcionado.


Con ello el script lo que podria hacer SI EL AMPLICON NO ES MUY GRANDE es
considerar es usar una expresion regular basada en la secuencia del amplicon
para definir que posicion en la referencia es la que ocupa la posicion 1 de
cada read. Con ello luego ya seria contar las Ts donde hubiera que haber Cs.


Si el amplicon fuera grande, igual lo suyo seria hacer una tabla hash y
dividir los reads en kmers


Si ves que puedes hacerlo Guillermo dinoslo, si no, si quieres Ricardo te
hacemos nosotros el script y lo ejecutais vosotros por comandos en el
servidor.


Ya me decis


Carlos