frag results from the CPAN

frag

App-Sandy

view release on metacpan or search on metacpan

lib/App/Sandy/Simulator.pm view on Meta::CPAN

				$std_id => $id
			);

			# It is necessary to catch gene -> transcript relation
			# # TODO: Make a hash tarit for indexed fasta
			if (defined $fields[1]) {
				my $pid = $fields[1];

				# Chomp and remove id version
				$self->_clean_fasta_id(\$pid);

				$fasta_rtree{$id} = $pid;
			}
		} else {
			die "Error reading fasta file '$fasta': Not defined id"
				unless defined $id;
			$indexed_fasta{$id}{seq} .= uc($_);
		}
	}

	for (keys %indexed_fasta) {
		$indexed_fasta{$_}{size} = length $indexed_fasta{$_}{seq};
	}

	unless (%indexed_fasta) {
		die "Error parsing '$fasta'. Maybe the file is empty\n";
	}

	$fh->close
		or die "Cannot close file $fasta: $!\n";

	$self->_set_fasta_rtree(%fasta_rtree) if %fasta_rtree;
	return \%indexed_fasta;
}

sub _build_fasta {
	my $self = shift;
	my $fasta = $self->fasta_file;

	log_msg ":: Indexing fasta file '$fasta' ...";
	my $indexed_fasta = $self->_index_fasta;

	# Validate genome about the read size required
	log_msg ":: Validating fasta file '$fasta' ...";
	# Entries to remove
	my @blacklist;

	unless ($self->truncate) {
		for my $id (keys %$indexed_fasta) {
			my $index_size = $indexed_fasta->{$id}{size};
			my $class = ref $self->seq;

			if ($class eq 'App::Sandy::Seq::SingleEnd') {
				my $read_mean = $self->seq->read_mean;
				if ($index_size < $read_mean) {
					log_msg ":: Parsing fasta file '$fasta': Seqid sequence length (>$id => $index_size) lesser than required read mean ($read_mean)";
					delete $indexed_fasta->{$id};
					push @blacklist => $id;
				}
			} elsif ($class eq 'App::Sandy::Seq::PairedEnd') {
				my $fragment_mean = $self->seq->fragment_mean;
				if ($index_size < $fragment_mean) {
					log_msg ":: Parsing fasta file '$fasta': Seqid sequence length (>$id => $index_size) lesser than required fragment mean ($fragment_mean)";
					delete $indexed_fasta->{$id};
					push @blacklist => $id;
				}
			} else {
				croak "Unknown option '$class' for sequencing type\n";
			}
		}
	}

	unless (%$indexed_fasta) {
		die sprintf "Fasta file '%s' has no valid entry\n" => $self->fasta_file;
	}

	# If fasta_rtree has entries
	unless ($self->_has_no_fasta_rtree) {
		# Remove no valid entries from id -> pid relation
		$self->_delete_fasta_rtree(@blacklist) if @blacklist;
	}

	return $indexed_fasta;
}

sub _populate_fasta_tree {
	my $self = shift;

	# If fasta_rtree has entries
	unless ($self->_has_no_fasta_rtree) {
		# Build parent -> child ids relation
		my %fasta_tree;

		# Reverse fasta_rtree to pid -> \@ids
		for my $pair ($self->_fasta_rtree_pairs) {
			my ($id, $pid) = (@$pair);
			push @{ $fasta_tree{$pid} } => $id;
		}

		$self->_set_fasta_tree(%fasta_tree);
	}
}

sub _build_fasta_blacklist {
	my $self = shift;

	log_msg ":: Index unidentified (NNN..) regions";

	my $indexed_fasta = $self->_fasta;
	my $class = ref $self->seq;

	my %fasta_blacklist;
	my @blacklist_id;

	for my $id (keys %$indexed_fasta) {
		my $seq = $indexed_fasta->{$id}{seq};
		my $seq_len = $indexed_fasta->{$id}{size};

		my $offset = 0;
		my ($st, $en) = (0, 0);
		my $pos = 0;
		my $init = 0;
		my $len = 0;
		my @stack;

		while (($pos = index $seq, "N", $offset) != -1) {
			unless ($init) {
				$st = $pos;
				$en = $pos;
				$init = 1;
			}

			if ($en < ($pos - 1)) {
				push @stack => [$st, $en];
				$len += $en - $st + 1;
				$st = $pos;
			}

			$en = $pos;
			$offset = $pos + 1;
		}

		if ($init) {
			push @stack => [$st, $en];
			$len += $en - $st + 1;
		}

		my $read_len = $class eq 'App::Sandy::Seq::SingleEnd'
			? $self->seq->read_mean
			: $self->seq->fragment_mean;

		if ($len > ($seq_len - $read_len)) {
			log_msg ":: Seqid '$id' has too much NNN and will not be used";
			delete $indexed_fasta->{$id};
			push @blacklist_id => $id;
		} else {
			$fasta_blacklist{$id} = \@stack;
		}
	}

	unless (%$indexed_fasta) {
		die sprintf "Fasta file '%s' has too much NNN regions\n" => $self->fasta_file;
	}

	# If fasta_rtree has entries
	unless ($self->_has_no_fasta_rtree) {
		# Remove no valid entries from id -> pid relation
		$self->_delete_fasta_rtree(@blacklist_id) if @blacklist_id;
	}

	return \%fasta_blacklist;
}

sub _retrieve_expression_matrix {
	my $self = shift;

	my $expression = App::Sandy::DB::Handle::Expression->new;
	my $indexed_file = $expression->retrievedb($self->expression_matrix);

	# Remove transcript id version
	for my $id (keys %$indexed_file) {
		my $old_id = $id;
		$self->_clean_fasta_id(\$id);
		$indexed_file->{$id} = delete $indexed_file->{$old_id}
			if $old_id ne $id;
	}

	return $indexed_file;
}

sub _build_seqid_raffle {
	my $self = shift;

	# Get the piece table
	my $piece_table = $self->_piece_table;

	# The builded function
	my $seqid_sub;

	if ($self->seqid_weight eq 'same') {
		my ($keys, $weights) = $self->_populate_key_weight($piece_table, sub { 1 });

		# If weight == 1 means that there are 2 keys for
		# the same seq_id.
		# If weight == 2 means that there is only one key
		# for the seq_id, so I double that key
		for (my $i = 0; $i < @$weights; $i++) {
			if ($weights->[$i] > 1) {
				push @$keys => $keys->[$i];
			}

lib/App/Sandy/Simulator.pm view on Meta::CPAN

	# Build piece table
	my %piece_table;

	# Let's construct the piece_table
	log_msg ":: Build piece table ...";

	while (my ($seq_id, $fasta_h) = each %$indexed_fasta) {
		my $seq = \$fasta_h->{seq};
		my $std_seq_id = $self->_get_seqname($seq_id);

		# Initialize piece tables for $seq_id ref
		$piece_table{$seq_id}{ref}{table} = App::Sandy::PieceTable->new(orig => $seq);

		# If there is indexed_snv for seq_id, then construct the piece table with it
		if (defined $indexed_snv && defined $indexed_snv->{$std_seq_id}) {
			my $snvs = $indexed_snv->{$std_seq_id};

			# Filter only the homozygotic snvs to feed reference seq_id
			my @snvs_homo = grep { $_->{plo} eq 'HO' } @$snvs;

			if (@snvs_homo) {
				# Populate reference seq_id
				$self->_populate_piece_table($piece_table{$seq_id}{ref}{table}, \@snvs_homo);
			}

			# Initialize piece tables for $seq_id alt
			$piece_table{$seq_id}{alt}{table} = App::Sandy::PieceTable->new(orig => $seq);

			# Populate alternative seq_id
			$self->_populate_piece_table($piece_table{$seq_id}{alt}{table}, $snvs);
		}
	}

	# Initialize the logical offsets and valodate the
	# new size due to the genomic variation

	my @blacklist;

	for my $seq_id (keys %piece_table) {
		my $type_h = delete $piece_table{$seq_id};

		for my $type (keys %$type_h) {
			my $table_h = delete $type_h->{$type};
			my $table = $table_h->{table};

			# Initialize the logical offset
			$table->calculate_logical_offset;

			# Get the new size
			my $new_size = $table->logical_len;

			unless ($self->truncate) {
				my $class = ref $self->seq;

				if ($class eq 'App::Sandy::Seq::SingleEnd') {
					if ($new_size < $self->seq->read_mean) {
						log_msg ":: Skip '$seq_id:$type': So many deletions resulted in a sequence lesser than the required read-mean";
						next;
					}
				} elsif ($class eq 'App::Sandy::Seq::PairedEnd') {
					if ($new_size < $self->seq->fragment_mean) {
						log_msg ":: Skip '$seq_id:$type': So many deletions resulted in a sequence lesser than the required fragment mean";
						next;
					}
				} else {
					die "No valid options for 'seq'";
				}
			}

			# If all's right
			$table_h->{size} = $new_size;
			$type_h->{$type} = $table_h;
		}

		# if there is at least one type,
		# then return it to the piece_table
		if (%$type_h) {
			$piece_table{$seq_id} = $type_h;

		# else, just remove it!
		} else {
			push @blacklist => $seq_id;
		}
	}

	unless (%piece_table) {
		die "All fasta entries were removed due to deletions. ",
			"Please, verify the genomic variation '$genomic_variation'\n";
	}

	# If fasta_rtree has entries
	unless ($self->_has_no_fasta_rtree) {
		# Remove no valid entries from id -> pid relation
		$self->_delete_fasta_rtree(@blacklist) if @blacklist;
	}

	# Make the id -> pid relationship
	$self->_populate_fasta_tree;

	# HASH -> SEQ_ID -> @(REF @ALT) -> @(TABLE SIZE)
	return \%piece_table;
}

sub _populate_piece_table {
	my ($self, $table, $snvs) = @_;

	for my $snv (@$snvs) {
		# If there is an ID, make sure that it is not a comma, colon
		# separated list. Else, make sure to keep the ref/alt length
		# to max 25+25+1=51
		my $annot = defined $snv->{id} && $snv->{id} ne '.'
			? sprintf "%d:%s" => $snv->{pos} + 1, (split(/[,;]/, $snv->{id}))[0]
			: sprintf "%d:%.25s/%.25s" => $snv->{pos} + 1, $snv->{ref}, $snv->{alt};

		# Insertion
		if ($snv->{ref} eq '-') {
			$table->insert(\$snv->{alt}, $snv->{pos}, $annot);

		# Deletion
		} elsif ($snv->{alt} eq '-') {
			$table->delete($snv->{pos}, length $snv->{ref}, $annot);

lib/App/Sandy/Simulator.pm view on Meta::CPAN

			close $fhs[$fh_idx];
		}

		# If it is a bam and it is the last loop, then
		# write a eof marker
		if ($output_format eq 'bam' && $tid == $number_of_jobs) {
			$self->gen_eof_marker($files_t[0]);
		}

		# Child exit
		log_msg "  => Job $tid is finished";
		$pm->finish(0, \%counter);
	}

	# Back to parent
	# Interlace parent/child(s) processes
	my $sig = App::Sandy::InterlaceProcesses->new(foreign_pid => \@child_pid);
	$pm->wait_all_children;

	if ($sig->signal_catched) {
		log_msg ":: Termination signal received!";
	}

	log_msg ":: Saving the work ...";

	# Concatenate all temporary files
	log_msg ":: Concatenate all temporary files";

	# Save time. Rename tmp_file (1,2)
	for my $file (@{ $files{$file_class} }) {
		my $tmp = shift @tmp_files;
		log_msg "  => Concatenating $tmp to $file ...";
		rename $tmp => $file
			or die "Cannot create '$file': $!\n";
	}

	# Append to renamed tmp files
	my @fh = map { $self->with_open_a($_) } @{ $files{$file_class} };

	for my $i (0..$#tmp_files) {
		my $fh_idx = $i % scalar @fh;

		log_msg "  => Concatenating $tmp_files[$i] to $files{$file_class}[$fh_idx] ...";
		cat $tmp_files[$i] => $fh[$fh_idx]
			or die "Cannot concatenate $tmp_files[$i] to $files{$file_class}[$fh_idx]: $!\n";

		# Clean up the mess
		unlink $tmp_files[$i]
			or die "Cannot remove temporary file '$tmp_files[$i]': $!\n";
	}

	# Close files
	log_msg ":: Writing and closing output file: @{ $files{$file_class} }";
	for my $fh_idx (0..$#fh) {
		close $fh[$fh_idx]
			or die "Cannot write file $files{$file_class}[$fh_idx]: $!\n";
	}

	if ($self->count_loops_by eq 'number-of-reads') {
		# It is necessary to correct the abundance according to
		# fragment sequencing end
		my $count_factor = ref($self->seq) eq 'App::Sandy::Seq::PairedEnd'
			? 2
			: 1;

		# Save transcripts
		log_msg ":: Saving transcripts count";
		my $fh = $self->with_open_w($count_file{transcripts}, 0);

		log_msg "  => Writing counts to $count_file{transcripts} ...";
		for my $id (sort keys %counters) {
			printf {$fh} "%s\t%d\n" => $id,
				int($counters{$id} / $count_factor);
		}

		# Close transcripts file
		log_msg ":; Writing and closing $count_file{transcripts} ...";
		close $fh
			or die "Cannot write file $count_file{transcripts}: $!\n";

		# Calculate 'gene' like expression
		my $parent_count = $self->_calculate_parent_count(\%counters);

		if (%$parent_count) {
			# Save genes
			log_msg ":: Saving genes count";
			my $fh = $self->with_open_w($count_file{genes}, 0);

			log_msg "  => Writing counts to $count_file{genes} ...";
			for my $id (sort keys %$parent_count) {
				printf {$fh} "%s\t%d\n" => $id,
					int($parent_count->{$id} / $count_factor);
			}

			# Close genes file
			log_msg ":; Writing and closing $count_file{genes} ...";
			close $fh
				or die "Cannot write file $count_file{genes}: $!\n";
		}
	} else {
		# Save coverage
		log_msg ":: Saving coverage count";
		my $fh = $self->with_open_w($count_file{coverage}, 0);

		log_msg "  => Writing counts to $count_file{coverage} ...";
		for my $id (sort keys %counters) {
			printf {$fh} "%s\t%d\n" => $id, $counters{$id};
		}

		# Close coverage file
		log_msg ":; Writing and closing $count_file{coverage} ...";
		close $fh
			or die "Cannot write file $count_file{coverage}: $!\n";
	}
}

__END__

=pod

=encoding UTF-8

( run in 0.820 second using v1.01-cache-2.11-cpan-b16cb0d3907 )