Programa


Objetivo

En este trabajo, hemos querido elaborar un programa que fuese capaz de mostrar el nivel de conservación de tránscritos de genes ortólogos de humano y ratón. Éste es fácilmente aplicable a otras especies, la única diferencia consistiría en aplicar distintos ficheros de base.

Este objetivo puede lograrse de dos maneras distintas, en función de las necesidades:

  1. analizar un fichero que contenga genes o tránscritos de una misma especie obteniendo como resultado final una relació,n de genes ortólogos, tránscritos de los genes ortólogos, alineamiento de tránscritos y valor de dicho alineamiento pudiendo sacar conclusiones referentes a su conservación
  2. analizar la conservación de un tránscrito concreto en la otra especie, de manera que podrímos obtener su gen correspondiente, el gen ortólogo, y los tránscritos del gen ortólogo con el correspondiente alineamiento y valores del alineamiento de los tránscritos

Para la elaboración de este trabajo, podemos distinguir dos partes bien diferenciadas: obtención y manipulación de archivos y ficheros, y diseño de algorismos para construir un programa capaz de trabajar con dichos archivos.

A continuación explicaremos con detenimiento las distintas partes de las que consta la elaboración de este trabajo.

Ficheros

Todos los ficheros que hemos utilizado para la confección de este programa han sido obtenidos a partir de EnsMart, un servicio de Ensembl, que nos permite obtener identificadores de genes, tránscritos, péptidos y proteínas, además de sus secuencias y/o estructuras entre otras informaciones. De este modo, y mediante la selección del formato adecuado, podemos conseguir los distintos archivos de texto, que podrán ser modificados posteriormente en función de la información que se requiera en cada momento.

El programa no está diseñado exclusivamente para poder procesar informació procedente de Ensembl, ya que podemos elaborar ficheros de base a partir de otros servidores o bases de datos y utilizarlos de este modo con nuestro programa. Aun así sí es imprescindible contar con los siguientes ficheros para obtener resultados a partir de este programa:

  1. un fichero que contenga los genes ortólogos de las dos especies a estudiar en formato fasta o texto
  2. otro fichero con la correspondencia de genes y tránscritos también en formato fasta o texto
  3. un fichero para cada una de las especies en formato gff (éste contendrá identificadores para tránscrito, gen y péptido, estructura exónica y secuencia peptídica entre otra información)


Ortólogos

Tras obtener todos los genes en humano con ortólogos en ratón y viceversa tal como hemos descrito en el punto anterior, hemos modificado ambos ficheros mediante comandos de UNIX de manera que nos quedara un solo fichero que consta de dos columnas: en una se muestran los genes en humano, y en la otra los ortólogos en ratón (este procedimiento se puede llevar a cabo para dos especies distintas, sin necesidad de que sean las que hemos utilizado para la realización de este trabajo) 

ENSG00000067646.1       ENSMUSG00000053211.1
ENSG00000169953.2       ENSMUSG00000045336.1
ENSG00000185275.1       ENSMUSG00000047139.1
ENSG00000184895.1       ENSMUSG00000043876.2
ENSG00000172468.2       ENSMUSG00000045336.1
ENSG00000156265.2       ENSMUSG00000025610.2
ENSG00000185808.1       ENSMUSG00000022940.2
ENSG00000171587.4       ENSMUSG00000050272.2
ENSG00000171587.4       ENSMUSG00000040619.1
ENSG00000171587.4       ENSMUSG00000051090.2
ENSG00000171587.4       ENSMUSG00000051267.1
ENSG00000182670.1       ENSMUSG00000040785.2
ENSG00000142173.3       ENSMUSG00000020241.1
Ver fichero

Tránscritos

A partir del fichero de genes ortólogos, separamos las dos columnas en dos ficheros distintos: uno para los genes de humano y otro para los genes de ratón. Una vez hecho esto, volveremos a utilizar EnsMart para obtener los tránscritos correspondientes a los distintos genes. De este modo podremos llegar a obtener dos ficheros distintos, en los que encontramos en una columna los tránscritos y en la otra los genes, uno para cada especie. Una vez tenemos los dos ficheros, los unificamos de manera que queden como a continuación: 

ENST00000264955.1 ENSG00000001952.1      
ENST00000226299.2 ENSG00000002549.2      
ENST00000002596.1 ENSG00000002587.1      
ENST00000265522.1 ENSG00000002746.2      
ENST00000265854.2 ENSG00000002822.2      
ENST00000004103.1 ENSG00000002933.1      
ENST00000262820.2 ENSG00000003096.3      
ENST00000322001.1 ENSG00000003096.3      
ENST00000265576.2 ENSG00000003147.4      
ENSMUST00000025370.2 ENSMUSG00000024493.2
ENSMUST00000025374.1 ENSMUSG00000024497.1
ENSMUST00000025375.2 ENSMUSG00000024498.2
ENSMUST00000052934.2 ENSMUSG00000024498.2
ENSMUST00000025377.2 ENSMUSG00000024500.2
ENSMUST00000069341.1 ENSMUSG00000024500.2
ENSMUST00000025379.1 ENSMUSG00000024501.2
ENSMUST00000070288.1 ENSMUSG00000024501.2
ENSMUST00000025381.1 ENSMUSG00000024503.1
Ver fichero

Péptidos y secuencias peptídicas

Utilizando una vez más los recursos que nos proporciona Ensembl, obtenemos dos ficheros, uno gff y otro fasta, para cada especie, a partir de las listas de identificadores de tránscritos de las que disponíamos para elaborar el fichero anterior 


HUMANO: estrucura exónica

Y       EnsEMBL exon    2448676 2448940 .       +       .       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00001334555.1";
Y       EnsEMBL exon    2467080 2467168 .       +       .       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000891589.3";
Y       EnsEMBL CDS     2467108 2467168 .       +       0       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000891589.3";
Y       EnsEMBL start_codon     2467108 2467110 .       +       .       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000891589.3";
Y       EnsEMBL exon    2474245 2474817 .       +       .       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000652223.1";
Y       EnsEMBL CDS     2474245 2474817 .       +       2       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000652223.1";
Y       EnsEMBL exon    2488266 2488415 .       +       .       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000652224.1";
Y       EnsEMBL CDS     2488266 2488415 .       +       2       gene_id "ENSG00000067646.1"; transcript_id "ENST00000155093.1"; exon_id "ENSE00000652224.1";
Ver fichero 

HUMANO: secuencia peptídica e identificadores

>ENST00000296205.3|ENSG00000163867.4|ENSP00000296205.3 assembly=NCBI34|chr=1|strand=reverse|peptide translation of coding sequence
QNMNFKYVGRYIKNIAYLFLKITVIQIFHSDLPMPNEKNDAELDSPPSKKKRLGFFQTYDTEYLKVGFIICPGSKESSPR
PQCVICGEILSSENMKPANLSHHLKTKHSELENKPVDFFEQKSLEMECQNSSLKKCLLVEKSLVKASYLIAFQTAASKKP
FSIAEELIKPYLVEMCSEVLGSSAGDKMKTIPLSNVTIQHRIDELSADIEDQLIQKVRESKWFALQIDESSEISNITLLL
CYIRFIDYDCRDVKEELLFCIEMPTQITGFEIFELINKYIDSKSLNWKHCVGLCTDGAASMTGRYSGLKAKIQEVAMNTA
AFTHCFIHRERLVAEKLSPCLHKILLQSAQILSFIKSNALNSRMLTILCEEMGSEHVSLPLHAEVRWISRGRMLKRLFEL
RHEIEIFLSQKHSDLAKYFHDEEWVGKLAYLSDIFSLINELNLSLQGTLTTFFNLCNKIDVFKRKLKMWLKRTQENDYDM
FPSFSEFSNSSGLNMTDITRIIFEHLEGLSQVFSDCFPPEQDLRSGNLWIIHPFMNHQNNNLTDFEEEKLTELSSDLGLQ
ALFKSVSVTQFWINAKTSYPELHERAMKFLLPFSTVYLCDAAFSALTESKQKNLLGSGPALRLAVTSLIPRIEKLVKEKE
*

>ENST00000311990.1|ENSG00000163867.4|ENSP00000311570.1 assembly=NCBI34|chr=1|strand=reverse|peptide translation of coding sequence
MNSSVGDLGVGGCSLWDDPARFIVVPAAYALALGLGLPANVAALAMFIRSGGRLGQALLLYLFNLALVDEFFTLTLQLWL
TYYLGLARRPPATRPGPPTTCPPMRRWSSPRSSACAAAASYAVPGPGRLPAWPGAYGAPRALPAPSPGWRAWPLPAWSTA
GQARGWPPPRWPSRPPSCWCSRPT*
Ver fichero 

RATÓN: estructura exónica

10_random_NT_078648     EnsEMBL exon    81887   82026   .       +       .       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309539.2";
10_random_NT_078648     EnsEMBL CDS     81887   82026   .       +       0       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309539.2";
10_random_NT_078648     EnsEMBL start_codon     81887   81889   .       +       .       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309539.2";
10_random_NT_078648     EnsEMBL exon    91738   91790   .       +       .       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309523.1";
10_random_NT_078648     EnsEMBL CDS     91738   91790   .       +       1       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309523.1";
10_random_NT_078648     EnsEMBL exon    94148   94241   .       +       .       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309503.1";
10_random_NT_078648     EnsEMBL CDS     94148   94241   .       +       2       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000309503.1";
10_random_NT_078648     EnsEMBL exon    95497   95630   .       +       .       gene_id "ENSMUSG00000035765.2"; transcript_id "ENSMUST00000039608.2"; exon_id "ENSMUSE00000401274.1";
Ver fichero 

RATÓN: secuencia peptídica e identificadores

>ENSMUST00000027195.1|ENSMUSG00000026034.1|ENSMUSP00000027195.1 assembly=NCBIM32|chr=1|strand=reverse|peptide translation of coding sequence
MRHSKRTYCPDWDERDWDYGTWRSSSSHKRKKRSHSSAREQKRCRYDHSKTTDSYYLESRSINEKAYHSRRYVDEYRNDY
MGYEPGHPYGEPGSRYQMHSSKSSGRSGRSSYKSKHRSRHHTSQHHSHGKSHRRKRSRSVEDDEEGHLICQSGDVLSARY
EIVDTLGEGAFGKVVECIDHKVGGRRVAVKIVKNVDRYCEAAQSEIQVLEHLNTTDPHSTFRCVQMLEWFEHRGHICIVF
ELLGLSTYDFIKENSFLPFRMDHIRKMAYQICKSVNFLHSNKLTHTDLKPENILFVKSDYTEAYNPKMKRDERTIVNPDI
KVVDFGSATYDDEHHSTLVSTRHYRAPEVILALGWSQPCDVWSIGCILIEYYLGFTVFPTHDSREHLAMMERILGPLPKH
MIQKTRKRRYFHHDRLDWDEHSSAGRYVSRRCKPLKEFMLSQDAEHELLFDLIGKMLEYDPAKRITLKEALKHPFFYPLK
KHT*

>ENSMUST00000040132.2|ENSMUSG00000038242.2|ENSMUSP00000046541.2 assembly=NCBIM32|chr=1|strand=forward|peptide translation of coding sequence
DPIQLLFYVNGQKVVEKNVDPEMMLLPYLRKNLRLTGTKYGCGGGGCGACTVMISRYNPSTKAIRHHPVNACLTPICSLH
GTAVTTVEGLGNTRTRLHPIQERIAKCHGTQCGFCTPGMVMSMYALLRNHPEPTLDQLTDALGGNLCRCTGYRPIIDACK
TFLVVTELSPENKRYYDDNFCHFLFLKTSPELFSEEEFLPLDPTQELIFPPELMRMAEESQNTVLTFRGERTTWIAPGTL
NDLLELKMKHPSAPLVIGNTYLGLHMKFTDVSYPIIISPARILELFVVTNTKQGLTLGAGLSLTQVKNVLSDVVSRLPKE
KTQIYCALLKQLKTLAGQQIRNVAVGGHIISRLPTSDLNPILGIGNCILNVASTDINGKAVQVTCLSFLGPNHSETWLFI
LKEWFLSPQREFVSAFRQAQCHQNALPDVNAGMRVLFREGTDVIEELSIAYGGVGPTTVSAQRSCQQLLGRRWNALMLDE
ACRLLLDEVSLPGSALGGKVEFRRTLIVSLFFKFYLEVLQELKADQKLPPESTVSALGDSDRCSHTRSWTHIDSHQPLQD
PVGRPIMHLSGLKHATGEAVFCDDIPRVDKELFMALVTSTRAHARIISIDSSEVLDLPGVVDVITAEDIPGNNGEEDDKL
LAVDKVLCVGQVICAVVAETDVQAKRATEKIKITYEDLKPVIFTIEVSWAAPCAPGVLLHLPLTGTVHVGGQEHFYMETQ
RVLVIPKTEDKELDMYVSTQDPAHVQKTVSSTLNIPISRITCHVKRVGGGFGGKVGRPAVFGAIAAVGAVKTGHPIRLVL
DREDDMLITGGRHPLFAKYKVSATSEVIKMPHFVTEFLVLKLENAYKIRNLRLRGRACMTNLPSNTAFRGFGFPQGALVT
ESCITAVAAKCGLPPEKIREKNMYKTVDKTIYKQAFNPDPLIRCWNECLDKSSFHIRRTRVDEFNKKSYWKKRGIAIVPM
KFSVGFAATSYHQAAALVHIYTDGSVLVAHGGNELGQGIHTKMLQVASRELKIPLSYLHICETSTTTVPNTIATAASVGA
DVNGRAVQNACQILLKRLEPVIKKNPEGTWRDWVEAAFEKRISLSATGYFRGYKAFMDWEKGEGDPFPYYVYGAACSEVE
IDCLTGAHKKIRTDIVMDACCSLNPAIDIGQIEGAFIQGMGLYTTEELLYSPEGVLYSRSPDKYKIPTVTDVPEQFNVSL
LPSSQTPLTLYSSKGLGESGMFLGSSVFFAIVDAVAAARRQRDIAEDFTVKSPATPEWVRMACADRFTDMV
Ver fichero

Elaboración del programa

Introducción de datos

El primer punto en la elaboración del programa pasa por introducir los datos iniciales y declarar las variables necesarias: 


print "enter the gene id or the transcript id: \n";
my $id = < STDIN >;
chomp ($id);

print "enter the organism: \n";
my $especie = < STDIN >;
chomp ($especie);

print "enter whether your sequence is a gene or a transcript: \n";
my $tog = < STDIN >;
chomp($tog);

En este caso, las variables que necesitaremos para iniciar la construcción del programa serán: el identificador del tránscrito o gen, la especie a la que pertenece, y el tipo de secuencia, es decir, si se trata de un tránscrito o de un gen.


Relación gen - tránscrito

El objetivo de este paso es obtener el identificador complementario al que hemos introducido inicialmente.

Para ello y, a partir del fichero "De tránscritos a genes ", realizaremos un hash. Este procedimiento nos permitirá obtener el tránscrito correspondiente al gen o viceversa, en función de la variable $tog (tránscrito o gen) que hayamos introducido.



if (!open(TRANSCRIPTSGENEID,"< transcripts_to_genes.gff")) {
    print "t_to_g.pl: no se puede abrir el fichero transcripts_to_genes.gff\n";
    exit;
}

my %transcripts;              # hash donde registraremos la matriz transcript-gene
my %genes;                    # hash donde registraremos la matriz gene-transcript
my @ids;                      # vector donde registraremos los ids (de transcript o gene)

while (< TRANSCRIPTSGENEID >){
    chomp;
    @ids = split " ", $_;

    if ($ids[0] =~ m/(ENST|ENSMUST)/ && $ids[1] =~ m/(ENSG|ENSMUSG)/ ) {
        push ( @{$transcripts{$ids[0]}}, $ids[1] );
        push ( @{$genes{$ids[1]}}, $ids[0] );
    }
}

my $gene_id;

print "tog= $tog\n";

if ( uc($tog) eq 'TRANSCRIPT' ) {
    my $string = join " ", @{$transcripts{$id}};
    print "---> $id es el transcript correspondiente al gen ",$string,"\n";
    $gene_id = $string;
}

if ( uc($tog) eq 'GENE' ){
    my $string = join " ", @{$genes{$id}};
    print "---> $id es el gen correspondiente al transcript ",$string,"\n";
   $gene_id = $id;
}

close TRANSCRIPTSGENEID;

Obtención del gen ortólogo

En este momento, ya tenemos un gen y sus tránscritos correspondientes, por lo que pasaremos a buscar el ortólogo del gen en la otra especie. Esto lo haremos utilizando otro hash, y basándonos en otro fichero, el archivo "Ortólogos"


if (!open(ORTOLOGOS,"< ortologos.fa")) {
  print "ortologos.pl: no se puede abrir el fichero ortologos.fa\n";
  exit;
}

my %Hortologos;              # hash donde registraremos la matriz de ortologos human-mouse
my %Mortologos;              # hash donde registraremos la matriz de ortologos mouse-human
my @genes;                   # vector donde registraremos los gene_id

while (< ORTOLOGOS >){
    chomp;

    @genes = split "\t", $_;

    if ($genes[0] =~ m/ENS/ && $genes[1] =~ m/ENSMUS/ ) {
        push ( @{$Hortologos{$genes[0]}}, $genes[1] );
        push ( @{$Mortologos{$genes[1]}}, $genes[0] );
    }
}

my $gene_id2;

if ( uc($especie) eq 'HUMAN' ) {
    $gene_id2 = join " ", @{$Hortologos{$gene_id}};
    print "---> $gene_id is ortholog to ",$gene_id2,"\n";
}

if ( uc($especie) eq 'MOUSE' ){
    $gene_id2 = join " ", @{$Mortologos{$gene_id}};
    print "---> $gene_id is ortholog to ",$gene_id2,"\n";
}

close(ORTOLOGOS);

Relación tránscrito - gen ortólogo

Una vez contamos con el gen ortólogo, el siguiente paso es obtener el tránscrito correspondiente. El proceso es otra vez un hash, que podemos aprovechar del primer paso, por lo que no tendremos que volver a realizarlo, sólo será necesario lo siguiente


my $string = join " ", @{$genes{$gene_id2}};

print "---> $gene_id2 es el gen correspondiente al transcript ",$string,"\n";

Pasos previos al alineamiento

Para que sea posible el alineamiento de secuencias peptídicas, necesitamos relacionar los identificadores de tránscritos con sus secuencias de aminoácidos correspondientes. Para ello usaremos dos funciones para leer los ficheros: get_gff y get_seq


        my $i = 0;
        my $j = 0;

        while ($i < scalar(@htrans)){
            while ($j < scalar(@mtrans)){
                my $hpid = $htrans[$i];
                my $mpid = $mtrans[$j];

                my $hgff = get_gff($hpid, "human_transcripts.gff");
                my $mgff = get_gff($mpid, "mouse_transcripts.gff");


                if ( scalar( @$hgff ) != 1 && scalar( @$mgff ) != 1 ){

                    my $hseq = get_seq($hpid,"/disc8/home/marbla04/human_peptides.fa");
                    my $mseq = get_seq($mpid,"/disc8/home/marbla04/mouse_peptides.fa");

                    open(FA,">pep.fa");
                    print FA ">$hpid\n";
                    print FA @$hseq;
                    print FA ">$mpid\n";
                    print FA @$mseq;
                    close(FA);



sub get_gff{
    my $t_id = $_[0];
    my $file = $_[1];
    my $gff;
    my @lines;

    open(GFF,"<$file");

    while(< GFF >){

        if ( /$t_id/ && /CDS/ ){
            my @entries = split;
            my $linea = $t_id."\tensembl\tCDS\t".
                $entries[3]."\t".$entries[4]."\t.\t".
                $entries[6]."\t".$entries[7]."\t".$t_id."\n";
            push(@lines, $linea);
        }
    }

    close(GFF);

    $gff = \@lines;
    return $gff;
}



sub get_seq{
    my $t_id = $_[0];
    my $file = $_[1];
    my $sequence;
    my @lines;

    open(PEP,"<$file");

  FILE:
    while(< PEP >){

        if ( /\>/ && /$t_id/ ){
        
            while(< PEP >){

                if ( !/\>/ ){
                    push( @lines, $_ );
                }
                else{
                    last FILE;
                }
            }
        }
    }
    $sequence = \@lines;

    close(PEP);

    return $sequence;
}

Alineamiento de las secuencias peptídicas

Para llevar a cabo este punto, lo único que tenemos que hacer es llamar a un programa de alineamiento de secuencias, en este caso clustalw en su versión 1.8, de la siguiente manera:


system("/disc8/bin/clustalw1.8/clustalw pep.fa");

open(HGFF,"> hum.gff");
print HGFF @$hgff;
close(HGFF);

open(MGFF,"> mus.gff");
print MGFF @$mgff;
close(MGFF);

Estudio de la conservación de tránscritos ortólogos

Para este último punto utilizaremos otro programa, exstral, de la siguiente manera:


                    open(EX, "/disc8/home/marbla04/exstral.pl  hum.gff mus.gff pep.aln | ");

                    while(< EX >){
                        print "-->".$_;
                        chomp;
                        next if (!/\#/);
                        my @result = split '\s+', $_;
                        if ($result[8] == 1){
                            print "OUTPUT TRAN $result[1] - $result[2] \t GENE $id - ",$gene_id2,"\n";
                        }
                    }
                }
                else{
                    if ( scalar( @$hgff ) == 1 ){
                        print "$hpid tiene un solo exon \n"}
                    if ( scalar( @$mgff ) == 1 ){
                        print "$mpid tiene un solo exon \n"}
                }
                $j = $j + 1;
            }
            $i = $i + 1;
        }
        $orto_index = $orto_index + 1;
    }  # fin del bucle ORTO_GENE



De este modo se podría aplicar el programa para estudiar la conservación de un tránscrito concreto. Nosotras hemos aplicado otra forma del programa que se basa en los mismos pasos, pero lo aplica a un fichero con un subconjunto de genes humanos para estudiar la conservación de los tránscritos en los ortólogos en ratón. En el link siguiente se muestra nuestro programa:

Ver programa