view release on metacpan or  search on metacpan

scripts/get_datasets.pl  view on Meta::CPAN

			if ( $formatter and looks_like_number($score) ) {
				$score = sprintf( $formatter, $score );
			}
			$row->value( $index, $score );
		}
	);
}

sub get_fractionated_dataset {
	my ( $dataset, $index ) = @_;

	# collect the scores from the dataset for this index
	$Data->iterate(
		sub {
			my $row = shift;

			# calculate length
			my $length = $row->length;

			# calculate new fractional start and stop positions
			# the fraction depends on the length
			# this depends on both feature orientation and the
			# relative position requested
			my $relative_start = int( ( $length * $fstart ) + 0.5 );
			my $relative_stop  = int( ( $length * $fstop ) + 0.5 );
			my ( $start, $stop );
			if ( $length >= $limit ) {

				# length exceeds our minimum limit
				# we can take a fractional length

				if ( $position == 5 and $row->strand >= 0 ) {

					# 5' end of forward strand
					$start = $row->start + $relative_start;
					$stop  = $row->start + $relative_stop;
				}
				elsif ( $position == 5 and $row->strand < 0 ) {

					# 5' end of reverse strand
					$start = $row->end - $relative_stop;
					$stop  = $row->end - $relative_start;
				}
				elsif ( $position == 3 and $row->strand >= 0 ) {

					# 3' end of forward strand
					$start = $row->end + $relative_start;
					$stop  = $row->end + $relative_stop;
				}
				elsif ( $position == 3 and $row->strand < 0 ) {

					# 3' end of reverse strand
					$start = $row->start - $relative_stop;
					$stop  = $row->start - $relative_start;
				}

				# midpoint is not accepted
			}
			else {
				# length doesn't meet minimum limit
				# simply take the whole fragment
				$start = $row->start;
				$stop  = $row->end;
			}

			# now collect score
			my $score = $row->get_score(
				'start'    => $start,
				'stop'     => $stop,
				'db'       => $ddb,
				'dataset'  => $dataset,
				'method'   => $method,
				'stranded' => $stranded,
			);
			if ( $formatter and looks_like_number($score) ) {
				$score = sprintf( $formatter, $score );
			}
			$row->value( $index, $score );
		}
	);
}

# subroutine to record the metadata for a dataset
sub add_new_dataset {
	my $dataset = shift;

	# generate column name
	my $column_name = simplify_dataset_name($dataset);

	# add new column
	my $index = $Data->add_column($column_name);

	# update metadata
	$Data->metadata( $index, 'dataset',        $dataset );
	$Data->metadata( $index, 'method',         $method );
	$Data->metadata( $index, 'strand',         $stranded );
	$Data->metadata( $index, 'extend',         $extend )     if defined $extend;
	$Data->metadata( $index, 'start',          $start_adj )  if defined $start_adj;
	$Data->metadata( $index, 'stop',           $stop_adj )   if defined $stop_adj;
	$Data->metadata( $index, 'fstart',         $fstart )     if defined $fstart;
	$Data->metadata( $index, 'fstop',          $fstop )      if defined $fstop;
	$Data->metadata( $index, 'limit',          $limit )      if defined $limit;
	$Data->metadata( $index, 'subfeature',     $subfeature ) if $subfeature;
	$Data->metadata( $index, 'forced_strand',  'yes' )       if $set_strand;
	$Data->metadata( $index, 'decimal_format', $format )     if defined $format;
	$Data->metadata( $index, 'total_reads',    $dataset2sum{$dataset} )
		if exists $dataset2sum{$dataset};

	if ( $position == 5 ) {
		$Data->metadata( $index, 'relative_position', "5'end" );
	}
	if ( $position == 3 ) {
		$Data->metadata( $index, 'relative_position', "3'end" );
	}
	elsif ( $position == 4 ) {
		$Data->metadata( $index, 'relative_position', 'middle' );
	}
	elsif ( $position == 9 ) {
		$Data->metadata( $index, 'relative_position', 'summit' );
	}

scripts/get_datasets.pl  view on Meta::CPAN

Prevent input annotation files from being automatically parsed into sequence 
features. Coordinates will be used as is and new data columns will be appended 
to the input file. 

=item --version

Print the version number.

=item --help

Display the POD documentation for this program.

=back

=head1 DESCRIPTION

This program will collect dataset values from a variety of sources, including 
features in a BioPerl Bio::DB::SeqFeature::Store database, binary wig files 
F<.wib> loaded in a database using Bio::Graphics::Wiggle, bigWig files, 
bigBed files, Bam alignment files, or a Bio::DB::BigWigSet database. 

The values are collected for a list of known database features (genes, 
transcripts, etc.) or genomic regions (defined by chromosome, start, and 
stop). The list may be provided as an input file or generated as a new 
list from a database. Output data files may be reloaded for additional 
data collection.

At each feature or interval, multiple data points within the genomic segment 
are combined statistically and reported as a single value for the feature. 
The method for combining datapoints may be specified; the default method is 
the mean of all datapoints.

The coordinates of the features may be adjusted in numerous ways, including 
specifying a specific relative start and stop, a fractional start and stop, 
an extension to both start and stop, and specifying the relative position 
(5' or 3' or midpoint).

Stranded data may be collected, if the dataset supports stranded information. 
Also, two or more datasets may be combined and treated as one. Note that 
collecting stranded data may significantly slow down data collection.

=head1 EXAMPLES

These are some examples of some common scenarios for collecting data.

=over 4

=item Simple mean scores

You want to collect the mean score from a bigWig file for each feature 
in a BED file of intervals.

  get_datasets.pl --in input.bed --data scores.bw

=item Collect normalized counts

You want to collect normalized read counts from multiple Bam files 
for each feature in a BED file. This will count alignment names (safe for 
paired-end alignments) over the intervals, and transform to Fragments (Reads) 
Per Million, a depth-normalizing function based on the total number of 
fragments counted in each dataset. 

  get_datasets.pl --in input.bed --method ncount --fpkm region *.bam 

=item Collect stranded RNASeq data

You have stranded RNASeq data, and you would like to determine the 
expression level for all genes from an annotation file. Use the 
C<ncount> method to count alignment names to avoid double counting 
alignments split over multiple exons.
  
  get_datasets.pl --in annotation.gtf --feature transcript --subfeature exon \
  --strand sense --method ncount --out expression.txt *.bam

=item Restrict to specific region

You have ChIPSeq enrichment scores in a bigWig file and you now want 
to score just the transcription start site of known transcripts in a 
L<Bio::DB::SeqFeature::Store> annotation database. Here you will 
restrict to 500 bp flanking the TSS.
  
  get_datasets.pl --db annotation.sqlite --feature mRNA --start=-500 \
  --stop=500 --pos 5 --data scores.bw --out tss_scores.txt

=item Avoid first and last 1 Kb of each interval

  get_datasets.pl --in file.bed --start=1000 \
  --stop=-1000 --pos 53 --data scores.bw --out file_scores.txt


=item Count intervals

You have identified all possible transcription factor binding sites in 
the genome and put them in a bigBed file. Now you want to count how 
many exist in each upstream region of each gene.
  
  get_datasets.pl --db annotation.gtf --feature gene --start=-5000 \
  --stop=0 --data tfbs.bb --method count --out tfbs_sums.txt

=item Many datasets at once

While you can provide multiple data source files as a space-delimited 
list at the end of the command, you can also treat a folder or directory 
of bigWig files as a special database, known as a BigWigSet. Each file 
becomes a database feature, and you can interactively choose one or more 
from which to collect. Each dataset is appended to the input file as a new 
column. Provide the folder as a data database C<--ddb> option.
  
  get_datasets.pl --in input.txt --ddb /path/to/bigwigset

=item Stranded BigWig data

You can generate stranded RNASeq coverage from a Bam file using the 
BioToolBox script bam2wig.pl, which yields rnaseq_f.bw and rnaseq_r.bw 
files. These are automatically interpreted as stranded datasets in a 
BigWigSet folder context; see above.
  
  get_datasets.pl --in input.txt --strand sense \
  --ddb /path/to/rnaseq/bigwigset 

=item Binned coverage across the genome