File Coverage

lib/App/SimulateReads/Read/PairedEnd.pm
Criterion Covered Total %
statement 31 32 96.8
branch 7 8 87.5
condition 2 3 66.6
subroutine 6 6 100.0
pod 0 2 0.0
total 46 51 90.2


line stmt bran cond sub pod time code
1             package App::SimulateReads::Read::PairedEnd;
2             # ABSTRACT: App::SimulateReads::Read subclass for simulate paired-end reads.
3              
4 6     6   36 use App::SimulateReads::Base 'class';
  6         12  
  6         82  
5 6     6   2647 use Math::Random 'random_normal';
  6         26200  
  6         556  
6              
7             extends 'App::SimulateReads::Read';
8              
9             our $VERSION = '0.06'; # VERSION
10              
11             use constant {
12 6         2506 NUM_TRIES => 1000
13 6     6   42 };
  6         12  
14              
15             has 'fragment_mean' => (
16             is => 'ro',
17             isa => 'My:IntGt0',
18             required => 1
19             );
20              
21             has 'fragment_stdd' => (
22             is => 'ro',
23             isa => 'My:IntGe0',
24             required => 1
25             );
26              
27             sub BUILD {
28 32     32 0 62 my $self = shift;
29 32 100       892 unless (($self->fragment_mean - $self->fragment_stdd) >= $self->read_size) {
30 1         32 croak sprintf "fragment_mean (%d) minus fragment_stdd (%d) must be greater or equal to read_size (%d)\n"
31             => $self->fragment_mean, $self->fragment_stdd, $self->read_size;
32             }
33             }
34              
35             sub gen_read {
36 701     701 0 11096 my ($self, $seq_ref, $seq_size, $is_leader) = @_;
37              
38 701 100       19026 if ($seq_size < $self->fragment_mean) {
39 1         28 croak sprintf "seq_size (%d) must be greater or equal to fragment_mean (%d)\n"
40             => $seq_size, $self->fragment_mean;
41             }
42              
43 700         1017 my $fragment_size = 0;
44 700         880 my $random_tries = 0;
45              
46 700   66     17250 until (($fragment_size <= $seq_size) && ($fragment_size >= $self->read_size)) {
47             # seq_size must be greater or equal to fragment_size and
48             # fragment_size must be greater or equal to read_size
49             # As fragment_size is randomly calculated, try out NUM_TRIES times
50 700 50       1364 if (++$random_tries > NUM_TRIES) {
51 0         0 croak sprintf
52             "So many tries to calculate a fragment. the constraints were not met:\n" .
53             "fragment_size <= seq_size (%d) and fragment_size >= read_size (%d)\n"
54             => $seq_size, $self->read_size;
55             }
56              
57 700         1472 $fragment_size = $self->_random_half_normal;
58             }
59              
60 700         2094 my ($fragment_ref, $fragment_pos) = $self->subseq_rand($seq_ref, $seq_size, $fragment_size);
61              
62 700         16552 my $read1_ref = $self->subseq($fragment_ref, $fragment_size, $self->read_size, 0);
63 700         15832 $self->update_count_base($self->read_size);
64 700         1973 $self->insert_sequencing_error($read1_ref);
65              
66 700         16304 my $read2_ref = $self->subseq($fragment_ref, $fragment_size, $self->read_size, $fragment_size - $self->read_size);
67 700         1850 $self->reverse_complement($read2_ref);
68 700         16069 $self->update_count_base($self->read_size);
69 700         1753 $self->insert_sequencing_error($read2_ref);
70              
71 700 100       3007 return $is_leader ?
72             ($read1_ref, $read2_ref, $fragment_pos, $fragment_size) :
73             ($read2_ref, $read1_ref, $fragment_pos, $fragment_size);
74             }
75              
76             sub _random_half_normal {
77 700     700   1030 my $self = shift;
78 700         18189 return abs(int(random_normal(1, $self->fragment_mean, $self->fragment_stdd)));
79             }
80              
81             __END__
82              
83             =pod
84              
85             =encoding UTF-8
86              
87             =head1 NAME
88              
89             App::SimulateReads::Read::PairedEnd - App::SimulateReads::Read subclass for simulate paired-end reads.
90              
91             =head1 VERSION
92              
93             version 0.06
94              
95             =head1 AUTHOR
96              
97             Thiago L. A. Miller <tmiller@mochsl.org.br>
98              
99             =head1 COPYRIGHT AND LICENSE
100              
101             This software is Copyright (c) 2017 by Teaching and Research Institute from Sírio-Libanês Hospital.
102              
103             This is free software, licensed under:
104              
105             The GNU General Public License, Version 3, June 2007
106              
107             =cut