Mercurial > repos > cpt > cpt_psm_prep

comparison lib/CPT/Analysis/PAUSE.pm @ 1:d724f34e671d draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload commit 94b0cd1fff0826c6db3e7dc0c91c0c5a8be8bb0c
author cpt
date Mon, 05 Jun 2023 02:50:07 +0000
parents
children
comparison
equal deleted inserted replaced
0:e4de0a0e90c8 1:d724f34e671d
1 package CPT::Analysis::PAUSE;
2
3 # ABSTRACT: Library for use in PAUSE analysis
4 use strict;
5 use warnings;
6 use Moose;
7 use List::Util qw(sum);
8 use Statistics::Descriptive;
9
10 sub max ($$) { shift; $_[ $_[0] < $_[1] ] }
11 sub min ($$) { shift; $_[ $_[0] > $_[1] ] }
12
13 sub derivative {
14 my ( $self, $data_ref ) = @_;
15 my @data = @{$data_ref};
16 my @new_data;
17 foreach ( my $i = 0 ; $i < scalar(@data) - 1 ; $i++ ) {
18 $new_data[ $i + 1 ] = $data[ $i + 1 ] - $data[$i];
19 }
20 return \@new_data;
21 }
22
23 sub find_peaks {
24 my ( $self, %data ) = @_;
25 use IPC::Run3;
26 use File::Temp qw/tempfile/;
27
28 # Store to CSV File
29 my @starts = @{ $data{data} };
30 my ( $fh0, $filename0 ) = tempfile('galaxy.pause.XXXXXXX');
31 printf $fh0 ( "%s,%s\n", 'position', 'count' );
32 for ( my $i = 0 ; $i < scalar(@starts) ; $i++ ) {
33 printf $fh0 "%d,%d\n", $i,
34 ( defined $starts[$i] ? $starts[$i] : 0 );
35 }
36 close($fh0);
37
38 my ( $fh, $filename ) = tempfile('galaxy.pause.XXXXXXX');
39 my @cmd = (
40 'Rscript', $data{location_of_rscript_file},
41 $filename0, $filename, $data{snr}
42 );
43 my ( $in, $out, $err );
44 run3 \@cmd, \$in, \$out, \$err;
45
46 # Read in R data
47 my @values;
48 while (<$fh>) {
49 chomp;
50 push( @values, $_ );
51 }
52 close($fh);
53
54 unlink($filename0);
55 unlink($filename);
56
57 return @values;
58 }
59
60 sub smooth {
61 my ( $self, $data_ref ) = @_;
62 my @data = @{$data_ref};
63 my @new_data;
64 my $length = scalar @data;
65 foreach ( my $i = 0 ; $i < $length ; $i++ ) {
66 my $avg =
67 sum( @data[ $i - 20 .. $i - 1, $i + 1 .. $i + 20 ] ) / 40;
68 $new_data[$i] = $avg;
69 }
70 return \@new_data;
71 }
72
73 sub histogram {
74 my ( $self, %data ) = @_;
75
76 my @coverage = @{ $data{data} };
77 my @return_coverage;
78 for ( my $i = 0 ; $i < scalar(@coverage) ; $i++ ) {
79 my $size = $coverage[$i];
80 unless ($size) { $size = 0 }
81 $return_coverage[$i] = [ $i, $size, "*" x $size ];
82 }
83 my %results = (
84 'Sheet1' => {
85 headers => [qw(Base Count Plot)],
86 data => \@return_coverage,
87 }
88 );
89 return %results;
90 }
91
92 sub getCoverageDensity {
93 my ( $self, %data ) = @_;
94
95 # Load the sam file
96 my $sam = Bio::DB::Sam->new(
97 -bam => $data{bam},
98 -fasta => $data{genome},
99 -autoindex => 1,
100 );
101
102 # Get all alignments to our indicated FASTA file
103 my @alignments = $sam->get_features_by_location(
104 -seq_id => $data{fasta_id},
105 -start => 1,
106 -end => $data{fasta_length}
107 );
108
109 # Set up some variables
110 my $coverage_density_max_value = 0;
111 my ( @coverage_density, @read_starts, @read_ends );
112
113 # including some for statistics
114 my $stat_start = Statistics::Descriptive::Sparse->new();
115 my $stat_end = Statistics::Descriptive::Sparse->new();
116
117 # Looping over alignments
118 for my $a (@alignments) {
119 my $start = $a->start;
120 my $end = $a->end;
121
122 # Increment the number of reads starting there
123 $read_starts[$start]++;
124 $read_ends[$end]++;
125
126 # And increment the coverage density
127 foreach ( $start .. $end ) {
128 $coverage_density[$_]++;
129 if ( $coverage_density[$_] >
130 $coverage_density_max_value )
131 {
132 $coverage_density_max_value =
133 $coverage_density[$_];
134 }
135 }
136 }
137 my @start_data_for_stats;
138 my @end_data_for_stats;
139 for ( my $i = 0 ; $i < $data{fasta_length} ; $i++ ) {
140 if ( $read_starts[$i] ) {
141 push( @start_data_for_stats, $read_starts[$i] );
142 }
143 if ( $read_ends[$i] ) {
144 push( @end_data_for_stats, $read_ends[$i] );
145 }
146 }
147 $stat_start->add_data(@start_data_for_stats);
148 $stat_end->add_data(@end_data_for_stats);
149
150 # Lots of data to return
151 use CPT::Analysis::PAUSE::ParsedSam;
152 my $psam = CPT::Analysis::PAUSE::ParsedSam->new(
153 coverage_density => \@coverage_density,
154 read_starts => \@read_starts,
155 read_ends => \@read_ends,
156 max => $coverage_density_max_value,
157 stats_start_max => $stat_start->max(),
158 stats_end_max => $stat_end->max(),
159 stats_start_mean => $stat_start->mean(),
160 stats_end_mean => $stat_end->mean(),
161 stats_start_standard_deviation =>
162 $stat_start->standard_deviation(),
163 stats_end_standard_deviation => $stat_end->standard_deviation(),
164 );
165 return $psam;
166 }
167
168 no Moose;
169 1;
170
171 __END__
172
173 =pod
174
175 =encoding UTF-8
176
177 =head1 NAME
178
179 CPT::Analysis::PAUSE - Library for use in PAUSE analysis
180
181 =head1 VERSION
182
183 version 1.96
184
185 =head1 AUTHOR
186
187 Eric Rasche <rasche.eric@yandex.ru>
188
189 =head1 COPYRIGHT AND LICENSE
190
191 This software is Copyright (c) 2014 by Eric Rasche.
192
193 This is free software, licensed under:
194
195 The GNU General Public License, Version 3, June 2007
196
197 =cut