Mercurial > repos > bioitcore > splicetrap

comparison bin/rmap2eland.pl @ 1:adc0f7765d85 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload
author bioitcore
date Thu, 07 Sep 2017 15:06:58 -0400
parents
children
comparison
equal deleted inserted replaced
0:d4ca551ca300 1:adc0f7765d85
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);