File: run_jellyfish.pl

package info (click to toggle)
trinityrnaseq 2.2.0%2Bdfsg-2
  • links: PTS, VCS
  • area: main
  • in suites: stretch
  • size: 212,452 kB
  • ctags: 5,067
  • sloc: perl: 45,552; cpp: 19,678; java: 11,865; sh: 1,485; makefile: 613; ansic: 427; python: 313; xml: 83
file content (71 lines) | stat: -rwxr-xr-x 1,406 bytes parent folder | download | duplicates (3)
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
#!/usr/bin/env perl

use strict;
use warnings;
use Cwd;
use FindBin;

my $usage = "\n\nusage: $0 reads.fa hash_size\n\n";

my $reads_file = $ARGV[0] or die $usage;
my $hash_size = $ARGV[1] or die $usage;


my $JELLYFISH_DIR = "/usr";
my $CPU = 4;
my $min_kmer_cov = 1;

unless ($reads_file =~ /^\//) {
    $reads_file = cwd() . "/$reads_file";
}


my $workdir = "H_" . ($hash_size/1e9) . "G";
mkdir($workdir) or die "Error, cannot mkdir $workdir";
chdir ($workdir) or die "Error, cannot cd to $workdir";


my $jelly_kmer_fa_file = "jellyfish.kmers.fa";
            
#    my $jelly_hash_size = int( ($max_memory - $read_file_size)/7); # decided upon by Rick Westerman
             
my $cmd = "$JELLYFISH_DIR/bin/jellyfish count -t $CPU -m 25 -s $hash_size ";

#        $cmd .= " --both-strands ";
            
$cmd .= " $reads_file";
            
&process_cmd($cmd);

my @kmer_db_files;

foreach my $file (<mer_counts_*>) {
    my $cmd = "$JELLYFISH_DIR/bin/jellyfish dump -L $min_kmer_cov $file >> $file.kmer_fa";
    &process_cmd($cmd);
    
    $cmd = "cat $file.kmer_fa >> $jelly_kmer_fa_file";
    &process_cmd($cmd);
    
}
            
$cmd = "grep '>' $jelly_kmer_fa_file | wc -l | tee kmer_count.txt";
&process_cmd($cmd);


exit(0);


####
sub process_cmd {
    my ($cmd) = @_;

    my $ret = system($cmd);

    if ($ret) {
        die "Error, cmd: $cmd died with ret $ret";
    }

    return;
}