Mercurial > repos > lionelguy > spades

comparison tools/spades_3_1_1/tools/spades_3_5_0/filter_spades_output.pl @ 12:85c6121d92a5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author takadonet
date Tue, 16 Dec 2014 12:57:49 -0500
parents
children
comparison
equal deleted inserted replaced
11:1027f26f52d5 12:85c6121d92a5
1 #!/usr/bin/env perl -w
2
3 =head1 SYNOPSIS
4
5 filter_spades_output.pl - Filters contigs or scaffolds based on contig length and coverage.
6
7 =head1 USAGE
8
9 filter_spades_output.pl [-c|--coverage-cutoff] [-l|--length-cutoff] [-o|--filtered-out out.fasta] -t|--tab stats.tab seqs.fasta
10
11 =head1 INPUT
12
13 =head2 [-c|--coverage-cutoff]
14
15 Mininum coverage. Contigs with lower coverage will be discarded. Default 10.
16
17 =head2 [-l|--length-cutoff]
18
19 Mininum coverage. Smaller ontigs will be discarded. Default 500.
20
21 =head2 -t|--tab stats.tab
22
23 A tabular file, with three columns: contig name, length, and coverage:
24
25 NODE_1 31438 24.5116
26 NODE_2 31354 2316.96
27 NODE_3 26948 82.3294
28
29 Such a file is produced by spades.xml. Contigs should be in the same order as in the fasta file.
30
31 =head2 [-o|--filtered-out out.fasta]
32
33 If specified, filtered out sequences will be written to this file.
34
35 =head2 seqs.fasta
36
37 Sequences in fasta format. Start of IDs must match ids in the tabular file.
38
39 =head1 OUTPUT
40
41 A fasta file on stdout.
42
43 =head1 AUTHOR
44
45 Lionel Guy (lionel.guy@icm.uu.se)
46
47 =head1 DATE
48
49 Thu Aug 29 13:51:13 CEST 2013
50
51 =cut
52
53 # libraries
54 use strict;
55 use Getopt::Long;
56 use Bio::SeqIO;
57
58 my $coverage_co = 10;
59 my $length_co = 500;
60 my $out_filtered;
61 my $tab_file;
62
63 GetOptions(
64 'c|coverage-cutoff=s' => \$coverage_co,
65 'l|length-cutoff=s' => \$length_co,
66 'o|filtered-out=s' => \$out_filtered,
67 't|tab=s' => \$tab_file,
68 );
69 my $fasta_file = shift(@ARGV);
70 die ("No tab file specified") unless ($tab_file);
71 die ("No fasta file specified") unless ($fasta_file);
72
73 ## Read tab file, discard rows with comments
74 open TAB, '<', $tab_file or die "$?";
75 my @stats;
76 while (<TAB>){
77 chomp;
78 push @stats, $_ unless (/^#/);
79 }
80
81 ## Read fasta
82 my $seq_in = Bio::SeqIO->new(-file => $fasta_file,
83 -format => 'fasta');
84 my $seq_out = Bio::SeqIO->new(-fh => \*STDOUT,
85 -format => 'fasta');
86 my $seq_out_filt = Bio::SeqIO->new(-file => ">$out_filtered",
87 -format => 'fasta') if ($out_filtered);
88 while (my $seq = $seq_in->next_seq){
89 my $stat = shift @stats;
90 die "Less rows in tab than sequences in seq file" unless $stat;
91 my ($id_tab, $length, $coverage) = split(/\t+/, $stat);
92 die "id, length or coverate not defined at $stat\n"
93 unless ($id_tab && $length && $coverage);
94 my $id_seq = $seq->id;
95 die "Unmatched ids $id_seq and $id_tab\n" unless ($id_seq =~ /^$id_tab/i);
96 if ($length >= $length_co && $coverage >= $coverage_co){
97 $seq_out->write_seq($seq);
98 } elsif ($out_filtered){
99 $seq_out_filt->write_seq($seq);
100 } else {
101 # do nothing
102 }
103 }
104 die "More rows in tab than sequences in seq file" if (scalar(@stats) > 0);
105 exit 0;
106