/usr/bin/hits_to_fastq is in seer 1.1.2-3.
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 | #!/usr/bin/perl -w
use strict;
use warnings;
use Getopt::Long;
my $usage_message = <<USAGE;
Usage: hits_to_fastq -k <significant_kmers.txt> -b <bonferroni_correction>
Creates a fastq for mapping from significant kmers
Options
Required
-k, --kmers Significant kmers, as output by seer
-b, --bonf Bonferroni correction to use (try 10e-8)
Optional
-h, --help Displays this message
USAGE
#****************************************************************************************#
#* Main *#
#****************************************************************************************#
#* gets input parameters
my ($kmer_file, $bonf, $help);
GetOptions ("kmers|k=s" => \$kmer_file,
"bonf|b=s" => \$bonf,
"help|h" => \$help
) or die($usage_message);
if (defined($help) || !defined($kmer_file) || !defined($bonf))
{
print $usage_message;
}
else
{
open(KMERS, $kmer_file) || die("Could not open kmer file: $kmer_file\n");
my $header = <KMERS>;
my $i = 1;
while (my $line_in = <KMERS>)
{
chomp $line_in;
my ($kmer, $maf, $p_unadj, $p_adj, @junk) = split("\t", $line_in);
my $corrected_p = $p_adj / $bonf;
if ($corrected_p > 1)
{
$corrected_p = 1;
}
my $ascii_val;
if ($corrected_p == 0)
{
$ascii_val = 126;
}
else
{
$ascii_val = int(-10*log($corrected_p)/log(10)) + 33;
if ($ascii_val < 33)
{
$ascii_val = 33;
}
elsif ($ascii_val > 126)
{
$ascii_val = 126;
}
}
my $qual = chr($ascii_val) x length($kmer);
print join("\n", "\@$i", $kmer, "+", $qual) . "\n";
$i++;
}
close KMERS
}
exit(0);
|