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
|
# This demo code shows how to write a custom perl code and use it in
# the -i/-e filtering expressions.
#
# In this example we want to take variant consequences generated by `bcftools csq`,
# rank them by severity, filter by the most severe, and print the list using the
# following command:
#
# bcftools query -f '%CHROM:%POS \t %CSQ\n' -i 'perl:misc/demo-flt.pl; perl.severity(INFO/CSQ) > 10' test/perl-flt.vcf
#
# There can be multiple subroutines in the same script and they
# can be referenced from the command line by prefixing "perl.subroutine_name()"
#
sub severity
{
# Arbitrary number of arguments can be provided.
my (@args) = @_;
# Note that string arrays are passed to perl in the form of a single
# comma-separated string, but numeric arrays must be dereferenced, for
# example as shown in this example:
#
# for my $arg (@args)
# {
# if ( ref($arg) eq 'ARRAY' ) { print "array: ".join(',',@$arg)."\n"; }
# else { print "scalar: $arg\n"; }
# }
#
# In our case there should be only one parameter passed to the subroutine;
# check if this is the case
if ( scalar @args != 1 ) { error("Incorrect use, expected one argument, got ".scalar @args."!\n"); }
# Create a lookup table from consequences to an arbitrary severity score
my %severity =
(
"intergenic" => 1,
"intron" => 2,
"non_coding" => 3,
"5_prime_utr" => 4,
"3_prime_utr" => 5,
"stop_retained" => 6,
"synonymous" => 7,
"coding_sequence" => 8,
"missense" => 9,
"splice_region" => 10,
"inframe_altering" => 11,
"inframe_deletion" => 12,
"inframe_insertion" => 13,
"splice_acceptor" => 14,
"splice_donor" => 15,
"stop_lost" => 16,
"stop_gained" => 17,
"start_lost" => 18,
"frameshift" => 19,
);
# Split multiple consequences into an array
my @csq = split(/,/, $args[0]);
# Find the most severe consequence. The consequence string may look like this:
# inframe_deletion|ABC|ENST00000000002|protein_coding|-|5YV>5T|144TAC>T+148TA>T
my $max = 0;
for my $csq (@csq)
{
my @fields = split(/\|/, $csq);
my $string = $fields[0];
my $score = exists($severity{$string}) ? $severity{$string} : 0;
if ( $max < $score ) { $max = $score; }
}
return $max;
}
sub error
{
my (@msg) = @_;
print STDERR @msg;
exit 1;
}
|
