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
126
127
128
|
#!/usr/bin/env perl
use strict;
use warnings;
use lib ("/usr/lib/trinityrnaseq/PerlLib");
use DelimParser;
use Getopt::Long qw(:config posix_default no_ignore_case bundling pass_through);
my $pair_stats_file;
my $max_cov;
my $min_cov=1;
my $max_CV;
my $usage = <<__EOUSAGE__;
#############################################################
#
# Required:
#
# --stats_file <string> : pairs.stats.sorted
#
# --max_cov <int> : maximum coverage
#
# --min_cov <int> : minimum coverage
#
# --max_CV <int> : maximum coeff. var.
#
#
#############################################################
__EOUSAGE__
;
my $help_flag;
&GetOptions ( 'h' => \$help_flag,
'stats_file=s' => \$pair_stats_file,
'max_cov=i' => \$max_cov,
'min_cov=i' => \$min_cov,
'max_CV=i' => \$max_CV,
);
unless ($pair_stats_file &&
$max_cov &&
defined($min_cov) &&
defined($max_CV)) {
die $usage;
}
main: {
srand(12345);
my $count_aberrant_and_discarded = 0;
my $count_selected = 0;
my $count_total = 0;
my $count_below_min_cov = 0;
open (my $fh, $pair_stats_file) or die $!;
my $delim_parser = new DelimParser::Reader($fh, "\t");
while (my $row = $delim_parser->get_row()) {
$count_total++;
my $core_acc = $row->{acc};
$core_acc =~ s|/[12]$||;
my $med_cov = ($row->{median_cov});
if ($med_cov < $min_cov) {
$count_below_min_cov++;
next;
}
my $sd = $row->{stdev};
my $u = $row->{mean_cov};
if ($u <= 0) {
$count_aberrant_and_discarded++;
next;
}
my $cv = $sd/$u;
if ($cv > $max_CV) {
$count_aberrant_and_discarded++;
next;
}
if (rand(1) <= $max_cov/$med_cov) {
print "$core_acc\n";
$count_selected++;
}
}
close $fh;
unless ($count_total) {
die "Error, no reads made it to the normalization process... ";
}
print STDERR "$count_selected / $count_total = " . sprintf("%.2f", $count_selected/$count_total*100) . "% reads selected during normalization.\n";
print STDERR "$count_aberrant_and_discarded / $count_total = " . sprintf("%.2f", $count_aberrant_and_discarded/$count_total*100) . "% reads discarded as likely aberrant based on coverage profiles.\n";
print STDERR "$count_below_min_cov / $count_total = " . sprintf("%.2f", $count_below_min_cov/$count_total*100) . "% reads discarded as below minimum coverage threshold=$min_cov\n";
exit(0);
}
|
