perl-for-reysenbach-lab/break_scaffolds_into_orfs.pl at master · duffee/perl-for-reysenbach-lab · GitHub

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
#!/usr/bin/perl -w
##########################
### Jennifer Meneghin  ###
### August 4, 2009     ###
##########################
#-----------------------------------------------------------------------------------------------------------------------------------------
#Deal with passed parameters
#-----------------------------------------------------------------------------------------------------------------------------------------
#If no arguments are passed, show usage message and exit program.
if ($#ARGV == -1) {
    &usage;
}
#get the names of the blast file (first argument) and query file (second argument)
$fasta_file = $ARGV[0];
unless ( open(IN, "$fasta_file") ) {
    print "Got a bad fasta fasta file: $fasta_file\n";
    &usage;
}
unless ( open(OUT, ">open_reading_frames.fasta") ) {
    print "Couldn't create new fasta file: open_reading_frames.fasta\n";
    &usage;
}
print "Parameters:\nfasta file = $fasta_file\nopen reading frames = open_reading_frames.fasta\n\n";
#-----------------------------------------------------------------------------------------------------------------------------------------
#The main event
#-----------------------------------------------------------------------------------------------------------------------------------------
%scaffolds = ();
$sequence = "";
$header = "";
while (<IN>) {
    chomp;
    if (/^>/) {
	if (length($sequence) > 0) {
	    $scaffolds{$header} = $sequence;
	}
	$header = $_;
	$sequence = "";
    }
    else {
	$sequence = $sequence . $_;
    }
}
if (length($sequence) > 0) {
    $scaffolds{$header} = $sequence;
}
$last_index = 0;
$bad_counter = 0;
$good_counter = 0;
foreach $i (sort keys %scaffolds) {
    $seq = $scaffolds{$i};
    @bases = split(//, $seq);
    $counter = 0;
    foreach $j (0..$#bases) {
	if ($j+2 <= $#bases) {
	    $codon = $bases[$j] . $bases[$j+1] . $bases[$j+2];
	    if ($codon eq "TAG" || $codon eq "TGA" || $codon eq "TAA") {
		if ($j-$last_index > 99) {
		    print ".";
		    $counter++;
		    $scaffold_id = $i;
		    $scaffold_id =~ s/>(.+?)\s.+$/$1/g;
		    print OUT ">${scaffold_id}_${counter}\n";
		    foreach $k ($last_index..$j-1) {
			print OUT "$bases[$k]";
		    }
		    print OUT "\n";
		    $last_index = $j+3;
		    $good_counter++;
		}
		else {
		    $bad_counter++;
		}
	    }
	}
    }
}
print "\nGOOD = $good_counter\n";
print "BAD = $bad_counter\n";

sub usage {
    print "\nUSAGE: break_scaffolds_into_orfs.pl <scaffold fasta file>\n\n";
    print "This program breaks a fasta file of scaffolds into open reading frames*,\n";
    print "and returns a fasta file of all ORFs found in each scaffold with a minimum length of 99bp.\n\n";
    print "* An open reading frame is defined as the nucleotides between to two stop codons.\n\n";
    print "Jennifer Meneghin\n";
    print "08/06/2009\n\n";
    exit;
}