Protocolo reproducible de ddRAD-seq + GWAS

Hola, una tarea frecuente en mejora de plantas es obtener secuencias genómicas de un panel de individuos, por ejemplo variedades de melocotonero o de cebada, para mapear zonas del genoma asociadas a caracteres de interés (GWAS). 

Una manera efectiva de hacerlo es secuenciar solamente una fracción del genoma, una que maximice la información obtenida, para reducir costes. Najla Ksouri, una joven investigadora de nuestro campus, acaba de publicar un protocolo para esta tarea. En este caso se basa en la metodología ddRAD-seq, que secuencia fragmentos de ADN obtenidos tras digerir ADN nuclear con una pareja de enzimas de restricción elegidas para representar todo el genoma y minimizar las regiones repetidas.  

El protocolo está disponible en https://github.com/najlaksouri/GWAS-Workflow y comprende las siguientes etapas, que a su vez requieren de scripts en el repositorio: 

El protocolo usa como ejemplo datos de melocotonero, publicados en https://doi.org/10.1186/s13007-025-01415-3, pero en principio se puede aplicar a otras especies diploides sin modificaciones. 

La siguiente figura muestra un ejemplo de picos significativos obtenidos por GWAS en varios cromosomas, en este caso para fecha de cosecha:

 

Ejemplo de picos de GWAS (A), bloque ligado (B), carácter asociado (C) y genes candidato (D), tomado de https://doi.org/10.1186/s13007-025-01415-3

Un saludo

líneas distintas entre ficheros con grep

Si usas el terminal de linux sabrás de la utilidad de grep para encontrar de manera eficiente cadenas de caracteres o expresiones regulares en ficheros. Si no, tienes ejemplos por ejemplo en nuestro material. En mi caso, aunque usuario habitual, me he tropezado en diferentes ocasiones con el siguiente problema: dame todas las líneas del ficheroA no encontradas en ficheroB. Encontré la mejor solución aquí.

Opciones frecuentes de grep, https://wizardzines.com/comics/grep

Explico el problema con un ejemplo común en bioinformática. Imagina un fichero con secuencias (ficheroA.fasta) y otro con resultados de análisis de esas secuencias (ficheroB.tsv). Ahora quieres averiguar qué secuencias de A no están presentes en B, por ejemplo para repetir el análisis o arreglar errores en el código.

El ficheroA contiene las siguientes líneas:

>100007_TR35452-c0_g1_i1
CAATTTACGCCTATCGTTATCCATTTCTA...

>10000_TR33868-c0_g1_i1
GGGGGACCTACTCAAATCCCCATCTCCC...

>10001_TR436-c0_g1_i1
GTTTCCAACCGGATGTTGAAACAGACAA...

El ficheroB en cambio puede ser un formato  TSV, por ejemplo:

#metadatos, nombres de columnas, etc
100007_TR35452-c0_g1_i1 chr5H   540009332       540009636
1000_TR868-c0_g2  chr4H   340992292       340995709
...

Resuelvo el problema en dos comandos en el terminal:

1) extraigo de A únicamente los nombres de las secuencias:

$ perl -lne 'if(/>(\S+)/){print $1}' ficheroA.fasta > ficheroA.nombres

2.i) busco las secuencias reportadas en B para luego 2.ii) buscar las secuencias de A que no están en B:

$ grep -Fo -f ficheroA.nombres ficheroB.tsv | grep -vFf - ficheroA.nombres

 

Un saludo

Descarga datasets del NCBI desde el terminal

Hola, hace unos días necesitaba obtener todas las secuencias de genes de pimiento (Capsicum annuum) de la colección 'gene' del NCBI (que tan mal lo está pasando en 2025). Para ello abrí el navegador y obtuve una lista de 48K identificadores (gene-id) en https://www.ncbi.nlm.nih.gov/gene?term=capsicum%20annuum%5BOrganism%5D

Descargué el fichero y lo hice no redundante de esta manera:

$ head -3 gene_result.txt | cut -f 1-6
tax_id	Org_name	GeneID	CurrentID	Status	Symbol
4072	Capsicum annuum	107859632	0	live	LOC107859632
4072	Capsicum annuum	107868427	0	live	LOC107868427 

$ cut -f 3 gene_result.txt | sort -u | grep -v Gene > pimiento.geneids.txt

Descubrí que obtener la lista de genes es fácil, pero no tanto sus secuencias. Por ejemplo no lo logré con https://www.ncbi.nlm.nih.gov/sites/batchentrez , me daba secuencias que no eran de pimiento, supongo que por esperar otro tipo de identificadores. Entonces pedí ayuda en https://support.nlm.nih.gov/support/create-case y tras unos días de espera me dirigieron amablemente a la documentación de la herramienta datasets del NCBI, y me compartieron la siguiente figura:

Me descargué el binario datasets para linux de https://www.ncbi.nlm.nih.gov/datasets/docs/v2/command-line-tools/download-and-install y lo utilicé de esta manera:

$ chmod +x datasets

# probamos primero con un gen de ejemplo 

$ ./datasets download gene gene-id 20217883 --include gene,protein
Collecting 1 gene record [================================] 100% 1/1
Downloading: ncbi_dataset.zip    3.24kB valid data package
Validating package files [================================] 100% 5/5
$ unzip ncbi_dataset.zip 
Archive:  ncbi_dataset.zip
  inflating: README.md               
  inflating: ncbi_dataset/data/gene.fna  
  inflating: ncbi_dataset/data/data_report.jsonl  
  inflating: ncbi_dataset/data/dataset_catalog.json  
  inflating: md5sum.txt              
$ head ncbi_dataset/data/gene.fna
>NC_024624.1:447314-449261 rrn18 [organism=Capsicum annuum] [GeneID=20217883] [chromosome=MT]
ATCATAGTCAAAAGAAGAGTTTGATCCTGGCTCAGAAGGAACGCTAGCTATATGCTTAACACATGCAAGT
CGAACGTTGTTTTCGGGGAGCTGGGCAGAAGGAAAAGAGGCTCCTAGCTAAAGGTAGCTTGTCTCGCCCA
GGAGGTGAGAAGAGTTGAGAACAAAGTGGCGAACGGGTGCGTAACGCGTGGGAATCTGCCGAACAGTTCG
GGCCAAATCCTGAAGAAAGCTAAAAAGCGCTGTTTGATGAGCCTGCGTAGTATTAGGTAGTTGGTCAGGT
AAAGGCTGACCAAGCCAATGATGCTTAGCTGGTCTTTTCGGATGATCAGCCACACTGGGACTGAGACACG
GCCCGGACTCCCACGGGGGGCAGCAGTGGGGAATCTTGGACAATGGGCGAAAGCCCGATCCAGCAATATC
GCGTGAGTGAAGAAGGGCAATGCCGCTTGTAAAGCTCTTTCGTCGAGTGCGCGATCATGACAGGACTCGA
GGAAGAAGCCCCGGCTAACTCCGTGCCAGCAGCCGCGGTAAGACGGGGGGGGCAAGTGTTCTTCGGAATG
ACTGGGCGTAAAGGGCACGTAGGCGGTGAATCGGGTTGAAAGTGAAAGCCGCCAAAAACTGGCGGAATGC
 

# ahora con gene-ids del fichero que preparamos antes, tarda, ncbi_dataset.zip > 100MB

$ ./datasets download gene gene-id --inputfile pimiento.geneids.txt --include gene,protein

$ unzip ncbi_dataset.zip 
Archive:  ncbi_dataset.zip
  inflating: README.md               
  inflating: ncbi_dataset/data/gene.fna  
  inflating: ncbi_dataset/data/protein.faa  
  inflating: ncbi_dataset/data/data_report.jsonl  
  inflating: ncbi_dataset/data/dataset_catalog.json  
  inflating: md5sum.txt              
 

$ head ncbi_dataset/data/gene.fna
>NW_025826840.1:c5280-2277 LOC124890618 [organism=Capsicum annuum] [GeneID=124890618] [chromosome=Un]
GGAGGTACTTAAAGCATGACTTTTAAAAGTTTGCATAGGGCAAGAAGCAGGAGTGTGACTAAACTGATTT
TTCTTTCTGGTTTTAAGCATGATGCTATTCCTCAGCGTCCTCAAAATGAGCAAATTGAAAAGCTCAAGAA
GTTCAAGGCTGTGTTGGAACGCATTCTGATTTTCTTGCAGCTCAATAAGCATGACATTCAGCTTACTCAC
AAGGAGAAGTTGTGTTCGGTTGAGAGGCACATAGGTTTCTTTCTTAGCAAGCCTACTTCTCCTCCTCTGC
AGGGGCAACTTCCTCAGTCTTCCATGCAGCTTCAGCAACCACAATCACTTGATGTTCAAACTAATCCACC
GATGCAACCTCAACTTCATCAGGCACTATCTTCGCAGGTACGTCATCAACATTTTAATCCACTATTATCA
TTTCTGGAGGCAATTCTACCAATTGTATGGTGCATCATGCTGGATTTACTAAATTTTGATACTATAAAAG
GTCTCTTGACAAACAGCTAGCCAAATGTGTCAGTCGTCATTGAAAGTTCTGCTACCGTTTAGTTTCTTTT
TCTCCAATGTCTTTTGTTACACTTGTTTTGTATATTACTATATTGTTGGCTCTTTTTCTTTCTTTTGATC

$ head ncbi_dataset/data/protein.faa 
>XP_047258369.1 LOC124890618 [organism=Capsicum annuum] [GeneID=124890618]
MTFKSLHRARSRSVTKLIFLSGFKHDAIPQRPQNEQIEKLKKFKAVLERILIFLQLNKHDIQLTHKEKLC
SVERHIGFFLSKPTSPPLQGQLPQSSMQLQQPQSLDVQTNPPMQPQLHQALSSQAQSTGALQTATLDSDS
TSQTGNADGADWQEELYQEIKTMREKNLPELNALYQKIASKVQQHDAIPQRPQNEQIEKLKMFKAVLERI
LIFLQVNKHDIQLTHKEKLCSVERHIGFFLSKPTSPPLQGQLPQSSMQLQQPQSLDVQTNPPMQPQLHQA
LSSQAQSTGALQTATLDSDSTSQTGNADGADWQEELYQEIKTMRDKNLPELNA
>XP_047259376.1 LOC124891834 [organism=Capsicum annuum] [GeneID=124891834]
MTSNITESLNSILRDEREYPVASIFNSIAPRFGEIFRKRYAEVDNSKTTFIPVAETILRENMTKGDKLYV
NNINESTNEFTVLGYGRSAKVNLSRQPCSCRKYDLVKLPCAYTMAALHLKHGDEYGTSIYKNPFQIYSKE
SYLLAYLEPICAAPLESEWSVAREYLEIQVLPPDVDPKHGRRKVKHVKGVLEPSRYKKRNKCSKCKRLGH

Hay muchas otras maneras de utilizarlo y ejemplos en https://www.ncbi.nlm.nih.gov/datasets/docs/v2/how-tos

Hasta pronto, Bruno