annotate micomplete.pl @ 0:23e768cddc7d draft

Initial commit of micomplete

1 #!/usr/bin/perl -w
2 
3 =head1 SYNOPSIS
4 
5 micomplete.pl - Given a series of HMM profiles and a query proteome, and optionally, weights for the profiles, returns the completeness of a genome.
6 
7 =head1 USAGE
8 
9 micomplete.pl -h|--hmm_profiles -p|--proteome [-w|--weights] [-e|--evalue] [-s|--save-result-tab] [-r|--save-hmm-results] [-t|--threads]
10 micomplete.pl --help
11 
12 =head1 INPUT
13 
14 =head2 -h|--hmm
15 
16 A single file containing the HMM profiles against which the proteome will be run.
17 
18 =head2 -p|--proteome
19 
20 A fasta file containing all protein sequences in a genome.
21 
22 =head2 [-c|--cut-off]
23 
24 The evalue cut-off to accept the presence of a HMM. 1e-10 by default.
25 
26 =head2 [-w|--weights]
27 
28 Optionally, the weights to apply to the HMMs, for the completeness report. A tab file containing (i) a column with names corresponding to the names in the hmm profiles and (ii) a column with with the weighths. 
29 
30 =head2 [-e|--evalue]
31 
32 Shows evalue to the best hit, instead of presence/absence.
33 
34 =head2 [-s|--save-result-tab]
35 
36 Save the result per marker
37 
38 =head2 [-r|--save-hmm-tab]
39 
40 Save hmmsearch tabular result for this file.
41 
42 =head2 [-a|--save-hmm-alignments]
43 
44 Save hmmsearch tabular result for this file.
45 
46 =head2 [-t|--threads]
47 
48 Number of threads to use in the hmmsearch. 15 by default.
49 
50 =head2 --help
51 
52 Displays a help message.
53 
54 =head1 OUTPUT
55 
56 A completeness report expressed in fraction, and a list of all hmms identified as positive
57 
58 =head1 AUTHOR
59 
60 Lionel Guy ([e-mail])
61 
62 =head1 DATE
63 
64 Mon Oct 14 14:09:46 CEST 2013
65 
66 =cut
67 
68 # libraries
69 use strict;
70 use Pod::Usage;
71 use Getopt::Long;
72 # For safe use of tempfiles
73 use File::Temp qw/ tempfile tempdir /; 
74 use File::Copy qw/ cp /;
75 
76 # Options
77 my ($hmm_file, $fasta_in, $weights_file, $display_eval);
78 my $cutoff = 1e-10;
79 my ($save_tab_file, $save_aln_file, $save_res_tab);
80 my $nthreads = 15;
81 my $help;
82 GetOptions(
83     'h|hmm=s'      => \$hmm_file,
84     'p|proteome=s' => \$fasta_in,
85     'c|cut-off=s'  => \$cutoff,
86     'w|weights=s'  => \$weights_file,
87     'e|evalue'     => \$display_eval,
88     'r|save-hmm-tab=s'        => \$save_tab_file,
89     'a|save-hmm-alignments=s' => \$save_aln_file,
90     's|save-result-tab=s'     => \$save_res_tab,
91     't|threads=s'  => \$nthreads,
92     'help'         => \$help,
93 );
94 
95 # Tempfile if not saving results
96 my $template = 'tempXXXXXXXX';
97 my ($fh_tab, $tab_file) = tempfile($template, DIR => '.', UNLINK => 1);
98 my ($fh_aln, $aln_file) = tempfile($template, DIR => '.', UNLINK => 1);
99 
00 # Check options
01 pod2usage(-verbose => 1) if $help;
02 pod2usage(-message => "No hmm profile file (-h|--hmm) defined\n", 
03 	  -exitstatus => 2 ) unless ($hmm_file);
04 pod2usage(-message => "No proteome file (-p|--proteome) defined\n", 
05 	  -exitstatus => 2 ) unless ($fasta_in);
06 
07 # Parse eventual weights
08 my %hmms;
09 if ($weights_file){
10     open WEIGHTS, '<', $weights_file or die "Cannot open $weights_file\n";
11     while (<WEIGHTS>){
12 	next if /^#/;
13 	my ($name, $weight) = split;
14 	$hmms{$name} = $weight;
15     }
16 }
17 
18 # Parse HMM, lists names
19 open HMMS, '<', $hmm_file or die "Cannot open $hmm_file\n";
20 while (<HMMS>){
21     next unless /^NAME/;
22     my ($bla, $name) = split;
23     if ($weights_file) {
24 	die "HMM $name doesn't have weight defined" unless $hmms{$name};
25     } else {
26 	$hmms{$name}++;
27     }
28 }
29 
30 # Get sum of weights
31 my $total;
32 foreach (values %hmms) {
33     $total += $_;
34 }
35 #print STDERR "Sum of weights: $total\n";
36 
37 # Run hmmsearch
38 my $cmd = "hmmsearch --tblout $tab_file --cpu $nthreads ";
39 $cmd .= "--noali " unless $save_aln_file;
40 $cmd .= "-E $cutoff $hmm_file $fasta_in > $aln_file";
41 system($cmd) == 0 or die ("Command $cmd failed: $?");
42 
43 # Read HMMER output. Assumes the first hit is always the best
44 open HMMER, '<', $tab_file or die $!;
45 my %seen;
46 my $sum;
47 while (<HMMER>){
48     next if /^#/;
49     my ($gid, $dash, $hname, $hid, $evalue) = split;
50     #my ($bla, $gi, $ref, $pid) = split(/\|/, $gid);
51     die ("No gid at $_\n") unless $gid;
52     die ("HMM not seen in the hmm file") unless $hmms{$hname};
53     # Skip if not first hit, or if evalue under cut-off
54     unless ($seen{$hname} || $evalue > $cutoff){
55 	$sum += $hmms{$hname};
56 	$seen{$hname} = $evalue;
57     }
58 }
59 
60 # Save results if required
61 cp($tab_file, $save_tab_file) if ($save_tab_file);
62 cp($aln_file, $save_aln_file) if ($save_aln_file);
63 
64 # Print results
65 printf "Completeness: %.4f", $sum/$total;
66 print "\nN markers found: ", scalar(keys %seen), " out of ", 
67    scalar(keys %hmms), "\n";
68 if ($save_res_tab){
69     open TAB, '>', $save_res_tab or die "Could not open $save_res_tab: $?\n";
70     print TAB "#marker\t"; 
71     if ($display_eval) {
72 	print TAB "evalue\n";
73     } else {
74 	print TAB "present\n";
75     }
76     foreach my $hmm (sort keys %hmms){
77 	print TAB "$hmm\t";
78 	if ($seen{$hmm}){
79 	    if ($display_eval){
80 		print TAB "$seen{$hmm}\n";
81 	    } else {
82 		print TAB "1\n";
83 	    }
84 	} else {
85 	    if ($display_eval){
86 		print TAB "NA\n";
87 	    } else {
88 		print TAB "0\n";
89 	    }
90 	}
91     }
92 }
93