Automatització i programes

Per tal d'automatitzar els programes que havíem de fer servir, hem creat uns sèrie de programes en llenguatge bash i en llenguatge perl per tal de facilitar el processament d'un gran volum de dades. Concretament hem creat 6 programes bash i 1 programa perl.

Llibreria.bash: Llibreria.bash és un programa senzill que permet obtenir una llibreria de selenoproteïnes de cada espècie. Aquest transforma una base de dades amb varies querys en format fasta a format tap, per tal de distribuir cada seqüència en fitxers individuals fasta. D'aquesta manera partim d'un fitxer fasta, que ens podem descarregar per exemple de SelenoDB, i obtenim llibreria cacacterística de cada espècie que conté totes les seqüències de cada query per separat.

Blast.bash: Aquest programa és capaç de realitzar-nos alineaments locals de seqüències de tipus tBlastn de forma automàtica. Realitza el blast de totes les querys en la base de dades format fasta contra el genoma seleccionat, en aquest cas, el de Chrysemys picta bellii. Es un programa necessari per executar el següent programa JGPC_Automatitzat.bash

JGPC_Automatitzat.bash: El programa JGPC_Automatitzat.bash inclou la majoria de programes necessaris per anotar una selenoproteïna en un genoma. Aquest realitza Fastaindex, Fastafetch, Fastasubseq, Exonerate, T-coffe i Blastp contre base de dades no redundant de NCBI, tot de forma automatitzada. Per al funcionament correcte d'aquest complex programa, és necessari el petit programa perl Canvi.pl i els demès softwares esmentats abans.

Genewise.bash: El programa Genewise.bas és similar al programa anterior JGPC_Automatitzat.bash però, en ves d'utilitzar Exonerate per a realitzar la predicció d'exons utilitza el software Genewise. Així doncs obtenim una nova i, possiblement diferent, predicció de selenoproteïna en el nostre genoma. Aquest inclou Genewise, T-coffee i Blastp contra el conjunt no-redundant de NCBI. Per la seva execució es necessita el programa perl GFFtoFA.pl i els fitxers fastasubseq creats per JGPC_Automatitzat.bash.

Secis.bash: Secis.bash es un petit programa bash que realitza la cerca d'elements SECIS en la subseqüència de cada scaffold, en totes les selenoproteïnes. És necessari el programa SECISearch.pl.

SelSearch.bash: SelSearch.bash és un programa que inclou tots els passos necessaris per anotar una selenoproteïna en el nostra genoma d'interés, incloent Blast, Fastaindex, Fastafetch, Fastasubseq, predicció d'exons mitjançant Exonerate i Genewise, T-coffee i Blastp contra el conjunt no-redundant de NCBI. Realitza l'anàlisi proteïna per proteïna, fent intereccionar a l'usuari per tal d'evitar errors o bé per tal de modificar certes variables. La principal finalitat d'aquest programa ha estat tornar a analitzar aquelles selenoprotïnes que amb els programes anteriors n'obteniem resultats dubtosos o inexistents.

Per descarregar els programes clica aquí

Llibreria.bash

Aquest programa ens permet generar un directori llibreria a partir d'un d'una base de dades de querys en format fasta, que serà d'utilitat per poder consultar les querys i realitzar altres programes com T-coffe.

Requeriments previs:

El programa Llibreria.bash

Les querys en un fitxer format fasta (Ex: selenoproteinsdb_gallus.fa)
(Podeu obtenir les querys de Homo sapiens, Gallus gallus i Danio reiro en l'apartat de Resultats )

Otorgar-li permís: chmod u+x Llibreria.bash

Funcionament bàsic

El programa s'executarà un cop seleccionem la espècie per la qual volem fer la llibreria. Un cop introduïda la espècie, el programa crea el directori corresponent. A partir del fitxer fasta es defineix la variable $selenoprotein (Ex: SPPSelH_human). Transforma el fitxer fasta en format tab, per tal de poder seleccionar i separar selenoproteïnes diferents, que es tornaran a desar individualment en fitxers fasta. El resultat és un directori anomenat Llibreria_selenoproteins_$specie (Ex:Llibreria_selenoproteins_$human).

Script

#!/bin/bash

# Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
# Contacte: julia.domingo01@estudiant.upf.edu, gloria.clua01@estudiant.upf.edu

# En primer lloc, creem un fitxer amb la seqüència de cada selenoproteïna anotada a selenoproteinsdb_$specie.fa i guardem cada seqüència en
 un fitxer en format fasta en una carpeta anomenada Llibreria_selenoproteins. Cada fitxer tindrà el nom de la selenoproteïna seguint aquest
 format "SPP(nom_especie).fa". Abans però necessitem convertir el fitxer fasta en un fitxer tab. A partir daquests fitxer tab podem obtenir
 diferents fitxers fasta amb cada una de les seqüències de les selenoproteïnes.

echo
echo -e "\033[1;34mHola usuari, benvingut a el nostre programa. Abans de començar, has creat la teva Llibreria de selenoproteïnes? Quina
 espècie t'interessa analitzar? Copia el nom exacte: human , gallus o danio\033[0m";
echo

read specie

mkdir -p ./Llibreria_selenoproteins_$specie

for selenoprotein in `grep '>SPP' selenoproteinsdb_$specie.fa | cut -c2- | cut -f 1 | sort -u`;
do
  if [ ! -f ./Llibreria_selenoproteins_$specie/$selenoprotein.fa ]; then
      awk 'BEGIN{RS=">"}NR>1{sub("\n","\t"); gsub("\n",""); print RS$0}' ./selenoproteinsdb_$specie.fa > ./selenoproteinsdb_$specie.tab
      grep $selenoprotein selenoproteinsdb_$specie.tab | awk -F'\t' '{print $1"\n"$2}' > ./Llibreria_selenoproteins_$specie/$selenoprotein.fa;
  fi;

done

echo
echo -e "\033[1;34mJa has creat la teva llibreria, pots visualitzar els fitxers que es troben dins el directori ./Llibreria_selenoproteins_$sp
ecie\033[0m"
echo

Tornar a dalt

Blast.bash

Blast.bash és un programa que realitza tBlastn de forma automàtica per cada selenoproteïna del fitxer fasta de querys i genera carpetes individuals de cada selenoproteïna amb el blast desat. El mateix programa inclou les PATHs necessàries per còrrer el blast.

Requeriments previs:

El programa Blast.bash

La base de dades de querys fasta utilitzada anteriorment per la llibreria

Donar permisos per ser executat: chmod u+x Blast.bash

Funcionament bàsic

Aquest programa també es selectiu per espècie. Per tant, un cop seleccionada l'espècie que volem analitzar, crea un subdirectori blast per a cada selenoproteïna. Un cop creats, realitza dos tipus de tBlastn. Un primer normal on es mostren els alineaments, i un segon amb opció -m8 que permetrà obtenir un output en format columnar (necessari per automatitzar la presa de dades dels següents programes). Només mostra els hit que tinguin un e-value inferior a 1e-5 (l'hem considerat el nostre llindar significatiu). En cas que es volgués reduir aquesta restricció per alguna selenoproteïna, es pot realitzar mitjançant un altre programa anomenat SelSearch.bash. Finalment aquells hits significatius es desen el subdirectori creat anteriorment, específicament per a cada selenoproteïna. Aquest programa es necessari a l'hora de realitzar Exonerate i Genewise, però al ser lent (triga unes 2h a mostrar resultats), s'ha separat del programa general que inclou gran part de les comandes utilitzades.

Script

#!/bin/bash

# Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
# Contacte: gloria.clua01@estudian.upf.edu, julia.domingo01@estudiant.upf.edu

##Per a executar aquest programa cal prèviament haver creat la llibreria de selenoproteïnes (amb el programa llibreria.bash) i poder
 accedir al genoma en el qual volem anotar les selenoproteïnes (el podem trobar en /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa)

#####BLAST#####

echo
echo -e "\033[1;34mEn quina query vols que es realitzi el Blast? (Copia exactament: human , gallus o danio)\033[0m"; 
echo

read specie

##Primer hem de crear un directori per a cada selenoproteïna i un subdirectori blast per a cadascuna.

for selenoprotein  in `grep 'SPP' selenoproteinsdb_$specie.fa | cut -c2- | cut -f1 | sort -u`; do
    mkdir -p ./$selenoprotein/blast
done;

##Realitzem el blast en si.

export PATH=/cursos/BI/soft/ncbiblast/blast-2.2.13/bin/:$PATH
cp /cursos/BI/soft/ncbiblast/blast-2.2.13/bin/blastall ~/
export PATH=/cursos/BI/bin:$PATH

echo
echo -e "\033[1;34mEsperant els resultats del Blast...\033[0m";
echo

if [ ! -f ./$selenoprotein/blast/$selenoprotein.tblastn ]; then

blastall -p tblastn -i ./selenoproteinsdb_$specie.fa -d /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa -o /homes/users/U6374
4/Treball_BI/blastoutput_$specie -e 1e-5

blastall -p tblastn -i ./selenoproteinsdb_$specie.fa -d /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa -o /homes/users/U6374
4/Treball_BI/blastoutputm8_$specie -e 1e-5 -m8

fi; 

echo -e "\033[1;34mEt mostrem el resultat del blast\033[0m";

head -20 ./blastoutputm8_$specie

##Els hits estadísticament significatius de cada selenoproteina els guardem en diferents fitxers cada un dins d'un directori corresponent
 a cada selenoproteina anomenat blast:

for selenoprotein in `grep 'SPP' selenoproteinsdb_$specie.tab | cut -c2- | cut -f 1 | sort -u`; do
  for hit in `grep 'gi|' blastoutputm8_$specie | cut -f1 | sort -u`; do
     if [ "$hit" == "$selenoprotein" ]; then
        grep $hit blastoutputm8_$specie > ./$selenoprotein/blast/${hit}.tblastn;
     fi;
  done;
done;

##Els hits significatius de cada selenoproteïna han estat guardats en /$selenoprotein/blast/$hit.tblastn

Tornar a dalt

JGPC_Automatitzat.bash

Aquest programa permet anotar en el genoma que volem cada selenoproteïna inclosa en la base de dades de querys. Inclou les comandes per fer fastaindex, fastafetch, fastasubseq, exonerate, t-coffee de la predicció i cerca d'homòlegs al conjunt de proteïnes no redundants de NCBI de forma automàtica. Permet escollir a l'usuari la espècie en la qual vol realitzar l'anàlisi. A partir d'aquest punt el programa està optimitzat per realitzar totes les comandes automàticament, facilitzant a la vegada la presentació de les dades resultants per ser analitzades.

Requeriments previs:

El programa JGPC_Automatitzat.bash

El programa Canvi.pl

Accés al directori /cursos/BI/bin i /cursos/BI/soft (Inclouen programes essencials: Fastaindex, Fastafetch, Fastasubseq, Exonerate, Fastaseqfromgff.pl, Fastatranslte, T-coffe i NCBI blast)

Els resultats en format .tblastn del programa Blast.bash

La llibreria amb totes les selenoproteïnes que utilitzem com a query

Donar permisos: chmod u+x JGPC_Automatitzat.bash

Donar permisos: chmod u+x Canvi.pl

Funcionament bàsic

En primer lloc, un cop donades les PATHs, el programa procedeix a generar un index del nostre genoma. En aquest cas l'ordena per scaffolds. Un cop fet el fastaindex, gràcies a la variable $selenoprotein selecciona una selenoproteïna i realitza tots els programes que inclou JGPC_Automatitzat.bash en aquesta. Un cop acaba torna a començar amb la següent. A més a més dins de cada selenoproteïna, el blast ens pot haver trobat més d'un hit significatiu en diferents scaffolds. Ja que aquests scaffolds es troben en diferents regions del genoma, vam voler analitzar els hits de cada scaffold per separat (per poder distinguir selenoproteïnes amb graus d'homologia elevats, dominis conservats o bé, duplicacions). És per això que vam introduir la variable $scaffold. Per tant un cop s'inicia el programa i se sel·lecciona la espècie d'origen de les querys, es realitza el fastafetch, fastasubseq, exonerate, t-coffee i Blast NCBI non-redundant els hits significatius de cada scaffold de cada selenoproteïna. Remarcar que, previ a la realització del t-coffe era necessari canviar la U de la query per una X. Vam crear un petit programa perl per solucionar aquesta limitació i eliminar qualsevol altre símbol que Exonerate no reconeguès.

Script

#!/bin/bash

## Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
## Contacte: julia.domingo01@estudiant.upf.edu, gloria.clua01@estudiant.upf.edu

echo -e "\033[1;34mHola usuari, a continuació iniciarem el programa JGPC_Automatitzat! D'aquesta manera es podran obtenir anotacions de
 les selenoproteïnes actualment descrites en el genoma de Chrysemys picta belli. Per quina espècie vols que es realitzi el programa
 automatitzat?\033[0m";

echo -e "\033[1;34mCopia exactament: human , gallus o danio\033[0m";
echo

read specie

###El primer pas consisteix en l'extracció de la regió genòmica que potencialment conté el gen que estem buscant. Amb aquest objectiu
 farem anar uns programes que formen part d'un software anomenat exonerate. Cal tenir instal·lats els següents programes; fastaseqfromGFF.pl
 i canvi.pl. També has de tenir el resultat del blast en la carpeta "Llibreria_selenoproteins" i una carpeta generada amb el nom de cada
 selenoproteïna on es guardaran tots els documents resultants en córrer el programa separats en diferents carpetes internes.A més cal poder
 accedir al genoma de c.picta que està guardat en /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa.  

export PATH=/cursos/BI/bin:$PATH
export PATH=/cursos/BI/soft/exonerate/i386/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH

### Fastaindex: et crea el document genome.index que serà la base de dades que requereix el fastafetch pel seu funcionament.

if [ ! -f ./genome.index ]; then
	fastaindex /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa ./genome.index;
fi;

#####EXONERATE#####

echo -e "\033[1;34mExonerate\033[0m";
echo -e "\033[1mA partir del document blastoutput_$specie, s'estimen les posicions cromosòmiques on comença i acaba la regió sencera que
 ocupa aquest gen i s'extreu la regió cromosòmica on està en un fitxer FASTA.\033[0m";

### Inici del programa

for selenoprotein in `grep 'SPP' ./selenoproteinsdb_$specie.tab | cut -c2- | cut -f1 | sort -u`; do
for scaffold in `grep 'gi|' ./$selenoprotein/blast/$selenoprotein.tblastn | cut -f2 |  sort -u`; do

	echo "$selenoprotein, $scaffold";
	echo

	scaf=`echo $scaffold|cut -f4 -d"|"`

### Programes per executar exonerate

### Fastafetch: et localitza la regió cromosòmica dins el genoma de C.picta on es localitzen els hits amb una valoració positiva.

	mkdir -p ./$selenoprotein/fastafetch/
	fastafetch /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa genome.index $scaffold > "./$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db"

	echo -e "\033[1;34mS'ha obtingut el resultat de fastafetch\033[0m";

### Fastasubseq: et permet localitzar i extendre la regió genòmica que possiblement codifica per la proteïna cercada.

mkdir -p ./$selenoprotein/fastasubseq/
        column9=`grep $scaf ~/Treball_BI/$selenoprotein/blast/$selenoprotein.tblastn | cut -f9 | sort -n | sed -n 1p`;
        column10=`grep $scaf ~/Treball_BI/$selenoprotein/blast/$selenoprotein.tblastn | cut -f10 | sort -n | sed -n 1p`;

        if [ "$column9" -lt "$column10" ];
            then
                start=$(($column9 - 20000));
		if [ $start -le 0 ]
                then
                    $start=0;
                else
                    end=`grep $scaf ~/Treball_BI/$selenoprotein/blast/$selenoprotein.tblastn | cut -f10 | sort -r | sed -n 1p`;
                    length=$(($end - $start + 20000));
                fi;

        else
                start=$(($column10 - 20000));
		if [ $start -le 0 ]
                then
                    $start=0;
                else
                    end=`grep $scaf ~/Treball_BI/$selenoprotein/blast/$selenoprotein.tblastn | cut -f9 | sort -r | sed -n 1p`;
                    length=$(($end - $start + 20000));
		    fi;

        fi;

	fastasubseq ./$selenoprotein/fastafetch/$selenoprotein.$scaf.db $start $length > ./$selenoprotein/fastasubseq/$selenoprotein.$scaf.subseq.fa 

###Exonerate

	mkdir -p ./$selenoprotein/exonerate

        ###Per si al fasta de sequencia hi ha Us, cal substituir-les per X, però només un cop!

        if [ "$int_UX"=="0" ]; then
        ./canvi.pl < ./Llibreria_selenoproteins_$specie/"${selenoprotein}.fa" > ./$selenoprotein/canviXU.fa;
        int_UX="1";
        fi;

	exonerate -m p2g --showtargetgff -q "./$selenoprotein/canviXU.fa" -t ./$selenoprotein/fastasubseq/"${selenoprotein}.${scaf}.subseq.fa" > ./$selenoprotein/exonerate/"${selenoprotein}.${scaf}.exonerate.gff"

	###El programa fastaseqfromGFF.pl el farem anar en dos passos, en un primer pas obtindrem en un fitxer apart
	###l'anotació en GFF del programa exonerate, i en un segon pas extreurem la seqüència.

	echo
	echo -e "\033[1;34mS'ha obtingut el fitxer gff\033[0m";
	echo

###FastaseqfromGFF.pl (extraccio de cDNA)

###Fastatranslate, de cDNA a proteina

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

        echo
        echo -e "\033[1;34mS'esta duent a terme el programa Fastatranslate que ens traduira el cDNA a proteina\033[0m";
        echo

    ##Per extreure el cDNA i traduir-lo a proteïnes

		echo -e "\033[1;34mPreparant i executant per t_coffee:\033[0m";

		exon=`grep -w exon ./$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff`
		echo "$exon" > ./$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.exon.gff
                if [ "$exon" != "" ]; then

                        mkdir -p ./$selenoprotein/cDNA

                        fastaseqfromGFF.pl ./${selenoprotein}/fastasubseq/"${selenoprotein}.${scaf}.subseq.fa" ./${selenoprotein}/exonerate/"${selenoprotein}.${scaf}.exonerate.exon.gff" > ./${selenoprotein}/cDNA/${selenoprotein}.${scaf}.cDNA

                        mkdir -p ./$selenoprotein/proteina
                        fastatranslate -F 1 "./$selenoprotein/cDNA/${selenoprotein}.${scaf}.cDNA" > "./$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa"

                fi

	echo

#####TCOFFEE#####

	mkdir -p ./$selenoprotein/t_coffee

	t_coffee "./Llibreria_selenoproteins_$specie/${selenoprotein}.fa" "./$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa" > ./$selenoprotein/t_coffee/"${selenoprotein}.${scaf}.aln"

mv *html ./$selenoprotein/t_coffee
mv *aln ./$selenoprotein/t_coffee

#####NCBI non-redundant#####

echo -e "\033[1:34m Cerca de proteïnes similars en la base de dades de proteïnes no redundant d'NCBI\033[0m"

mkdir -p ./$selenoprotein/NCBI

export PATH=/cursos/BI/bin/netblast/bin:$PATH
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

blastall -p blastp -i "./$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa" -d /cursos/BI/soft/selenoprofiles_2_installation/libraries/nr -o ./$selenoprotein/NCBI/${selenoprotein}.${scaf}.ncbi 

   done
done

Canvi.pl

#!/usr/bin/perl
use strict;

my @v;
my $string;
my $counter = 0;

while ()
{
    if ($counter == 0)
    {
        print $_;
    }
    if ($counter > 0)
    {
        @v=split(//,$_);

        while (<@v>)
        {
            $_=~ tr/\$#%@~!&*()[];.,:?^ `\\\///d;
	    if ($_ eq "U")
            {
                $string.= "X";
            }
            else
            {
                $string.= $_;
            }
	    if ($_ eq "@")
	    {
		$string.= " ";
	    }
        }
    }
    $counter++;
    $string.="\n";
}

Tornar a dalt

Genewise.bash

Genewise.bash ens permet tornar anotar en el nostre genoma cada selenoprotina, partint base de querys existent, mitjançant el software Genewise. És un programa lent (acaba de córrer entre 1 i 3 hores segons el tamany de la llibreria), raó per la qual no es troba integrat en el programa JGPC_Automatitzat.bash

Requeriments previs:

El programa Genewise.bash

El programa GFFtoFA.pl

Accés al directori /cursos/BI/bin i /cursos/BI/soft

Els tBlastn generats per JGPC_Automatitzat.bash

La llibreria amb les querys

Permisos: chmod u+x Genewise.bash

Permisos: chmod u+x GFFtoFA.pl

Funcionament bàsic

El programa s'inicia i ens permet triar de quina procedència volem que siguin les nostres querys. Un cop selecionada la espècie, defineix la variable $selenoprotein i $scaffold (de la mateixa manera que ho fa JGPC_Automatitzat.bash). Gràcies a això, aconseguim realitzar genewise, t-coffee de la proteïna predita per aquest úlim i Blastp NCBI non-redundant els hits significatius de cada scaffold de cada selenoproteïna. Aquest programa incorpora l'us del programa GFFtoFA.pl. Aquest es va modificar a partir del programa creat pel grup 4c del curs 2011-2012. La seva finalitat és transformar el fitxer creat per genewise en format gff a format fasta per poder ser utilitzat per T-coffee.

Script

# Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
# Contacte: gloria.clua01@estudian.upf.edu, julia.domingo01@estudiant.upf.edu

#!/bin/bash

export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

echo -e "\033[1;34mHola usuari, a continuació iniciarem el programa genewise.bash! D'aquesta manera es podran obtenir anotacions de les
 selenoproteïnes actualment descrites en el genoma de Chrysemys picta belli. Per quina espeie  vols que es realitzi el programa?
 (Copia exactament: human , gallus o danio)\033[0m";
echo

read specie

for selenoprotein in `grep 'SPP' ./selenoproteinsdb_$specie.tab | cut -c2- | cut -f1 | sort -u`; do
for scaffold in `grep 'gi|' ./$selenoprotein/blast/$selenoprotein.tblastn | cut -f2 |  sort -u`; do

	    scaf=`echo $scaffold|cut -f4 -d"|"`

#####GENEWISE#####

##Fem un genewise tenint en compte si es una seqüencia és forward o reverse)

        mkdir -p ./$selenoprotein/genewise

        column9=`grep $scaf ~/Treball_BI/${selenoprotein}/blast/${selenoprotein}.tblastn | cut -f9 | sort -n | sed -n 1p`;
        column10=`grep $scaf ~/Treball_BI/${selenoprotein}/blast/${selenoprotein}.tblastn | cut -f10 | sort -n | sed -n 1p`;

      if [ "$column9" -lt "$column10" ]; then

	    genewise -pep -pretty -cdna -gff "./$selenoprotein/canviXU.fa" ./$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa >
./$selenoprotein/genewise/${selenoprotein}.${scaf}.genewise.gff

	else
	
	    genewise -pep -pretty -cdna -gff "./$selenoprotein/canviXU.fa" ./$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa -trev > ./$selenoprotein/genewise/${selenoprotein}.${scaf}.genewise.gff
	    
	fi;
	
mkdir -p ./$selenoprotein/genewise_proteina
echo
echo -e "\033[1;34mS'està obtenint la seqüència peptídica a partir del programa GFFtoFA.pl\033[0m";

./GFFtoFA.pl < "./$selenoprotein/genewise/${selenoprotein}.${scaf}.genewise.gff" > "./$selenoprotein/genewise_proteina/${selenoprotein}.${scaf}.genewise.fa"


#####T_COFFEE del genewise######

echo
echo -e "\033[1;34mRealitzant t_coffee\033[0m";
echo

mkdir ./$selenoprotein/t_coffee_genewise

t_coffee "./Llibreria_selenoproteins_$specie/${selenoprotein}.fa" ./$selenoprotein/genewise_proteina/${selenoprotein}.${scaf}.genewise.fa >./$selenoprotein/t_coffee_genewise/${selenoprotein}.${scaf}.genewise.tcoffee.aln

mv *html ./$selenoprotein/t_coffee_genewise
mv *aln ./$selenoprotein/t_coffee_genewise

echo -e "\033[1;34mRealitzant NCBI nr\033[0m";

###Busquem la proteïna nova en la base de dades de proteïnes no redundant de l'NCBI

mkdir ./$selenoprotein/NCBI_genewise

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

blastall -p blastp -i "./$selenoprotein/genewise_proteina/${selenoprotein}.${scaf}.genewise.fa" -d /cursos/BI/soft/selenoprofiles_2_installation/libraries/nr -o ./$selenoprotein/NCBI_genewise/${selenoprotein}.${scaf}.ncbi

echo
echo -e "\033[1mJa s'ha obtingut el resultat del genewise, t_coffee del genewise i NCBI del genewise\033[0m"
echo 

   done
done

GFFtoFA.pl

#!/usr/bin/perl -w 

use strict;
my $i;

$/=">gi";
$i = 0;

print ">gi";
while (){
    if ($i==1){
	chomp;
	$_=~ s/\///g; 
	print;
    }
    $i = $i + 1;
}

Tornar a dalt

Secis.bash

Secis.bash realitza la busqueda d'elements SECIS en cada regió genòmica seleccionada prèviament. Ho fa mitjançant l'execució del programa SECISearch.pl, disponible internet o al directori /cursos/BI/bin.

Requeriments previs:

El programa Secis.bash

Descarregar el programa SECISearch.pl (també disponible al directori /cursos/BI/bin)

Els fitxers del fastasubseq creats per JGPC_Automatitzat.bash

Donar permisos: chmod u+x Secis.bash

Donar permisos: chmod u+x SECISearch.pl

Funcionament bàsic

Secis.bash bàsicament permet l'execcució del programa SECISerch.pl. El punt clau recau en definir les variables $selenoprotein i $scaffold, ja que un cop definides, realitzarà la cerca d'elements secis en tot el subseq de cada scaffold determinat en cada una de les selenoproteïnes.

Script

# Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
# Contacte: julia.domingo01@estudiant.upf.edu, gloria.clua01@estudiant.upf.edu

#!/bin/bash

echo -e "\033[1;34mHola usuari, a continuació iniciarem el programa SECIS! D'aquesta manera es podran obtenir les seqüències de les regions
 SECIS trobades en el genoma. Per quina espècie vols cercar els SECIS?\033[0m";

echo -e "\033[1;34mCopia exactament: human , gallus o danio\033[0m";
echo

read specie

for selenoprotein  in `grep 'SPP' selenoproteinsdb_$specie.fa | cut -c2- | cut -f1 | sort -u`; do
	for scaffold in `grep 'gi|' ./$selenoprotein/blast/$selenoprotein.tblastn | cut -f2 |  sort -u`; do

		scaf=`echo $scaffold|cut -f4 -d"|"`

mkdir -p ./$selenoprotein/secis

./SECISearch.pl < ./$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa > ./$selenoprotein/secis/${selenoprotein}.${scaf}.html

	done
done

Tornar a dalt

SelSearch.bash

Aquest programa automatitza tot el protocol descrit a materials i mètodes permetent a l'usuari manipular els paràmetres manualment en els punts crítics evitant possibles errors del procés d'automatització. Fa tBLASTn, fastaindex, fastafetch, fastasubseq, exonerate, genewise, t-coffee d'ambdues prediccions i cerca d'homòlegs al conjunt de proteïnes no redundants de NCBI.

Requeriments previs:

El programa SelSearch.bash

El programa Canvi.pl

El programa GFFtoFA.pl

Les bases de querys en format fasta

Les llibreries de selenoproteïnes (Executar previment Llibreria.bash)

Accés al directori /cursos/BI/bin

Executar els permisos per els programes: chmod u+x SelSearch, chmod u+x Canvi.pl i chmod u+x GFFtoFA.pl

Script

#!/bin/bash

# Programa fet pel grup 3a de 4t de biologia del curs 2012-2013  de la Universitat Pompeu Fabra
# Contacte: gloria.clua01@estudiant.upf.edu, julia.domingo01@estudiant.upf.edu

#####SelSearch#####

##Consisteix en anotar una selenoproteïna de manera maunal però concatenant tots els programes requerits en un sol programa, d'aquesta
 manera es podran repetir o repassar aquelles selenoproteïnes que hagin sortit subtoses o calgui retocar la seva cerca.

echo
echo -e "\033[1;34mHola usuari, a continuació realitzarem el programa SelSearch. L'objectiu és anotar una selenoproteïna de manera manual
en el genoma de Chrysemys picta bellii. Et recomanem que obris una altra pestanya amb el Luke obert per realitzar alguna modificació 
manual si fos convenient. Doncs comencem!\033[0m"

### Fastaindex

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

### Fastaindex: et crea el document genome.index que serà la base de dades que requereix el fastafetch pel seu funcionament.

if [ ! -f ./genome.index ]; then
        fastaindex /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa ./genome.index;
fi;

### Fastafetch d'una selenoproteïna d'interès

### Escullim la selenoproteïna que volem cercar en el nostre organisme

echo -e "\033[1mQuina selenoproteïna volem cercar en el nostre genoma? Primer indica'ns l'espècie que t'interessa (copia exactament una de
les tres)\033[0m"
echo
echo -e "\033[1mhuman gallus danio\033[0m"
echo
read specie
echo

echo -e "\033[1mQuina selenoproteïna vols cercar exactament?\033[0m"
ls ./Llibreria_selenoproteins_$specie/
echo

echo -e "\033[1mIntrodueix el nombre exacte de la selenoproteïna (copia el nom que comença per "SPP" i l'enganxes directament a continuació):\033[0m"
read selenoprotein

mkdir -p ./SelSearch/$selenoprotein/blast

##Realitzem el blast en si.

export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/
export PATH=/cursos/BI/bin:$PATH

echo
echo -e "\033[1;34mEsperant els resultats del Blast...requereix aproximadament uns 3 minuts\033[0m";
echo

echo -e "\033[1mQuin e-value consideres significatiu?. Copia directament: 1e-5 , 1e-4 , 1e-3 , 1e-2 , 1e-1\033[0m";
read evalue

echo
echo 

blastall -p tblastn -i ./Llibreria_selenoproteins_$specie/$selenoprotein.fa -d /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa -o /homes/users/U63744/Treball_BI/SelSearch/$selenoprotein/blast/blastoutput_${selenoprotein}_$specie -e $evalue

blastall -p tblastn -i ./Llibreria_selenoproteins_$specie/$selenoprotein.fa -d /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa -o /homes/users/U63744/Treball_BI/SelSearch/$selenoprotein/blast/blastoutputm8_${selenoprotein}_$specie -e ${evalue} -m8

echo
echo -e "\033[1;34mEt mostrem el resultat del blast\033[0m";
echo -e "\033[1m Quin d'aquests hits resultants del Blast creus que es correspon a la selenoproteïna que busques? \033[0m"
echo -e "\033[1m Copia directament la regió escollida, amb aquest format gi|...|gb|...| :  \033[0m"
echo
more ~/Treball_BI/SelSearch/$selenoprotein/blast/blastoutputm8_${selenoprotein}_$specie
echo
echo
read scaffold
echo
echo

scaf=`echo $scaffold|cut -f4 -d"|"`

mkdir -p ./SelSearch/$selenoprotein/fastafetch/

fastafetch /cursos/BI/genomes/project_2013/Chrysemys_picta_bellii/genome.fa genome.index $scaffold > ./SelSearch/$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db

###Fastasubseq, per realitzar el fastasubseq és necessari modificar la regió genòmica. A partir d'aquestes posicions cromosòmiques, 
estimeu les posicions cromosòmiques on comença i acaba la regió sencera que ocupa aquest gen i feu-les anar per extreure la regió
cromosomica on està en un fitxer FASTA apart mitjançant el programa fastasubseq (què es part dels programes que acompanyen exonerate)

echo -e "\033[1mAra posarem en marxa el programa Fastasubseq ...\033[0m"
echo
more ~/Treball_BI/SelSearch/$selenoprotein/blast/blastoutputm8_${selenoprotein}_$specie
echo
echo
echo -e "\033[1mPrimer localitza la posició cromosòmica on comença la regió que ocupa aquest gen. Per poder analitzar la seqüència acuradament
cal que extenguis la posició cromosòmica inicial uns 20000 nucleòtids. D'aquesta manera t'assegures que estàs agafant tot el gen. 
Si l'inici es troba pròxim al començament de l'scaffold (a una llargada menor a 20000nt) posa 1 com a inici de la seqüència :  \033[0m"
echo
read inici
echo
echo -e "\033[1mIntrodueix la longitut de la seqüència, extén uns 20000 nucleotids la posició cromosòmica final tal i com has fet en el 
pas anterior:  \033[0m"
echo
read longitut
echo

mkdir -p ./SelSearch/$selenoprotein/fastasubseq/

fastasubseq ./SelSearch/$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db $inici $longitut > ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa 2> error.temp

error=`cat error.temp | egrep "Subsequence must end before end of"`;
if [ "$error" != "" ]; then

echo -e "\033[1mLa longitut indicada és més gran que la seqüència de l'scaffold, en aquest cas el programa detecta un error. Cal que redueixis 10000 nucleòtids. Torna a indicar la longitut:\033[0m"
        echo
        read longitut
        echo
fastasubseq ./SelSearch/$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db $inici $longitut > ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa 2> error1.temp
fi;

error1=`cat error1.temp | egrep "Subsequence must end before end of"`;
if [ "$error1" != "" ]; then
echo -e "\033[1mLa longitut indicada continua essent més gran que la seqüència de l'scaffold, en aquest cas el programa detecta un error.
Cal que tornis a indicar la longitut però considerant l'últim nucleòtid mostrat pel blast com a final de seqüència:\033[0m"
	echo
	read longitut
	echo

	fastasubseq ./SelSearch/$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db $inici $longitut > ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa
fi;

echo -e "\033[1mEt mostrem la seqüència seleccionada. Si estàs d'acord pressiona Q per sortir del more\033[0m"

head ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa
echo
echo
read -p "A continuació tecleja [Enter] per continuar"

echo

###A partir de la selenoproteïna escollida i de la regió genòmica que hem extret anteriorment generarem una anotació del gen que dóna lloc a
aquesta proteina mitjançant el programa exonerate. El directori /exonerate guardat en /cursos/BI/bin conté el programa fastaseqfromGFF.pl,
essencial per poder realitzar el programa Exonerate.\033[0m" 

###Exonerate

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

mkdir -p ./SelSearch/$selenoprotein/exonerate

echo -e "\033[1;34mExonerate\033[0m"
echo

./canvi.pl < ./Llibreria_selenoproteins_$specie/"${selenoprotein}.fa" > ./SelSearch/$selenoprotein/canviXU.fa;

exonerate -m p2g --exhaustive yes --showtargetgff -q ./SelSearch/$selenoprotein/canviXU.fa -t ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa > ./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff

more  ./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff

read -p "Si no estàs content amb el resultat obtingut bé perquè no t'ha alineat la selenocisteïna o perquè necessites allargar la seqüència,
pots modificar l'inici i la longitud de la seqüència per realitzar un altre Exonerate un cop més. Vols realitxar aquesta modificació? Copia 
SI (en majúscules) en cas afirmatiu si no, copia NO (en majúscules) per continuar;    "

if [[ "$REPLY" == "SI" ]]; then
echo -e "\033[1mObre el resultat de la taula del blast ./SelSearch/$selenoprotein/blast/blastoutputm8_$selenoprotein_$specie i allarga 
l'inici i la longitud per a poder tornar a realitzar l'Exonerate. Indica'ns el nou inici: \033[0m" 
read inici
echo -e "\033[1m Indica'ns la nova longitut desitjada: \033[0m"
read longitut

fastasubseq ./SelSearch/$selenoprotein/fastafetch/${selenoprotein}.${scaf}.db $inici $longitut > ./SelSearch/$selenoprotein/fastasubseq/$selenoprotein.${scaffold}.subseq.fa

echo -en "\033[1m Et mostrem la seqüència seleccionada. Si estàs d'acord pressiona Q per sortir del more\033[0m"
echo
read -p "Tecleja [Enter] per conitnuar"
more  ./SelSearch/$selenoprotein/fastasubseq/$selenoprotein.$scaffold.subseq.fa
echo

exonerate -m p2g --exhaustive yes --showtargetgff -q ./SelSearch/$selenoprotein/canviXU.fa -t ./SelSearch/$selenoprotein/fastasubseq/${selenoprotein}.${scaf}.subseq.fa > ./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff

more  ./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff

fi;
echo
echo
echo -e "\033[1m A continuació guardarem els exons en un document per separat anomenat ${selenoprotein}.${scaf}.exonerate.exon.gff. El 
trobaràs dins el directori ./SelSearch/exonerate\033[0m"

echo
echo -e "\033[1mObtenint el cDNA\033[0m"
echo

###El programa fastaseqfromGFF.pl el farem anar en dos passos, en un primer pas obtindrem en un fitxer apart l'anotació en GFF del programa 
exonerate, i en un segon pas extreurem la seqüència.

###FastaseqfromGFF.pl (extraccio de cDNA)

mkdir -p ./SelSearch/$selenoprotein/cDNA

export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

exon=`grep -w exon ./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff`
echo "$exon" > "./SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.exon.gff";

if [ "$exon" != "" ]; then

	fastaseqfromGFF.pl ./SelSearch/${selenoprotein}/fastasubseq/"${selenoprotein}.${scaf}.subseq.fa" ./SelSearch/${selenoprotein}/exonerate/"${selenoprotein}.${scaf}.exonerate.exon.gff" > ./SelSearch/${selenoprotein}/cDNA/${selenoprotein}.${scaf}.cDNA

fi;

###Fastatranslate, de cDNA a proteina

echo
echo -e "\033[1mS'està duent a terme el programa Fastatranslate que ens traduirà el cDNA a proteïna\033[0m"
echo

mkdir -p ./SelSearch/$selenoprotein/proteina

fastatranslate -F 1 "./SelSearch/$selenoprotein/cDNA/${selenoprotein}.${scaf}.cDNA" > "./SelSearch/$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa"
echo

echo -e "\033[1mS'ha creat un fitxer amb la seqüència aminoacídica de la proteïna predita a 
"./SelSearch/proteina/${selenoprotein}.${scaf}.traslate.fa". Obre l'arxiu en la pestanya que prèviament has obert i escull un dels 6 marcs 
de lectura que et mostra el document. Deixa'l en format fasta i guarda teclejant C-X C-S i per sortir C-X C-C. Si tens dubtes del patró a 
escollir pots mirar el resultat de l'exonerate /SelSearch/$selenoprotein/exonerate/${selenoprotein}.${scaf}.exonerate.gff.\033[0m"

###T_coffee

echo
echo -e "\033[1mRealitzant t_coffee. Si t'indica que escriguis en email pots teclejar Enter per continuar.\033[0m"
echo

mkdir -p  ./SelSearch/$selenoprotein/t_coffee

t_coffee ./Llibreria_selenoproteins_$specie/${selenoprotein}.fa ./SelSearch/$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa > ./SelSearch/$selenoprotein/t_coffee/"${selenoprotein}.${scaf}.aln"

mv *html ./SelSearch/$selenoprotein/t_coffee 
mv *aln ./SelSearch/$selenoprotein/t_coffee  

###Hauríeu de poder observar que l'anotació generada anteriorment dóna lloc exactament a la mateixa proteïna a partir de la qual hem buscat 
l'anotació del gen que la codifica. Observeu també la forma en què TCOFFEE mostra l'alineament de la selenocisteina.

###Comparar el teu hit amb la NRdatabase per comporvar si existeix alguna proteina coneguda codificada dins els teu hit

echo -e "\033[1mComparant el hit resultant amb la base de dades no redundant d'NCBI.\033[0m"

mkdir -p ./SelSearch/$selenoprotein/NCBIex

export PATH=/cursos/BI/bin/netblast/bin:$PATH
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

blastall -p blastp -i "./SelSearch/$selenoprotein/proteina/${selenoprotein}.${scaf}.translate.fa" -d /cursos/BI/soft/selenoprofiles_2_installation/libraries/nr -o ./SelSearch/$selenoprotein/NCBIex/${selenoprotein}.${scaf}.ncbi

###Anàlogament podem realitzar el programa Genewise per confirmar els nostres resultats. A partir de la selenoproteina escollida i de la 
regio genòmica que hem extret anteriorment del resultat del blast, generarem una anotació del gen que dóna lloc a aquesta proteïna 
mitjançant el programa GeneWise.

echo -e "\033[1mVols realitzar el programa Genewise? Aquest et permetrà comporvar els resultats obtinguts amb Exonerat per confirmar la teva
cerca. Indica SI o NO (en majúsucles)    \033[0m" 
read resposta

if [ "$resposta" = "SI" ]; then
echo
echo -e "\033[1mRealitzant el Genewise si la teva seqüència és invertida, cal que escriguis -trev a continuació. Si és directa tecleja Enter
per continuar  \033[0m"
read trev
echo

mkdir -p ./SelSearch/$selenoprotein/genewise

export PATH=/cursos/BI/bin:$PATH
export WISECONFIGDIR=/cursos/BI/soft/wise-2.2.0/wisecfg

echo
genewise -pep -pretty -cdna -gff ./Llibreria_selenoproteins_$specie/${selenoprotein}.fa ./SelSearch/${selenoprotein}/fastasubseq/${selenoprotein}.${scaf}.subseq.fa $trev > ./SelSearch/$selenoprotein/genewise/${selenoprotein}.${scaf}.genewise.gff

echo
echo -e "\033[1mObtenint seqüència peptídica a partir del fitxar obtingut en format .gff. Cal tenir el programa GFFtoFA.pl instal·lat en el
directori ./Treball_BI. Si ja teniu el programa guardat en la carpeta Treball_BI podeu prémer a l'enter directament\033[0m"
read si

chmod u+x ./GFFtoFA.pl

mkdir -p ./SelSearch/$selenoprotein/genewise_proteina/

./GFFtoFA.pl < ./SelSearch/$selenoprotein/genewise/$selenoprotein.${scaf}.genewise.gff > ./SelSearch/$selenoprotein/genewise_proteina/$selenoprotein.${scaf}.genewise.proteina.fa

###Igual que amb el resultat de l'Exonerate, també realitzem un Tcoffee del resultat de Genewise

echo
echo -e "\033[1mRealitzant el t_coffee\033[0m"
echo
mkdir -p ./SelSearch/$selenoprotein/genewise_TCOFFEE

t_coffee ./Llibreria_selenoproteins_$specie/${selenoprotein}.fa ./SelSearch/${selenoprotein}/genewise_proteina/${selenoprotein}.${scaf}.genewise.proteina.fa > ./SelSearch/$selenoprotein/genewise_TCOFFEE/"${selenoprotein}.${scaf}.aln"

mv *html ./SelSearch/$selenoprotein/genewise_proteina/
mv *aln ./SelSearch/$selenoprotein/genewise_proteina/

fi;

###NCBI_genewise

###Comparar el teu hit amb la NRdatabase per comporvar si existeix alguna proteina coneguda codificada dins els teu hit

echo -e "\033[1mComparant el hit resultant amb la base de dades no redundant d'NCBI.\033[0m"

mkdir -p ./SelSearch/$selenoprotein/NCBI

export PATH=/cursos/BI/bin/netblast/bin:$PATH
export PATH=/cursos/BI/bin/ncbiblast/bin:$PATH
cp /cursos/BI/bin/ncbiblast/.ncbirc ~/

blastall -p blastp -i "./SelSearch/$selenoprotein/genewise_proteina/${selenoprotein}.${scaf}.genewise.proteina.fa" -d /cursos/BI/soft/selenoprofiles_2_installation/libraries/nr -o ./SelSearch/$selenoprotein/NCBI/${selenoprotein}.${scaf}.ncbi

echo
echo -e "\033[1mFI.Esperem que que el programa hagi funcionat correctament i hagis pogut anotar la selenoproteïna en el genoma de C.picta bellii\033[0m"
echo 

echo
echo -e "\033[1mIMPORTANT! Recorda que pots guardar tots els resultats obtinguts dins el directori de la selenoproteina anotada en un altre
directori de l'ordinador. T'aconsellem que copiis directament el directori en un directori extern en el shell personal de l'usuari que
s'anomeni $selenoprotein_$scafold_$especie. Un cop et trobis dins el directori, pots copiar els resultats utilitzant aquesta comanda desde
EL SHELL DEL VOSTRE ORDINADOR:( $scp U63739@luke.upf.edu:/homes/users/UXXXXX/Treball_BI/SelSearch/$selenoprotein .) Moltes gràcies!\033[0m"
echo

Tornar a dalt

Recerca de selenoproteïnes en Chrysemys picta bellii

Automatització i programes

Llibreria.bash

Blast.bash

JGPC_Automatitzat.bash

Canvi.pl

Genewise.bash

GFFtoFA.pl

Secis.bash

SelSearch.bash

Recerca de selenoproteïnes
en Chrysemys picta bellii