This file is indexed.

/usr/bin/maq_plot.pl is in maq 0.7.1-5.

This file is owned by root:root, with mode 0o755.

The actual contents of the file can be viewed below.

  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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
#!/usr/bin/perl -w

# Author: lh3

# This script will call "gnuplot" almost everywhere.

use strict;
use warnings;
use Getopt::Std;

&usage if (@ARGV < 1);

my $version = "0.1.1";
my $command = shift(@ARGV);
my %func = (abpair=>\&abpair, depth=>\&depth_plot, mixed=>\&mixed);
die("Unknown command \"$command\".\n") if (!defined($func{$command}));
&{$func{$command}}();
exit(0);

sub depth_plot
{
	my %opts;
	getopts('c:D:ax:', \%opts);
	die qq(Usage: maq_plot.pl depth [-c chr] [-D 100] [-a] [-x 2] <prefix> <cns2win-output>\n) if (@ARGV < 2); 
	my $chr = $opts{c};
	my $max_depth = $opts{D}? $opts{D} : 100;
	my $xtics = $opts{x}? $opts{x} : 2;
	my $back = 5;
	my $prefix = shift(@ARGV);
	my ($fh1, $fh2, $fh3, $fh, $fhin);
	open($fhin, $ARGV[0]) || die;
	my ($mean, $var, $n) = (0, 0, 0);
	my $col = (defined $opts{a})? 5 : 4;
	while (<$fhin>) {
		my @t = split;
		$chr = $t[0] if (!$chr);
		if ($chr eq $t[0] && $t[7] >= 1e-6) { # then, not N
			++$n;
			$mean += $t[$col];
		}
	}
	$mean /= $n;
	open($fh1, ">$prefix.depth") || die;
	open($fh2, ">$prefix.outlier") || die;
	open($fh3, ">$prefix.gc") || die;
	seek($fhin, 0, 0);
	$n = 0;
	while (<$fhin>) {
		my @t = split;
		if ($chr eq $t[0]) {
			if ($t[$col] < 2 * $mean && $t[$col] > 0.5 * $mean) {
				if ($t[7] >= 1e-6) {
					++$n;
					$var += ($t[$col] - $mean) * ($t[$col] - $mean);
				}
				$fh = ($t[$col] < 1.5 * $mean)? $fh1 : $fh2;
			} else { $fh = $fh2; }
			my $depth = ($t[$col] >= $max_depth-$back)? $max_depth - $back * (1.0 - rand()) : $t[$col];
			print $fh "$t[1]\t$depth\t$t[7]\n";
			if (rand() < 0.1) {
				print {$fh3} "$depth\t$t[7]\t$t[1]\n";
			}
		}
	}
	$var = sqrt($var / $n);
	close($fh1); close($fh2); close($fh3); close($fhin);
	open($fh, "| gnuplot") || die;
	my $label = sprintf("Depth (%.2f +/- %.2f)", $mean, $var);
	print $fh qq(
set t po eps co so 24
set grid;
set ylab "Depth";
set xlab "Percent GC";
set out "$prefix.gc.eps";
plot "$prefix.gc" u 2:1 w d t "";
set t po eps co so lw 2 18
set xlab "Coordinate (Mbp)";
set ylab "Depth";
set y2ran [0:100];
set y2tic 10;
set xtics $xtics;
set yran [-20:$max_depth];
set size 5, 1;
set grid;
set out "$prefix.depth.eps";
set style line 3 lt 3 lw 3;
plot "$prefix.depth" u 1:2 t "$label" w d 1, "$prefix.outlier" u 1:2 t "Outlier" w p 4, "$prefix.gc" u 3:2 smo be axis x1y2 ls 3;
exit;
);
	close($fh);
}

sub abpair
{
	my %opts;
	getopts('m:M:c:', \%opts);
	die qq(Usage: maq_plot.pl abpair [-m 250] [-M 1000000] [-c chr] <prefix> <abpair-output>\n) if (@ARGV < 2);
	my $min_dist = (defined $opts{m})? $opts{m} : 250;
	my $max_dist = (defined $opts{M})? $opts{M} : 1000000;
	my $prefix = $ARGV[0];
	my $chr = $opts{c};
	my ($fh1, $fh2, $fh, $fhin);
	open($fhin, "sort +3 -nr $ARGV[1] |") || die;
	open($fh1, ">$prefix.badori") || die;
	open($fh2, ">$prefix.baddist") || die;
	my (%hash1, %hash2);
	while (<$fhin>) {
		my @t = split;
		shift(@t);
		$chr = $t[0] if (!$chr);
		my $dist = $t[4] - $t[3] + 1;
		if ($t[0] eq $chr && $dist >= $min_dist && $dist <= $max_dist) {
			next if ($hash1{$t[3]} || $hash2{$t[4]});
			$hash1{$t[3]} = $hash2{$t[4]} = 1;
			$fh = ($t[1] == 2)? $fh2 : $fh1;
			my $q = $t[2] + rand() - 0.5;
			$q = 0 if ($q < 0);
			print $fh $t[3]/2000000.0 + $t[4]/2000000.0, "\t$q\t$dist\n";
		}
	}
	close($fhin); close($fh1); close($fh2);
	# plot
	open($fh, "| gnuplot") || die("Fail to execute gnuplot");
	print $fh qq(
set t po eps co so;
set xlab "$chr coordinate (Mbp)";
set ylab "MapQ of abnormal pairs";
set out "$prefix.abpair.eps";
set size 5, 1;
set xtics 2;
set grid;
plot "$prefix.badori" t "bad orientation" w p 2, "$prefix.baddist" t "bad distance" w p 4;
exit;
);
	close($fh);
}

sub mixed
{
	my %opts = (m=>300, M=>1000000, D=>100);
	getopts('c:D:m:M:', \%opts);
	die qq(
Usage:   maq_plot.pl mixed [options] <prefix> <cns2win-output> <abpair-output>

Options: -c STR        target chromosome [the first one]
         -D INT        maximum depth [100]
         -m INT        minumum pair distance [300]
         -M INT        maximum pair distance [1000000]
\n) if (@ARGV < 3); 
	my $chr_cmd = $opts{c}? "-c $opts{c}" : '';
	system(qq{$0 depth $chr_cmd -D $opts{D} $ARGV[0] $ARGV[1]}) && die;
	system(qq{$0 abpair $chr_cmd -m $opts{m} -M $opts{M} $ARGV[0] $ARGV[2]}) && die;
	my $prefix = $ARGV[0];
	my $fh;
	open($fh, "| gnuplot") || die;
	print $fh qq(
set t po eps co so;
set xlab "Coordinate (Mbp)";
set ylab "MapQ of abnormal pairs";
set out "$prefix.eps";
set size 5, 1;
set y2ran [-20:$opts{D}];
set y2lab "Depth";
set y2tics 10;
set xtics 2;
set grid;
plot "$prefix.badori" t "bad orientation" w p 2, "$prefix.outlier" t "Outlier" axis x1y2 3, "$prefix.baddist" t "bad distance" w p 4;
exit;
);
	close($fh);
}

sub usage
{
	die qq(
Usage:   maq_plot.pl <command> <arguments>

Command: abpair     plot abnormal pairs along the chromosome
         depth      depth plot
         mixed      mixed plot
\n);
}