Bio::Tools::Alignment::Consed provides methods and objects to dealwith the output from the Consed software suite. Specifically,takes an .ace file and provides objects for the results.A word about doublets: This module was written to accommodate a largeEST sequencing operation. In this case, EST's were sequenced from the3' and from the 5' end of the EST. The objective was to find aconsensus sequence for these two reads. Thus, a contig of two is whatwe wanted, and this contig should consist of the forward and reversereads of a getn clone. For example, for a forward designator of "F"and a reverse designator of "R", if the two reads chad1F and chad1Rwere in a single contig (for example Contig 5) it will be determinedthat the consensus sequence for Contig 5 will be the sequence forclone chad1.Doublets are good!This module parses .ace and related files. A detailed list of methodscan be found at the end of this document.I wrote a detailed rationale for design that may explain the reasonswhy some things were done the way they were done. That document isbeyond the scope of this pod and can probably be found in thedirectory from which this module was 'made' or athttp://www.dieselwurks.com/bioinformatics/consedpm_documentation.pdf.Note that the POD in that document might be old but the originalrationale still stands.

Title : new(-acefile => $path_to_some_acefile, -verbose => "1") Usage : $o_consed = Bio::Tools::Alignment::Consed-> new(-acefile => $path_to_some_acefile, -verbose => "1"); Function: Construct the Bio::Tools::Alignment::Consed object. Sets verbosity for the following procedures, if necessary: 1. Construct a new Bio::Tools::Alignment::Trim object, to handle quality trimming 2. Read in the acefile and parse it
Returns : A reference to a Bio::Tools::Alignment::Consed object.
Args : A hash. (-acefile) is the filename of an acefile. If a full path
is not specified "./" is prepended to the filename and used from
instantiation until destruction. If you want
Bio::Tools::Alignment::Consed to be noisy during parsing of
the acefile, specify some value for (-verbose).

This method for setting verbosity has largely been superseded by asub-by-sub way, where for every sub you can provide a (-verbose)switch. I am doing converting this bit-by-bit so do not be surprisedif some subs do not honour this.

Title : get_filename() Usage : $o_consed->get_filename(); Function: Returns the name of the acefile being used by theBio::Tools::Alignment::Consed object. Returns : A scalar containing the name of a file. Args : None.

Title : count_sequences_with_grep() Usage : $o_consed->count_sequences_with_grep(); Function: Use /bin/grep to scan through the files in the ace project dir and count sequences in those files. I used this method in the development of this module to verify that I was getting all of the sequences. It works, but it is (I think) unix-like platform dependent. Returns : A scalar containing the number of sequences in the ace project directory. Args : None.

If you are on a non-UNIX platform, you really do not have to usethis. It is more of a debugging routine designed to address veryspecific problems.This method was reimplemented to be platform independent with a pureperl implementation. The above note can be ignored.

Title : get_quality_array($contig_keyname) Usage : $o_consed->get_quality_array($contig_keyname); Function: Returns the quality for the consensus sequence for the given contig as an array. See get_quality_scalar to get this as a scalar. Returns : An array containing the quality for the consensus sequence with the given keyname. Args : The keyname of a contig. Note: This is a keyname. The key would normally come from get_contigs.

Returns an array, not a reference. Is this a bug? thinking No.Well, maybe. Why was this developed like this? I was using FreezeThawfor object persistence, and when it froze out these arrays it took along time to thaw it. Much better as a scalar.See get_quality_scalar()

Title : get_quality_scalar($contig_keyname) Usage : $o_consed->get_quality_scalar($contig_keyname); Function: Returns the quality for the consensus sequence for the given contig as a scalar. See get_quality_array to get this as an array. Returns : An scalar containing the quality for the consensus sequence with the given keyname. Args : The keyname of a contig. Note this is a _keyname_. The key would normally come from get_contigs.

Why was this developed like this? I was using FreezeThaw for objectpersistence, and when it froze out these arrays it took a coon's ageto thaw it. Much better as a scalar.See get_quality_array()

Title : freeze_hash() Usage : $o_consed->freeze_hash();
Function: Use Ilya's FreezeThaw module to create a persistent data
object for this Bio::Tools::Alignment::Consed data
structure. In the case of AAFC, we use
Bio::Tools::Alignment::Consed to pre-process bunches of
sequences, freeze the structures, and send in a harvesting
robot later to do database stuff.
Returns : 0 or 1;
Args : None.

Title : get_members($contig_keyname) Usage : $o_consed->get_members($contig_keyname); Function: Return the _names_ of the reads in this contig. Returns : An array containing the names of the reads in this contig. Args : The keyname of a contig. Note this is a keyname. The keyname would normally come from get_contigs.

Title : get_members_by_name($some_arbitrary_name) Usage : $o_consed->get_members_by_name($some_arbitrary_name); Function: Return the names of the reads in a contig. This is the name given to $contig{key} based on what is in the contig. This is different from the keys retrieved through get_contigs(). Returns : An array containing the names of the reads in the contig with this name. Args : The name of a contig. Not a key, but a name.

Highly inefficient. use some other method if possible.See get_contigs()

Title : get_contig_number_by_name($some_arbitrary_name) Usage : $o_consed->get_contig_number_by_name($some_arbitrary_name); Function: Return the names of the reads in a contig. This is the name given to $contig{key} based on what is in the contig. This is different from the keys retrieved through get_contigs(). Returns : An array containing the names of the reads in the contig with this name. Args : The name of a contig. Not a key, but a name.

Title : get_sequence($contig_keyname) Usage : $o_consed->get_sequence($contig_keyname); Function: Returns the consensus sequence for a given contig. Returns : A scalar containing a sequence. Args : The keyname of a contig. Note this is a key. The key would normally come from get_contigs.

Title : set_final_sequence($name,$some_sequence) Usage : $o_consed->set_final_sequence($name,$some_sequence); Function: Provides a manual way to set the sequence for a given key in the contig hash. Rarely used. Returns : 0 or 1; Args : The name (not the keyname) of a contig and an arbitrary string.

A method with a questionable and somewhat mysterious origin. May raisethe dead or something like that.

Title : set_reverse_designator($some_string) Usage : $o_consed->set_reverse_designator($some_string); Function: Set the designator for the reverse read of contigs in thisBio::Tools::Alignment::Consed object. Used to determine if contigs containing two reads can be named. Returns : The value of $o_consed->{reverse_designator} so you can check to see that it was set properly. Args : An arbitrary string.

Title : set_forward_designator($some_string) Usage : $o_consed->set_forward_designator($some_string); Function: Set the designator for the forward read of contigs in thisBio::Tools::Alignment::Consed object. Used to determine if contigs containing two reads can be named. Returns : The value of $o_consed->{forward_designator} so you can check to see that it was set properly. Args : An arbitrary string.

Title : set_trim_points_singlets_and_singletons() Usage : $o_consed->set_trim_points_singlets_and_singletons(); Function: Set the trim points for singlets and singletons based on quality. Uses the Bio::Tools::Alignment::Trim object. Use at your own risk because the Bio::Tools::Alignment::Trim object was designed specifically for me and is mysterious in its ways. Every time somebody other then me uses it a swarm of locusts decends on a small Central American village so do not say you weren't warned. Returns : Nothing. Args : None.

Title : set_trim_points_doublets() Usage : $o_consed->set_trim_points_doublets(); Function: Set the trim points for doublets based on quality. Uses the Bio::Tools::Alignment::Trim object. Use at your own risk because the Bio::Tools::Alignment::Trim object was designed specifically for me and is mysterious in its ways. Every time somebody other then me uses it you risk a biblical plague being loosed on your city. Returns : Nothing. Args : None. Notes : Working on exceptions here.

Title : set_doublets() Usage : $o_consed->set_doublets(); Function: Find pairs that have similar names and mark them as doublets and set the {name}. Returns : 0 or 1. Args : None.

A complicated subroutine that iterates over theBio::Tools::Alignment::Consed looking for contigs of 2. If the forwardand reverse designator are removed from each of the reads in{'member_array'} and the remaining reads are the same, {name} is setto that name and the contig's class is set as "doublet". If any ofthose cases fail the contig is marked as a "pair".

Title : set_singlets Usage : $o_consed->set_singlets(); Function: Read in a singlets file and place them into theBio::Tools::Alignment::Consed object. Returns : Nothing. Args : A scalar to turn on verbose parsing of the singlets file. Notes :

Title : set_quality_by_name($name,$quality) Usage : $o_consed->set_quality_by_name($name,$quality); Function: Deprecated. Make the contig with {name} have {'quality'} $quality. Probably used for testing. Returns : Nothing. Args : The name of a contig and a scalar for its quality. Notes : Deprecated.

Title : set_singlet_quality() Usage : $o_consed->set_singlet_quality(); Function: For each singlet, go to the appropriate file in phd_dir and read in the phred quality for that read and place it into {'quality'} Returns : 0 or 1. Args : None. Notes : This is the next subroutine that will receive substantial revision in the next little while. It really should eval the creation of Bio::Tools::Alignment::Phred objects and go from there.

Title : get_all_members() Usage : @all_members = $o_consed->get_all_members(); Function: Return a list of all of the read names in the Bio::Tools::Alignment::Consed object. Returns : An array containing all of the elements in all of the {'member_array'}s. Args : None. Notes :

Title : sum_lets($total_only) Usage : $statistics = $o_consed->sum_lets($total_only); Function: Provide numbers for how many sequences were accounted for in theBio::Tools::Alignment::Consed object. Returns : If a scalar is present, returns the total number of sequences accounted for in all classes. If no scalar passed then returns a string that looks like this: Singt/singn/doub/pair/mult/total : 2,0,1(2),0(0),0(0),4 This example means the following: There were 1 singlets. There were 0 singletons. There were 1 doublets for a total of 2 sequences in this class. There were 0 pairs for a total of 0 sequences in this class. There were 0 multiplets for a total of 0 sequences in this class. There were a total of 4 sequences accounted for in theBio::Tools::Alignment::Consed object. Args : A scalar is optional to change the way the numbers are returned. Notes:

Title : write_stats() Usage : $o_consed->write_stats(); Function: Write a file called "statistics" containing numbers similar to those provided in sum_lets(). Returns : Nothing. Write a file in $o_consed->{path} containing something like this:
0,0,50(100),0(0),0(0),100
Where the numbers provided are in the format described in the
documentation for sum_lets().
Args : None.
Notes : This might break platform independence, I do not know.

Title : _get_contig_name(\@array_containing_reads) Usage : $o_consed->_get_contig_name(\@array_containing_reads); Function: The logic for the set_doublets subroutine. Returns : The name for this contig. Args : A reference to an array containing read names. Notes : Depends on reverse_designator. Be sure this is set the way you intend.

Title : dump_hash() Usage : $o_consed->dump_hash(); Function: Use dumpvar.pl to dump out the Bio::Tools::Alignment::Consed object to STDOUT. Returns : Nothing. Args : None. Notes : I used this a lot in debugging.

Title : get_phreds() Usage : @phreds = $o_consed->get_phreds(); Function: For each doublet in the Bio::Tools::Alignment::Consed hash, go and get the phreds for the top and bottom reads. Place them into {top_phreds} and {bottom_phreds}. Returns : Nothing. Args : Nothing.

Requires parse_phd() and reverse_and_complement(). I realize that itwould be much more elegant to pull qualities as required but therewere certain "features" in the acefile that required a bit moredetailed work be done to get the qualities for certain parts of theconsensus sequence. In order to make _sure_ that this was doneproperly I wrote things to do all steps and then I used dump_hash()and checked each one to ensure expected behavior. I have never changedthis, so there you are.

Title : parse_phd($read_name) Usage : $o_consed->parse_phd($read_name); Function: Suck in the contents of a .phd file. Returns : A reference to an array containing the quality values for the read. Args : The name of a read. Notes : This is a significantly weak subroutine because it was always intended that these functions, along with the functions provided by get_phreds() be put into the Bio::SeqIO:phd module. This is done now but the Bio::Tools::Alignment::Consed module has not be rewritten to reflect this change.

my($self) = @_;
my($line,$in_contig,$in_quality,$contig_number,$top);
# make it easier to type $fhlwhile(defined($line=$self->_readline())){chomp$line;
# check if there is anything on this line# if not, you can stop gathering consensus sequenceif(!$line){# if the line is blank you are no longer to gather consensus # sequence or quality values$in_contig = 0;
$in_quality = 0;
}# you are currently gathering consensus sequenceelsif($in_contig){if($in_contig == 1){$self->debug("Adding $line to consensus of contig number $contig_number.\n");
$self->{'contigs'}->{$contig_number}->{'consensus'} .= $line;
}}elsif($in_quality){if(!$line){$in_quality = undef;
}else{# I wrote this in here because acefiles produced by# cap3 do not have a leading space like the acefiles# produced by phrap and there is the potential to have# concatenated quality values like this: 2020 rather# then 20 20 whre lines collide. Thanks Andrew for# noticing.if($self->{'contigs'}->{$contig_number}->{'quality'} &&
!($self->{'contigs'}->{$contig_number}->{'quality'} =~ m/\ $/)) { $self->{'contigs'}->{$contig_number}->{'quality'} .= " ";}$self->{'contigs'}->{$contig_number}->{'quality'} .= $line;
}}elsif($line =~ /^BQ/){$in_quality = 1;
}# the line /^CO/ like this:# CO Contig1 796 1 1 U# can be broken down as follows:# CO - Contig!# Contig1 - the name of this contig# 796 - Number of bases in this contig# 1 - Number of reads in this contig# 1 - number of base segments in this contig# U - Uncomplementedelsif($line =~ /^CO/){$line =~ m/^CO\ Contig(\d+)\ \d+\ \d+\ \d+\ (\w)/;$contig_number = $1;
if($2eq"C"){$self->debug("Contig $contig_number is complemented!\n");
}$self->{'contigs'}->{$contig_number}->{'member_array'} = [];
$self->{'contigs'}->{$contig_number}->{'contig_direction'} = "$2";
$in_contig = 1;
}# 000713# this BS is deprecated, I think.# haha, I am really witty. <ew>elsif($line =~ /^BSDEPRECATED/){$line =~ m/^BS\s+\d+\s+\d+\s+(.+)/;my$member = $1;
$self->{'contigs'}->{$contig_number}->{$member}++;
}# the members of the contigs are determined by the AF line in the ace fileelsif($line =~ /^AF/){$self->debug("I see an AF line here.\n");
$line =~ /^AF\ (\S+)\ (\w)\ (\S+)/;
# push the name of the current read onto the member array for this contigpush @{$self->{'contigs'}->{$contig_number}->{'member_array'}},$1;
# the first read in the contig will be named the "top" readif(!$top){$self->debug("\$top is not set.\n");
if($self->{'contigs'}->{$contig_number}->{'contig_direction'}eq"C"){$self->debug("Reversing the order of the reads. The bottom will be $1\n");
# if the contig sequence is marked as the# complement, the top becomes the bottom and$$self->{'contigs'}->{$contig_number}->{'bottom_name'} = $1;
$self->{'contigs'}->{$contig_number}->{'bottom_complement'} = $2;
$self->{'contigs'}->{$contig_number}->{'bottom_start'} = $3;
}else{$self->debug("NOT reversing the order of the reads. ".
"The top_name will be $1\n");
# if the contig sequence is marked as the# complement, the top becomes the bottom and$$self->{'contigs'}->{$contig_number}->{'top_name'} = $1;
$self->{'contigs'}->{$contig_number}->{'top_complement'} = $2;
$self->{'contigs'}->{$contig_number}->{'top_start'} = $3;
}$top = 1;
}else{# if the contig sequence is marked as the complement,# the top becomes the bottom and the bottom becomes# the topif($self->{'contigs'}->{$contig_number}->{'contig_direction'}eq"C"){$self->debug("Reversing the order of the reads. The top will be $1\n");
$self->{'contigs'}->{$contig_number}->{'top_name'} = $1;
$self->{'contigs'}->{$contig_number}->{'top_complement'} = $2;
$self->{'contigs'}->{$contig_number}->{'top_start'} = $3;
}else{$self->debug("NOT reversing the order of the reads. The bottom will be $1\n");
$self->{'contigs'}->{$contig_number}->{'bottom_name'} = $1;
$self->{'contigs'}->{$contig_number}->{'bottom_complement'} = $2;
$self->{'contigs'}->{$contig_number}->{'bottom_start'} = $3;
}$top = undef;
}}}return 0;

my($self) = @_;
# set the designators in the Bio::Tools::Alignment::Trim object$self->{'o_trim'}->set_designators($self->{'reverse_designator'},
$self->{'forward_designator'});
foreachmy$key_contig(sortkeys %{$self->{'contigs'}}){# if there is a member array (why would there not be? This should be a die()able offence# but for now I will leave itif($self->{'contigs'}->{$key_contig}->{'member_array'}){# if there are two reads in this contig # i am pretty sure that this is wrong but i am keeping it for reference# if (@{$self->{'contigs'}->{$key_contig}->{'member_array'}} == 2 || !$self->{'contigs'}->{$key_contig}->{'class'}) {# <seconds later># <nod> WRONG. Was I on crack?if(@{$self->{'contigs'}->{$key_contig}->{'member_array'}} == 2){$self->{'contigs'}->{$key_contig}->{'num_members'} = 2;
$self->debug("\tThere are 2 members! Looking for the contig name...\n");
my$name = _get_contig_name($self,$self->{'contigs'}->{$key_contig}->{'member_array'});
$self->debug("The name is $name\n")ifdefined$name;
if($name){$self->{'contigs'}->{$key_contig}->{'name'} = $name;
$self->{'contigs'}->{$key_contig}->{'class'} = "doublet";
}else{$self->debug("$key_contig is a pair.\n");
$self->{'contigs'}->{$key_contig}->{'class'} = "pair";
}}# this is all fair and good but what about singlets?# they have one reads in the member_array but certainly are not singletonselsif(@{$self->{'contigs'}->{$key_contig}->{'member_array'}} == 1){# set the name to be the name of the read$self->{'contigs'}->{$key_contig}->{name} = @{$self->{'contigs'}->{$key_contig}->{'member_array'}}[0];
# set the number of members to be one$self->{'contigs'}->{$key_contig}->{num_members} = 1;
# if this was a singlet, it would already belong to the class "singlet"# so leave it alone# if it is not a singlet, it is a singleton! lablel it appropriatelyunless($self->{'contigs'}->{$key_contig}->{'class'}){$self->{'contigs'}->{$key_contig}->{'class'} = "singleton";
}}# set the multiplet characteristicselsif(@{$self->{'contigs'}->{$key_contig}->{'member_array'}} >= 3){$self->{'contigs'}->{$key_contig}->{'num_members'} = @{$self->{'contigs'}->{$key_contig}->{'member_array'}};
$self->{'contigs'}->{$key_contig}->{'class'} = "multiplet";
}$self->{'contigs'}->{$key_contig}->{'num_members'} = @{$self->{'contigs'}->{$key_contig}->{'member_array'}};
}}$self->{'doublets_set'} = "done";
return 0;
}# end set_doublets

# returns an array of singlet names# singlets have "singlet"=1 in the hashmy$self = shift;
if(!$self->{singlets_set}){$self->debug("You need to set the singlets before you get them. Doing that now.");
$self->set_singlets();
}my(@singlets,@array);
foreachmy$key(sortkeys %{$self->{'contigs'}}){# @array = @{$Consed::contigs{$key}->{'member_array'}};# somethimes a user will try to get a list of singlets before the classes for the rest of the# contigs has been set (see t/test.t for how I figured this out. <bah> # so either way, just return class=singletsif(!$self->{'contigs'}->{$key}->{'class'}){# print("$key has no class. why?\n");}elsif($self->{'contigs'}->{$key}->{'class'}eq"singlet"){push@singlets,$key;
}}return@singlets;

# returns an array of pair contig names# a pair is a contig of two where the names do not matchmy$self = shift;
my(@pairs,@array);
foreachmy$key(sortkeys %{$self->{'contigs'}}){if($self->{'contigs'}->{$key}->{'member_array'}){if(@{$self->{'contigs'}->{$key}->{'member_array'}} == 2 &&
$self->{'contigs'}->{$key}->{'class'}eq"pair"){push@pairs,$key;
}}}return@pairs;

my$self = shift;
if(!$self->{doublets_set}){$self->warn("You need to set the doublets before you can get them. Doing that now.");
$self->set_doublets();
}my@doublets;
foreach(sortkeys %{$self->{'contigs'}}){if($self->{'contigs'}->{$_}->{name} && $self->{'contigs'}->{$_}->{'class'}eq"doublet"){push@doublets,$_;
}}return@doublets;
}# end get_doublets

my$self = shift;
my$dumper = new Dumpvalue;
$self->debug("Bio::Tools::Alignment::Consed::dump_hash - ".
"The following is the contents of the contig hash...\n");
$dumper->dumpValue($self->{'contigs'});

User feedback is an integral part of the evolution of this and otherBioperl modules. Send your comments and suggestions preferably to oneof the Bioperl mailing lists. Your participation is much appreciated.

Please direct usage questions or support issues to the mailing list:bioperl-l@bioperl.orgrather than to the module maintainer directly. Many experienced and reponsive experts will be able look at the problem and quickly address it. Please include a thorough description of the problem with code and data examples if at all possible.

Title : reverse_recurse(\@source,\@destination) Usage : $o_consed->reverse_recurse(\@source,\@destination); Function: A recursive routine to reverse and complement an array of phred data. Returns : A reference to an array containing reversed phred data. Args : A reference to a source array and a reverence to a destination array.

Title : show_missing_sequence(); Usage : $o_consed->show_missing_sequence(); Function: Used by set_trim_points_doublets() to fill in quality values where consed (phrap?) set them to 0 at the beginning and/or end of the consensus sequences. Returns : Nothing. Args : None.

Acts on doublets only. Really very somewhat quite ugly. A disgustingkludge. insert pride here It was written stepwise with no real planbecause it was not really evident why consed (phrap?) was doing this.