Recerca de selenoproteïnes
en Equus Przewalskii

Finalment, aquestes permeten executar els programes T-Coffee i Genewise respectivament.

Per agilitzar el procés es pot utilitzar el script exports.sh, executant-lo abans d'iniciar el procés amb la següent ordre:

$ source exports.sh

El Basic Local Alignment Search Tool (BLAST) és un programa que ens permet comparar seqüències biològiques, ja siguin proteiques o de DNA. Es tracta d'un algoritme que troba similituds entre dues seqüències, realitzant alineaments locals de forma heurística entre aquestes. S'ha de tenir en compte que el fet d'utilitzar un algoritme heurístic pot comportar la pèrdua de hits reals que no presenten una similitud molt elevada.

Dels diferents tipus de BLAST, en aquest cas particular ens interessa el tblastn, que permet comparar una proteïna query amb una seqüència de nucleòtids o una base de dades de nucleòtids (genome.index). El nostre objectiu és, mitjançant aquesta eina, localitzar les regions on es troben els gens d'interès que probablement codificarant per les nostres querys.

Com a resultat obtindrem els alineaments possibles de la query i el genoma. Per a cada alineament se'ns presentarà el hit, és a dir, la regió que s'alinea amb la notsra seqüència problema, amb un determinat E-value i un score. L'E-value ens indica la probabilitat d'obtenir aquell alineament únicament per atzar, de manera que, com menor sigui el valor de l'E-value, menys probabilitat hi ha que aquell alineament s'hagi produit per atzar. S'acostumen a considerar significatius valors d'E-value a partir de 10^-4. Pel que fa al score en canvi, valors més grans indicaran un millor alineament. Tot i això, a l'hora d'avaluar un hit, ens basarem principalment en paràmetres estadístics com l'E-value que ens informin de la seva significança.

Per a executar el BLAST utilitzem la següent ordre des del terminal:

$ blastall -p tblastn -e 1e-3 -m 9 -i query.fa -d /cursos/BI/genomes/vertebrates/2014/Equusprzewalskii/genome.fa -o query.blast

On:

• -p fa referència al tipus de blast.
• -e indica el valor de E-value que posem com a llindar. El fet de treballar amb espècies tan properes permet utilitzar llindars més restrictius, per exemple, es podrien considerar només els hits amb un E-value inferior a 1e-10.
• -m indica com s'organitzarà l'output del programa. En el nostre cas i per facilitar la seva interpretació i posterior utilització, la informació obtinguda del BLAST estarà organitzada en columnes, amb una primera línia indicant el contingut de cadascuna de les columnes.

El següent pas un cop realitzat el BLAST és obtenir els hits, és a dir, la regió genòmica amb la qual ha estat alineada la nostra query. Per fer-ho, però, primer és necessari tenir un arxiu amb la seqüència genòmica per cadascun dels scaffolds d'interès d'Equus przewalskii. Per obtenir aquests arxius farem servir el programa Fastafetch mitjançant la següent comanda:

$ fastafetch /cursos/BI/genomes/vertebrates/2014/Equusprzewalskii/genome.fa genome.index "identificador" > scaffold.fa

En aquest cas utilitzarem el fitxer "genome.index", on el genoma d'Equus przewalskiies troba organitzat en scaffolds, facilitant així el procés.

A continuació, l'objectiu és, com ja hem dit, extreure les regions genòmiques que probablement codifiquen per les nostres querys, i que per tant voldrem avaluar. El primer pas és doncs delimitar el hit per tal d'obtenir la regió el més acotada possible que contingui el gen d'interès. Aquest punt serveix per facilitar el procés, ja que ens evita treballar amb seqüències innecessàriament llargues.

L'output del BLAST ens indica, per cada alineament realitzat, la posició d'inici i final d'aquest, tant a nivell d'aminoàcids de la query com a nivell de nucleòtids del genoma d'Equus przewalski. Cal tenir en compte que l'alineament no és perfecte, i que és possible que part del gen no es trobi en la regió alineada amb la proteïna. Per aquest motiu, és molt important agafar posicions del hit upstream i downstream. D'aquesta manera, ens assegurem que el gen d'interés es trobi a l'interior de la regió delimitada.

La grandària d'aquesta expansió dependrà de la llargada del scaffold, i de les posicions de l'alineament. Agafarem un màxim de 1000 nucleòtids upstream, i 1000 nucleòtids downstream. En cas que, en restar 1000 a la posició d'inici de l'alineament (start position), el seu valor esdevingui negatiu, agafarem la primera posició del scaffold com a inici del hit (start). D'altra banda, en cas que, en sumar 1000 a la posició final de l'alineament, aquesta sigui superior a la última posició del scaffold, s'agafarà aquesta com a última posició del hit.

Podeu trobar el codi en el següent enllaç: codi.pl (previsualitza). Cal executar-lo des del terminal dins la carpeta "Inputs" utilitzant la comanda en bash: codi.sh.

Aquest script presenta certes limitacions. D'una banda, el plantejament inicial, en què es va considerar que tenir els inputs organitzats en famílies facilitaria la interpretació, va suposar un problema a l'hora d'executar el programa Exonerate. És per això que es va corregir de manera que, a partir dels inputs per famílies, es realitza el blast i que, en funció de la query que s'analitza, s'agafi el corresponent fitxer (aquell que tingui el mateix ID) de la carpeta "Inputs_separats".

D'altra banda, degut al fet d'assignar per a cada hit la seva millor query perdem resultats en casos que dues proteïnes de la mateixa família es trobin en un mateix scaffold. Degut a que aquestes acostumen a presentar una gran homologia entre elles, en l'output del blast obtindrem hits tant en la regió del scaffold on realment es localitza aquella proteïna com per les regions on es trobi l'altra proteïna de la famïlia. Aixó fa que el subseq es generi des d'na start position corresponent a la regió més pròxima a l'extrem 5' del scaffold i una stop position corresponent a la regió situada més pròxima al 3'. El fet d'obtenir de manera incorrecta el subseq impedeix una correcta predicció d'ambdues proteïnes.

Degut a restriccions de temps, no hem pogut corregir aquesta limitació i el que hem fet, en canvi, és que en cada cas particular en que els resultats obtinguts eren erronis o incomplets hem dut a terme el procés manualment. Tot i això, creiem aquest podría ser un bon punt a tenir en compte per a futures millores en el script.

Una opció en cas de trobar-se fora de la FCSV és utilitzar el cluster de docència proporcionat pels professors de l'assignatura, i format per dos ordinadors amb un total de 16 processadors encesos en tot moment. En aquest equipament informàtic anomenat "Luke" trobarem tots els recursos necessaris per dur a terme el treball. Per tal de poder accedir al cluster, s'ha d'introduir la següent ordre, al terminal:

$ ssh UXXXXX@luke.upf.edu (on la U representa l'identificació de cada estudiant).

En cas de voler utilitzar el luke a la FCSV per tal de poder-hi emmagatzemar el material necessari per treballar fora d'aquesta, només caldrà introduir l'anterior ordre una vegada. A partir d'aquell moment, per accedir al cluster, la següent ruta serà suficient:

/cursos/BI/bin/maphome

Cada cop que iniciem el luke, ens demanarà el nom d'usuari (la U) i la contrasenya proporcionada pels professors, que podrem canviar mitjançant l'ordre ja esmentada.

Cal tenir en compte que el cluster de docència presenta una arquitectura diferent a la dels ordinadors de la FSCV. Per tant, per poder utilitzar els diferents programes citats en les punts anteriors, caldrà introduir unes comandes lleugerament diferents a les vistes en l'apartat Exports:

$ export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH

$ cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

$ export PATH=/cursos/BI/bin:$PATH

$ export PATH=/cursos/BI/soft/exonerate/x86_64/bin:$PATH

$ export PATH=/cursos/BI/soft/t_coffee/x86_64/bin:$PATH

$ export PATH=/cursos/BI/soft/genewise/x86_64/bin:$PATH

$ export WISECONFIGDIR=/cursos/BI/soft/genewise/x86_64/wise2.2.0/wisecfg/