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
|
#!/usr/bin/perl
use strict;
$0 = rindex($0, "/") > -1 ? substr($0, rindex($0, "/")+1) : $0;
my (@lines, @filt_lines);
my ($line, $line_in, $type);
my $mode = ($ARGV[0] eq "1" ? "M1" : ($ARGV[0] eq "2" ? "M2" : die("$0: Invalid base genome argument (expected 1 or 2)")));
die("$0: LAGAN_DIR not defined. Stopped") unless defined $ENV{"LAGAN_DIR"};
while (<STDIN>) {
$line_in = $_;
if ($line_in =~ /^\=.*(DM|M1|M2)$/) {
$type = $1; $line .= $line_in;
$lines[$#lines+1] = $line if $type eq "DM" or $type eq $mode;
undef $line; undef $type;
} else {
$line .= $line_in;
}
}
foreach my $line (@lines) {
if ($mode eq "M2") {
$line =~ /(\>[^\s\n]+\s([\+\-])[^\n]+)\n(.+)\n(\>[^\s\n]+\s([\+\-])[^\n]+)\n(.+)\n(\=.+?)\n/s;
# $line =~ /(\>[^\s\n]+\s([\+\-])[^\n]+)\n([^\n]+)\n(\>[^\s\n]+\s([\+\-])[^\n]+)\n([^\n]+)\n(\=.+?\n)/s;
my ($head1, $strand1, $seq1, $head2, $strand2, $seq2, $foot) = ($1, $2, $3, $4, $5, $6, $7);
die if $strand1 ne $strand2;
if ($strand1 eq "-") {
$seq1 =~ s/\n//g;
$seq2 =~ s/\n//g;
$seq1 = reverse($seq1);
$seq2 = reverse($seq2);
$seq1 =~ s/(.{80})/$1\n/g;
$seq2 =~ s/(.{80})/$1\n/g;
}
$line = $head2."\n".$seq2."\n".$head1."\n".$seq1."\n".$foot."\n";
}
push @filt_lines, $line;
}
open(OUT, "> tmp.xmfa");
foreach my $line (@filt_lines) { print OUT $line; }
close OUT;
system($ENV{"LAGAN_DIR"}."/utils/Glue tmp.xmfa > glue.out 2> glue.err");
open(IN, "< glue.out");
my @glue_out = <IN>;
close IN;
open(IN, "< glue.err");
my @glue_err = <IN>;
close IN;
unlink("tmp.xmfa");
unlink("glue.out");
unlink("glue.err");
print STDOUT @glue_out;
print STDERR @glue_err;
|
