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 89 90 91 92 93 94 95 96 97 98 99 100 101 102
|
#!/usr/bin/env perl
use strict;
use warnings;
use lib ("/usr/lib/trinityrnaseq/PerlLib");
use SAM_reader;
use SAM_entry;
my $usage = "\n\nusage: name_sorted_paired_reads.sam\n\n";
my $sam_file = $ARGV[0] or die $usage;
main: {
my $prev_read_name = "";
my $prev_scaff_name = "";
my @reads;
my $sam_reader = new SAM_reader($sam_file);
while ($sam_reader->has_next()) {
my $read = $sam_reader->get_next();
my $scaff_name = $read->get_scaffold_name();
my $core_read_name = $read->get_core_read_name();
if ($scaff_name ne $prev_scaff_name || $core_read_name ne $prev_read_name) {
if (@reads) {
&process_pairs(@reads);
@reads = ();
}
}
push (@reads, $read);
$prev_read_name = $core_read_name;
$prev_scaff_name = $scaff_name;
}
&process_pairs(@reads);
exit(0);
}
####
sub process_pairs {
my (@reads) = @_;
my @left_reads;
my @right_reads;
foreach my $read (@reads) {
if ($read->is_first_in_pair()) {
push (@left_reads, $read);
}
elsif ($read->is_second_in_pair()) {
push (@right_reads, $read);
}
}
unless (@left_reads && @right_reads) {
die "Error, dont have pairs!";
}
my $left_read = shift @left_reads;
my $aligned_pos = $left_read->get_aligned_position();
my $right_read = undef;
while ((! defined($right_read)) && @right_reads) {
my $read = shift @right_reads;
if ($read->get_mate_scaffold_position() == $aligned_pos) {
$right_read = $read;
last;
}
}
unless ($left_read && $right_read) {
die "Error, couldn't match pairs.";
}
print $left_read->toString() . "\n";
print $right_read->toString() . "\n";
return;
}
|