view release on metacpan or  search on metacpan

lib/Bio/ToolBox/db_helper/hts.pm  view on Meta::CPAN

					}
				}
				else {
					$scores = $coverage;
				}
			}
		}
	}

	## Return collected data
	if ( $param->[RETT] == 2 ) {
		return wantarray ? %pos2data : \%pos2data;
	}
	else {
		return wantarray ? @{$scores} : $scores;
	}
}

### Determine total number of alignments in a bam file
sub sum_total_bam_alignments {

	# Passed arguments;
	my $sam_file = shift;
	my $cpu      = shift || 2;    # number of forks to execute in parallel
	$cpu = 1 unless ($parallel);
	unless ($sam_file) {
		carp 'ERROR: no Bam file or bam db object passed!';
		return;
	}

	# Open Bam file if necessary
	my $bam;
	my $bam_ref = ref $sam_file;
	if ( $bam_ref =~ /Bio::DB::HTS/x ) {

		# we have an opened sam db object
		$bam = $sam_file;
	}
	else {
		# we have a name of a sam file
		$bam = open_bam_db($sam_file);
		return unless ($bam);
	}

	# prepare the alignment iterator callback
	my $callback = sub {
		my ( $a, $n ) = @_;

		# check alignment flags
		my $flag = $a->flag;
		return if $flag & 0xf04;    # UNMAP,SECONDARY,QCFAIL,DUP,SUPPLEMENTARY

		# check mapping quality based on global variable
		return if $a->qual < $MAPQ;

		# count properly paired alignments only once by skipping reverse alignment
		if ( $flag & 0x01 and $flag & 0x2 ) {
			return if ( $flag & 0x10 );
		}

		# count this fragment
		${$n}++;
	};

	# prepare the counting subroutine
	my $counter = sub {
		my $tid    = shift;
		my $number = 0;
		$bam->hts_index->fetch( $bam->hts_file, $tid, 0, $bam->target_len($tid),
			$callback, \$number );
		return $number;
	};

	# Count the alignments on each chromosome
	my $total_read_number = 0;
	if ( $cpu > 1 ) {

		# count each chromosome in multiple parallel threads to speed things up

		# generate relatively equal lists of chromosome for each process based on length
		my @chromosomes = map { $_->[0] }
			sort { $b->[1] <=> $a->[1] }
			map { [ $_, $bam->target_len($_) ] } ( 0 .. $bam->n_targets - 1 );
		my @list;    # array of arrays, [process][chromosome id]
		my $i = 1;
		while (@chromosomes) {
			push @{ $list[$i] }, shift @chromosomes;
			$i++;
			$i = 1 if $i > $cpu;
		}

		# we will use Parallel ForkManager for convenience
		my $pm = Parallel::ForkManager->new($cpu);
		$pm->run_on_finish(
			sub {
				my ( $pid, $exit_code, $ident, $exit_signal, $core_dump, $count ) =
					@_;
				$total_read_number += ${$count};
			}
		);

		# Count the chromosomes in parallel processess
		foreach my $n ( 1 .. $cpu ) {
			$pm->start and next;

			### In child
			$bam->clone;    # to make it fork safe
			my $count = 0;
			foreach ( @{ $list[$n] } ) {

				# count each chromosome in this process list
				$count += &{$counter}($_);
			}
			$pm->finish( 0, \$count );
		}
		$pm->wait_all_children;
	}

	else {
		# loop through all the chromosomes in one execution thread
		for my $tid ( 0 .. $bam->n_targets - 1 ) {

lib/Bio/ToolBox/db_helper/hts.pm  view on Meta::CPAN

if not, generate one. The L<Bio::DB::HTS> module uses the samtools 
style index extension, F<.bam.bai>, as opposed to the Picard style 
extension, F<.bai>. If a F<.bai> index is present, it will copy the 
file as F<.bam.bai> index. Unfortunately, a F<.bai> index cannot be 
used directly.

This method is called automatically prior to opening a bam database. 

=item write_new_bam_file

This subroutine will open a new empty Bam file. Pass the name of the 
new file as the argument. It will return a L<Bio::DB::HTSfile> object to 
which you can write a header followed by alignments. Be sure you know 
what to do before using this method! 

=item collect_bam_scores

This subroutine will collect only the data values from a binary bam file 
for the specified database region. The positional information of the 
scores is not retained.

Collected data values may be restricted to strand by specifying the desired 
strandedness (sense, antisense, or all), 
depending on the method of data collection. Collecting scores, or basepair 
coverage of alignments over the region of interest, does not currently support 
stranded data collection (as of this writing). However, enumerating 
alignments (count method) and collecting alignment lengths do support 
stranded data collection. Alignments are checked to see whether their midpoint 
is within the search interval before counting or length collected. 

As of version 1.30, paired-end bam files are properly handled with regards 
to strand; Strand is determined by the orientation of the first mate. However, 
pairs are still counted as two alignments, not one. To avoid this, use the 
value_type of 'ncount' and count the number of unique alignment names. 
(Previous versions treated all paired-end alignments as single-end alignments, 
severely limiting usefulness.)

The subroutine is passed a parameter array reference. See 
L</"Data Collection Parameters Reference"> below for details.

The subroutine returns an array or array reference of the requested dataset 
values found within the region of interest. 

=item collect_bam_position_scores

This subroutine will collect the score values from a binary bam file 
for the specified database region keyed by position. 

The subroutine is passed a parameter array reference. See 
L</"Data Collection Parameters Reference"> below for details.

The subroutine returns a hash or hash reference of the defined dataset values 
found within the region of interest keyed by position. The feature midpoint 
is used as the key position. When multiple features are found at the same 
position, a simple mean (for length data methods) or sum 
(for count methods) is returned. The ncount value type is not supported 
with positioned scores.

=item sum_total_bam_alignments

This subroutine will sum the total number of mapped fragments 
in a bam file. Properly paired alignments are counted once, and
supplemental, secondary, and marked duplicate alignments are ignored.
Alignments may be filtered by mapping quality. Use the
L<Bio::ToolBox::db_helper> C<use_minimum_mapq()> function to set
the minimum allowed mapping quality.
 
Pass the subroutine one or optionally two arguments in the following order. 

=over 4

=item 1. Bam file path or object 

The name of the Bam file which should be counted. Alternatively, an 
opened L<Bio::DB::HTS> object may also be given. Required.

=item 2. Number of forks (integer)

Optionally pass the number of parallel processes to execute 
when counting alignments. Walking through a Bam file is 
time consuming but can be easily parallelized. The module 
L<Parallel::ForkManager> is required, and the default is a 
conservative two processes when it is installed.

=back
       
The subroutine will return the number of alignments.

=back

=head2 Data Collection Parameters Reference

The data collection subroutines are passed an array reference of parameters. 
The recommended  method for data collection is to use the
L<Bio::ToolBox::db_helper/get_segment_score> method. 

The parameters array reference includes these items:

=over 4

=item 1. chromosome

=item 1. start coordinate

=item 3. stop coordinate

=item 4. strand

Should be standard BioPerl representation: -1, 0, or 1.

=item 5. strandedness

A scalar value representing the desired strandedness of the data 
to be collected. Acceptable values include "sense", "antisense", 
or "all". Only those scores which match the indicated 
strandedness are collected.

=item 6. score method

Acceptable values include score, count, pcount, and ncount.