1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
|
#! perl
use strict;
use warnings;
use Test::More;
use lib '.';
use t::TestUtils;
use Grinder;
my ($factory, $lib, $nof_libs, $nof_reads, $read);
# Multiple shotgun libraries
ok $factory = Grinder->new(
-reference_file => data('shotgun_database.fa'),
-read_dist => 48 ,
-num_libraries => 4 ,
-total_reads => 100 ,
), 'Multiple shotgun libraries';
$nof_libs = 0;
while ( $lib = $factory->next_lib ) {
$nof_libs++;
$nof_reads = 0;
while ( $read = $factory->next_read ) {
$nof_reads++;
ok_read($read, undef, $nof_reads, $nof_libs);
};
is $nof_reads, 100;
}
is $nof_libs, 4;
# Multiple mate pair libraries
ok $factory = Grinder->new(
-reference_file => data('shotgun_database.fa'),
-total_reads => 100 ,
-read_dist => 48 ,
-num_libraries => 4 ,
-insert_dist => 250 ,
), 'Multiple mate pair libraries';
$nof_libs = 0;
while ( $lib = $factory->next_lib ) {
$nof_libs++;
$nof_reads = 0;
while ( $read = $factory->next_read ) {
$nof_reads++;
ok_mate($read, undef, $nof_reads, $nof_libs);
};
is $nof_reads, 100;
}
is $nof_libs, 4;
done_testing();
sub ok_read {
my ($read, $req_strand, $nof_reads, $nof_libs) = @_;
isa_ok $read, 'Bio::Seq::SimulatedRead';
my $source = $read->reference->id;
my $strand = $read->strand;
if (not defined $req_strand) {
$req_strand = $strand;
} else {
is $strand, $req_strand;
}
my $letters;
if ( $source eq 'seq1' ) {
$letters = 'a';
} elsif ( $source eq 'seq2' ) {
$letters = 'c';
} elsif ( $source eq 'seq3' ) {
$letters = 'g';
} elsif ( $source eq 'seq4' ) {
$letters = 't';
} elsif ( $source eq 'seq5' ) {
$letters = 'atg';
}
if ( $req_strand == -1 ) { # Take the reverse complement
$letters = Bio::PrimarySeq->new( -seq => $letters )->revcom->seq;
}
like $read->seq, qr/[$letters]+/;
is $read->id, $nof_libs.'_'.$nof_reads;
is $read->length, 48;
}
sub ok_mate {
my ($read, $req_strand, $nof_reads, $nof_libs) = @_;
isa_ok $read, 'Bio::Seq::SimulatedRead';
my $source = $read->reference->id;
my $strand = $read->strand;
if (not defined $req_strand) {
$req_strand = $strand;
} else {
is $strand, $req_strand;
}
my $letters;
if ( $source eq 'seq1' ) {
$letters = 'a';
} elsif ( $source eq 'seq2' ) {
$letters = 'c';
} elsif ( $source eq 'seq3' ) {
$letters = 'g';
} elsif ( $source eq 'seq4' ) {
$letters = 't';
} elsif ( $source eq 'seq5' ) {
$letters = 'atg';
}
if ( $req_strand == -1 ) { # Take the reverse complement
$letters = Bio::PrimarySeq->new( -seq => $letters )->revcom->seq;
};
like $read->seq, qr/[$letters]+/;
my $id = $nof_libs.'_'.round($nof_reads/2).'/'.($nof_reads%2?1:2);
is $read->id, $id;
is $read->length, 48;
}
|
