File Coverage

blib/lib/Microarray/GEO/SOFT/GDS.pm

Criterion	Covered	Total	%
statement	12	12	100.0
branch			n/a
condition			n/a
subroutine	4	4	100.0
pod			n/a
total	16	16	100.0

package Microarray::GEO::SOFT::GDS;

use List::Vectorize qw(!table);
use Carp;
use strict;

use base "Microarray::GEO::SOFT";

1;


sub new {

	my $invocant = shift;
	my $class = ref($invocant) || $invocant;
	my $self = { "file" => "",
	             "verbose" => 1,
	             "use_identifier" => 0,
	             @_ };
	bless($self, $class);

	return $self;

}

sub parse {

	my $self = shift;

	my $fh;
	if(! List::Vectorize::is_glob_ref($self->{file})) {

		open F, $self->{file} or croak "cannot open $self->{file}.\n";
		$fh = \*F;
	}
	else {
		$fh = $self->{file};
	}

	$self->_parse_dataset($fh);

	return 1;
}

sub _parse_dataset {

	my $self = shift;

	my $fh = shift;

	Microarray::GEO::SOFT::_set_fh($self->{verbose});

	my $accession;
	my $title;
	my $platform;
	my $field_explain;
	my $table_colnames = [];
	my $table_rownames = [];
	my $table_colnames_explain_hash = {};
	my $table_matrix = [];

	while(my $line = <$fh>) {

		chomp $line;
		if($line =~/^\^DATASET = (GDS\d+)$/) {
			$accession = $1;
		}

		if($line =~/^!dataset_title = (.*?)$/) {
			$title = $1;
		}

		if($line =~/^!dataset_platform = (GPL\d+)$/) {
			$platform = $1;
		}

		if($line =~/^#(GSM\d+) = (Value for GSM\d+: )?(.*?)$/) {
			$table_colnames_explain_hash->{$1} = $3;
		}

		if($line =~/^!dataset_table_begin$/) {

			$line = <$fh>;
			chomp $line;

			@$table_colnames = split "\t", $line, -1;
			shift(@$table_colnames);
			shift(@$table_colnames);

			while($line = <$fh>) {

				if($line =~/^!dataset_table_end$/) {
					last;
				}

				chomp $line;
				my @tmp = split "\t", $line, -1;

				my $uid = shift(@tmp);

				# the second column in the matrix is identifier
				my $identifier = shift(@tmp);

				# do not recommond to use identifier
				# it is better to convert IDs using id_convert subroutine
				if($self->{use_identifier}) {
					push(@$table_rownames, $identifier);
				}
				else {
					push(@$table_rownames, $uid);
				}
				push(@$table_matrix, [@tmp]);

			}


		}
		if($line =~/^!dataset_table_end$/) {
			last;
		}

	}

	my $n_row = len($table_rownames);
	my $n_col = len($table_colnames);

	my $table_colnames_explain = [];
	for (@$table_colnames) {
		push(@$table_colnames_explain, $table_colnames_explain_hash->{$_});
	}

	print "Dataset info:\n";
	print "  Accession: $accession\n";
	print "  Title: $title\n";
	print "  Rows: $n_row\n";
	print "  Columns: $n_col\n";
	print "\n";

	#my $table_rownames_sorted = sort_array($table_rownames, sub {$_[0] cmp $_[1]});
	#my $table_rownames_sorted_index = order($table_rownames, sub {$_[0] cmp $_[1]});
	#my $table_matrix_sorted = subset($table_matrix, $table_rownames_sorted_index);

	$self->set_meta( accession => $accession,
	                 title     => $title,
					 platform  => $platform );
	$self->set_table( rownames => $table_rownames,
	                  colnames => $table_colnames,
					  colnames_explain => $table_colnames_explain,
					  matrix   => $table_matrix );

	Microarray::GEO::SOFT::_set_to_std_fh();

	return $self;

}

sub get_subset {

	my $self = shift;
	my $arg = {"byrow" => rep(1, len($self->rownames)),
	           "bycol" => rep(1, len($self->colnames)),
			   @_};

	$arg->{byrow} = sapply($arg->{byrow}, sub{($_[0] != 0)+0});
	$arg->{bycol} = sapply($arg->{bycol}, sub{($_[0] != 0)+0});

	if(len($arg->{byrow}) != len($self->rownames)
	   or len($arg->{bycol}) != len($self->colnames)) {

		croak "ERROR: Do not fit the dimension of the matrix";
	}

	if(sum($arg->{bycol}) == len($self->colnames)) {
		$self->set_table( rownames => subset($self->rownames, which($arg->{byrow})),
		                  matrix => subset($self->matrix, which($arg->{byrow})) );
	}
	else {
		my $new_matrix = sapply($self->matrix, sub { subset($_[0], which($arg->{bycol})) });
		$new_matrix = subset($new_matrix, which($arg->{byrow}));
		$self->set_table( rownames => subset($self->rownames, which($arg->{byrow})),
		                  colnames => subset($self->colnames, which($arg->{bycol})),
		                  colnames_explain => subset($self->colnames_explain, which($arg->{bycol})),
		                  matrix => $new_matrix );

	}

	return $self->soft2exprset;

}

sub id_convert {

	my $self = shift;
	my $gpl = shift;
	my $to_id = shift;

	if($self->{use_identifier}) {
		croak "ERROR: You are not suitable to use 'id_convert' since you have set use_identifier to TURE. You can only use 'soft2expr' for further analysis.";
	}

	my $platform_id = $gpl->accession;

	if(ref($to_id) eq "Regexp") {
		my @match = grep {/$to_id/} @{$gpl->colnames};
		if(!scalar(@match)) {
			croak "ERROR: Cannot find ID ($to_id) in $platform_id\n";
		}
		elsif(scalar(@match) > 1) {
			carp "WARNING: Find more than one matched ID types (".(join ", ", @match).") under $to_id, only take the first one ($match[0])";
			$to_id = $match[0];
		}
		else {
			$to_id = $match[0];
		}
	}
	else {
		if(! is_element($to_id, $gpl->colnames)) {
			croak "ERROR: Cannot find ID ($to_id) in $platform_id\n";
		}
	}

	my $new_rownames = $gpl->_mapping($to_id, $self->rownames);
	my $eset = Microarray::ExprSet->new;
	$eset->set_feature($new_rownames);
	$eset->set_phenotype($self->colnames_explain);
	$eset->set_matrix($self->matrix);

	return $eset;

}

# change SOFT to ExprSet
# in fact, it is only for GDS
sub soft2exprset {

	my $self = shift;

	my $eset = Microarray::ExprSet->new;
	$eset->set_feature($self->rownames);
	$eset->set_phenotype($self->colnames_explain);
	$eset->set_matrix($self->matrix);

	return $eset;

}



__END__