Skip to content

GitLab

Commit 50c50e02 authored by Javier Santoyo-Lopez's avatar Javier Santoyo-Lopez
Browse files

changed write_acedb so that transcripts 

can be written on the acedb complementary strand
parent 46a446ac
Hide whitespace changes
Inline Side-by-side
Showing with 48 additions and 34 deletions
modules/Bio/EnsEMBL/DB/ContigI.pm
View file @ 50c50e02
...@@ -523,21 +523,20 @@ sub get_AnnSeq { ...@@ -523,21 +523,20 @@ sub get_AnnSeq {
Title : write_acedb Title : write_acedb
Usage : $contig->write_acedb(\*FILEHANDLE); Usage : $contig->write_acedb(\*FILEHANDLE);
$contig->write_acedb(\*FILEHANDLE,$ace_seq_name, $type, $supp_evid); $contig->write_acedb(\*FILEHANDLE,$ace_seq_name, $type, $supp_evid, $revcom);
$ace_seq_name refers to the name of the aceDB-clone name Function: Dumps exons, transcript and gene objects of a contig in acedb format
$type refers to the type of gene in the ensEMBL database
$supp_evid, indicates supporting evidences
by default: $ace_seq_name is the accession number
$type is ensembl
$supp_evid is 0
Function: Dumps exon, transcript and gene objects in acedb format
Returns : Returns :
Args : Args : $ace_seq_name, refers to the name of the aceDB-clone name
$type refers, to the type of gene in the ensEMBL database
$supp_evid, indicates supporting evidences
$revcom, if true writes the sequence coordinates for the complementary strand
By default: $ace_seq_name is the accession number, $type is 'ensembl',
$supp_evid is 0.
=cut =cut
sub write_acedb { sub write_acedb {
my ($self,$fh,$seqname,$type,$supp_evid) = @_; my ($self, $fh, $seqname, $type, $supp_evid, $revcom) = @_;
my $contig_id = $self->id(); my $contig_id = $self->id();
...@@ -556,7 +555,7 @@ sub write_acedb { ...@@ -556,7 +555,7 @@ sub write_acedb {
GENE: GENE:
foreach my $gene ( @genes ){ foreach my $gene ( @genes ){
my $gene_id=$gene->id; my $gene_id = $gene->id;
# get all the transcripts of this gene. # get all the transcripts of this gene.
my @trans_in_gene = $gene->each_Transcript; my @trans_in_gene = $gene->each_Transcript;
...@@ -570,7 +569,7 @@ sub write_acedb { ...@@ -570,7 +569,7 @@ sub write_acedb {
# for each transcript # for each transcript
TRANSCRIPT: TRANSCRIPT:
foreach my $trans ( @trans_in_gene ) { foreach my $trans ( @trans_in_gene ) {
my $trans_id=$trans->id; my $trans_id = $trans->id;
# get all exons of this transcript # get all exons of this transcript
my @exons = $trans->each_Exon; my @exons = $trans->each_Exon;
...@@ -578,57 +577,72 @@ sub write_acedb { ...@@ -578,57 +577,72 @@ sub write_acedb {
# get transcript exons which belong to the contig # get transcript exons which belong to the contig
my @exons_in_contig; my @exons_in_contig;
EXON:
foreach my $exon ( @exons ) { foreach my $exon ( @exons ) {
if ( $exon->contig_id eq $contig_id ) { if ( $exon->contig_id eq $contig_id ) {
push ( @exons_in_contig,$exon ); push ( @exons_in_contig, $exon );
} else { next EXON; } }
} }
if (@exons_in_contig) { if (@exons_in_contig) {
# check the strand and get the coordinates my $tstart;
my $tend;
my $ace_tstart;
my $ace_tend;
my $tstrand = $exons_in_contig[0]->strand; my $tstrand = $exons_in_contig[0]->strand;
my ($tstart,$tend);
# check the strand and get the coordinates
if( $tstrand == 1 ) { if( $tstrand == 1 ) {
$tstart = $exons_in_contig[0]->start; $tstart = $exons_in_contig[0]->start;
$tend = $exons_in_contig[$#exons_in_contig]->end; $tend = $exons_in_contig[$#exons_in_contig]->end;
} else { } else {
$tstart = $exons_in_contig[0]->end; $tstart = $exons_in_contig[0]->end;
$tend = $exons_in_contig[$#exons_in_contig]->start; $tend = $exons_in_contig[$#exons_in_contig]->start;
}
unless ($revcom){
$ace_tstart = $tstart;
$ace_tend = $tend;
}else{
# remaping transcript coordinates in the complementary strand
my $contig_length = $self->length;
$ace_tstart = $contig_length - $tstart;
$ace_tend = $contig_length - $tend;
} }
# start .ace file printing... # start .ace file printing...
print $fh "Sequence $seqname\n";
print $fh "Subsequence $trans_id $tstart $tend\n\n";
# print coordinates of the transcript relative to the contig # print coordinates of the transcript relative to the contig
print $fh "Sequence $trans_id\nSource $seqname\nCDS\nCDS_predicted_by EnsEMBL\nMethod EnsEMBL\n"; print $fh "Sequence $seqname\n";
print $fh "Subsequence $trans_id $ace_tstart $ace_tend\n\n";
# print coordinates of each exon relative to the transcript # print coordinates of each exon relative to the transcript
print $fh "Sequence $trans_id\nSource $seqname\n";
print $fh "CDS\nCDS_predicted_by EnsEMBL\nMethod EnsEMBL\n";
foreach my $exon ( @exons_in_contig ) { foreach my $exon ( @exons_in_contig ) {
if( $tstrand == 1 ) { if( $tstrand == 1 ) {
print $fh "Source_Exons ", ($exon->start - $tstart + 1)," ",($exon->end - $tstart +1), "\n"; print $fh "Source_Exons ", ($exon->start - $tstart + 1),
" ",($exon->end - $tstart +1), "\n";
} else { } else {
print $fh "Source_Exons ", ($tstart - $exon->end +1 ), " ",($tstart - $exon->start+1),"\n"; print $fh "Source_Exons ", ($tstart - $exon->end +1 ),
" ",($tstart - $exon->start+1), "\n";
} }
} }
# indicate end or start not found for transcript across several contigs # indicate end or start not found for transcript across several contigs
if ($exons[0] != $exons_in_contig[0]) { if ($exons[0] != $exons_in_contig[0]) {
print $fh "Start_not_found\n"; print $fh "Start_not_found\n";
} }
if ($exons[$#exons] != $exons_in_contig[$#exons_in_contig]) { if ($exons[$#exons] != $exons_in_contig[$#exons_in_contig]) {
print $fh "End_not_found\n"; print $fh "End_not_found\n";
} }
# URL object tag
print $fh "Web_location Sanger-EnsEMBL\n";
print $fh "\n\n";
print $fh "\n\n";
} else { } else {
print STDERR "'$trans_id' has no exons in '$seqname'\n"; print STDERR "'$trans_id' has no exons in '$seqname'\n";
next TRANSCRIPT; next TRANSCRIPT;
} }
}
}
} }
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment