1
|
1 #!/usr/bin/perl
|
|
2
|
|
3 use strict;
|
|
4 my %read1;
|
|
5 my %read2;
|
|
6 my %qual1;
|
|
7 my %qual2;
|
|
8
|
|
9 while (my $line=<STDIN>) {
|
|
10 chomp();
|
|
11 my ($name, $flag, $chr1, $pos1, $mq, $cigar, $chr2, $pos2, $dist, $seq, $quals, $stuff)=split(/\s+/, $line);
|
|
12 if($flag & 16) { #reverse-complement it
|
|
13 $seq=reverse($seq);
|
|
14 $seq=~tr/ACGT/TGCA/;
|
|
15 $quals=reverse($quals);
|
|
16 }
|
|
17 if($flag & 64) { #first in pair
|
|
18 $read1{$name}=$seq;
|
|
19 $qual1{$name}=$quals;
|
|
20 }
|
|
21 else {
|
|
22 $read2{$name}=$seq;
|
|
23 $qual2{$name}=$quals;
|
|
24 }
|
|
25 }
|
|
26
|
|
27 open(FQ1, ">$ARGV[0]") or die "Can't open $ARGV[0]";
|
|
28 open(FQ2, ">$ARGV[1]") or die "Can't open $ARGV[1]";
|
|
29 open(FQ3, ">$ARGV[2]") or die "Can't open $ARGV[2]";
|
|
30
|
|
31 foreach my $key(keys %read1) {
|
|
32
|
|
33 if(exists $read2{$key}) {
|
|
34 print FQ1 "\@$key\n$read1{$key}\n+\n$qual1{$key}\n";
|
|
35 print FQ2 "\@$key\n$read2{$key}\n+\n$qual2{$key}\n";
|
|
36 }
|
|
37 else {
|
|
38 print FQ3 "\@$key\n$read1{$key}\n+\n$qual1{$key}\n";
|
|
39 }
|
|
40 }
|
|
41 foreach my $key(keys %read2) {
|
|
42
|
|
43 if(!(exists $read1{$key})) {
|
|
44 print FQ3 "\@$key\n$read2{$key}\n+\n$qual2{$key}\n";
|
|
45 }
|
|
46 }
|
|
47 close(FQ1);
|
|
48 close(FQ2);
|
|
49 close(FQ3);
|
|
50
|
|
51
|
|
52
|
|
53 |