|
7
|
1 use strict;
|
|
|
2
|
|
|
3 my $rmapfilename=$ARGV[0];
|
|
|
4 my $readsfilename=$ARGV[1];
|
|
|
5 my $elandfilename=$ARGV[2];
|
|
|
6
|
|
|
7 my $detectformat=`head -c 1 $readsfilename`;
|
|
|
8
|
|
|
9 #system("grep \"$detectformat\" $readsfilename |sort >$readsfilename.sort");
|
|
|
10 system("awk 'NR%2==1' $readsfilename |sort >$readsfilename.sort");
|
|
|
11 system("sort -k4,4 $rmapfilename >$rmapfilename.sort");
|
|
|
12
|
|
|
13
|
|
|
14 open(readsfile, $readsfilename.".sort");
|
|
|
15
|
|
|
16
|
|
|
17
|
|
|
18 #$looplinenumbers=2 if ($detectformat eq ">");
|
|
|
19 open(rmapfile, $rmapfilename.".sort");
|
|
|
20 open(elandfile, ">".$elandfilename);
|
|
|
21
|
|
|
22 while(my $rmapline=<rmapfile>)
|
|
|
23 {
|
|
|
24 chomp($rmapline);
|
|
|
25 my ($mapped_id, $start, $end, $rmapreadname, $mismatch, $strand)=split("\t",$rmapline);
|
|
|
26 while(my $readline=<readsfile>)
|
|
|
27 {
|
|
|
28 if($readline=~/^$detectformat/)
|
|
|
29 {
|
|
|
30 chomp($readline);
|
|
|
31 my $readname=substr($readline, 1, length($readline)-1);
|
|
|
32
|
|
|
33
|
|
|
34 if($readname ne $rmapreadname)
|
|
|
35 {
|
|
|
36 print elandfile $readname,"\tNA\tNM\n";
|
|
|
37 next;
|
|
|
38 }
|
|
|
39 else
|
|
|
40 {
|
|
|
41 my @mapped_ids=();
|
|
|
42 my @mapped_pos=();
|
|
|
43 my @mapped_strand=();
|
|
|
44 push(@mapped_ids, $mapped_id);
|
|
|
45 push(@mapped_pos,$start);
|
|
|
46 push(@mapped_strand,$strand);
|
|
|
47 while(1)
|
|
|
48 {
|
|
|
49 $rmapline=<rmapfile>;
|
|
|
50 chomp($rmapline);
|
|
|
51 ($mapped_id, $start, $end, $rmapreadname, $mismatch, $strand)=split("\t",$rmapline);
|
|
|
52 if( $rmapreadname eq $readname )
|
|
|
53 {
|
|
|
54 push(@mapped_ids, $mapped_id);
|
|
|
55 push(@mapped_pos,$start);
|
|
|
56 push(@mapped_strand,$strand);
|
|
|
57 }
|
|
|
58 else
|
|
|
59 {
|
|
|
60 seek(rmapfile, -1*length($rmapline)-1,1);
|
|
|
61 print elandfile $readname,"\t";
|
|
|
62 print elandfile "NA\t";
|
|
|
63 print elandfile scalar(@mapped_ids),":0:0\t";
|
|
|
64 for(my $i=0;$i<@mapped_ids;$i++)
|
|
|
65 {
|
|
|
66 print elandfile "/",$mapped_ids[$i];
|
|
|
67 print elandfile ":",$mapped_pos[$i]+1;
|
|
|
68 if($mapped_strand[$i] eq "+")
|
|
|
69 {
|
|
|
70 print elandfile "F0,";
|
|
|
71 }
|
|
|
72 else
|
|
|
73 {
|
|
|
74 print elandfile "R0,";
|
|
|
75 }
|
|
|
76
|
|
|
77 }
|
|
|
78 print elandfile "\n";
|
|
|
79 last;
|
|
|
80 }
|
|
|
81 }
|
|
|
82 last;
|
|
|
83
|
|
|
84 }
|
|
|
85 }
|
|
|
86 }
|
|
|
87 }
|
|
|
88
|
|
|
89 while(my $readline=<readsfile>)
|
|
|
90 {
|
|
|
91 if($readline=~/^$detectformat/)
|
|
|
92 {
|
|
|
93 chomp($readline);
|
|
|
94 my $readname=substr($readline, 1, length($readline)-1);
|
|
|
95 print elandfile $readname,"\tNA\tNM\n";
|
|
|
96 }
|
|
|
97 }
|
|
|
98
|
|
|
99 close(elandfile);
|
|
|
100 close(rmapfile);
|
|
|
101
|
|
|
102
|
|
|
103 close(readsfile);
|
