Mercurial > repos > bigrna > gpsrna

comparison filterReadsByCount.pl @ 0:87fe81de0931 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author bigrna
date Sun, 04 Jan 2015 02:47:25 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:87fe81de0931
1 #!/usr/bin/perl -w
2 #Filename:
3 #Author: Tian Dongmei
4 #Email: tiandm@big.ac.cn
5 #Date: 2010-01
6 #Modified:
7 #Description:
8 my $version=1.00;
9
10 use strict;
11 use Getopt::Long;
12 use File::Basename;
13
14 my %opts;
15 GetOptions(\%opts,"i=s","o=s","mark:s","h");
16 if (!(defined $opts{i} and defined $opts{o}) || defined $opts{h}) { #necessary arguments
17 &usage;
18 }
19
20 my $mark=defined $opts{'mark'} ? $opts{'mark'} : "Sample";
21 my @mark=split /\#/,$mark;
22
23 open OUT,">$opts{o}";
24 open IN,"<$opts{i}";
25 my %hash;my %reads;
26 while (my $aline=<IN>) {
27 chomp $aline;
28 my $seq=<IN>;
29 chomp $seq;
30 if($aline=~/:([\d|_]+)_x(\d+)$/){
31 if ($2>3) {
32 my @ss=split/_/,$1;
33 for (my $i=0;$i<@ss;$i++) {
34 $hash{length($seq)}[$i]++ if($ss[$i]>0);
35 $hash{length($seq)}[$i] +=0 if($ss[$i]==0);
36 $reads{length($seq)}[$i]+=$ss[$i];
37 }
38 print OUT "$aline\n$seq\n";
39 }
40 }
41 }
42 close IN;
43 close OUT;
44
45 my $dir=dirname($opts{'o'});
46 chdir $dir;
47 my $lengthfile=$dir."/reads_length_distribution_after_count_filter.txt";
48 open OUT, ">$lengthfile";
49 open R,">$dir/length_distribution_after_count_filter.R";
50
51 print OUT "Tags length\t@mark\n";
52
53 my $samNo=@mark;
54 my $avalue="";
55 my @length=sort{$a<=>$b} keys %hash;
56 foreach (@length) {
57 print OUT $_,"\t@{$hash{$_}}\n";
58 my $vv=join ", ",@{$hash{$_}};
59 $avalue .="$vv,";
60 }
61 $avalue =~s/,$//;
62 my $lengths=join ",",@length;
63 my $marks=join "\",\"",@mark;
64
65 print R "a<-c($avalue)
66 b<-matrix(a,ncol=$samNo,byrow=T)
67 cl<-colors()
68 names=c($lengths)
69 legends=c(\"$marks\")
70 png(\"Tags_length_after_count_filter.png\",width=800,height=600)
71 barplot(t(b),beside=TRUE,col=cl[1:$samNo],main=\"Tags Length Distribution After Count Filter\",names.arg=names,ylim=c(0,max(a)),legend.text=legends,args.legend=\"topleft\")
72 abline(h=0)
73 dev.off()
74
75 ";
76 $avalue="";
77 print OUT "\nReads length\t@mark\n";
78 foreach (@length) {
79 print OUT $_,"\t@{$reads{$_}}\n";
80 my $vv=join ", ", @{$reads{$_}};
81 $avalue .= "$vv,";
82 }
83 $avalue =~s/,$//;
84
85 print R "a<-c($avalue)\n
86 b<-matrix(a,ncol=$samNo,byrow=T)
87
88 png(\"Reads_length_after_count_filter.png\",width=800,height=600)
89 barplot(t(b),beside=TRUE,col=cl[1:$samNo],main=\"Reads Length Distribution After Count Filter\",names.arg=names,ylim=c(0,max(a)),legend.text=legends,args.legend=\"topleft\")
90 abline(h=0)
91 dev.off()
92
93 ";
94 close OUT;
95 close R;
96
97 system ("R CMD BATCH $dir/length_distribution_after_count_filter.R");
98
99 #system ("rm $dir/length_distribution.R");
100 #system ("rm $dir/length_distribution.Rout");
101 #system ("rm $dir/.RData");
102 sub usage{
103 print <<"USAGE";
104 Version $version
105 Usage:
106 $0 -i -o -min -max -mark
107 options:
108
109 -i input file
110 -o output file
111 -mark string #sample name eg: samA#samB#samC
112 -h help
113 USAGE
114 exit(1);
115 }
116