Mercurial > repos > big-tiandm > mirplant2

comparison collapseReads2Tags.pl @ 26:c90e7d84d521 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Thu, 31 Jul 2014 03:03:37 -0400
parents ea2fdf667620
children
comparison
equal deleted inserted replaced
25:10df3c84d54a 26:c90e7d84d521
1 #!/usr/bin/perl -w
2 #Filename:
3 #Author: Tian Dongmei
4 #Email: tiandm@big.ac.cn
5 #Date: 2014-3-20
6 #Modified:
7 #Description: fastq file form reads cluster(the same sequence in the same cluster)
8 my $version=1.00;
9
10 use strict;
11 use Getopt::Long;
12
13 my %opts;
14 GetOptions(\%opts,"i:s@","format=s","mark:s","qual:s","qv:i","o=s","h");
15 if (!(defined $opts{o} and defined $opts{'format'}) || defined $opts{h}) { #necessary arguments
16 &usage;
17 }
18 my @filein=@{$opts{i}} if(defined $opts{i});
19 my $name=defined $opts{'mark'} ? $opts{'mark'} : "seq";
20 my $fileout=$opts{'o'};
21 my $pq=defined $opts{'qv'} ? $opts{'qv'} : 33;
22 my %hash;##分块存放原始序列
23
24 my $format=$opts{'format'};
25 if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") {
26 die "Parameter -format is error!\n";
27 }
28
29 my ($qualT,$qualV);
30 if (defined $opts{'qual'} && ($format eq "fastq" || $format eq "fq")) { #quality filter
31 my @temp=split /:/,$opts{'qual'};
32 $qualT=$temp[0];
33 $qualV=$temp[1];
34
35 for (my $i=0;$i<@filein;$i++) {
36 open IN,"<$filein[$i]";
37 while (my $aline=<IN>) {
38 my $seq=<IN>;
39 my $n=<IN>;
40 my $qv=<IN>;
41 my $tag=&qvcheck($qv,$qualT,$qualV);
42 next if(!$tag);
43 my $str=substr($seq,0,6);
44 $hash{$str}[$i].=$seq;
45 }
46 close IN;
47 }
48 }
49 elsif($format eq "fastq" || $format eq "fq"){ ### do not filter low quality reads
50 for (my $i=0;$i<@filein;$i++) {
51 open IN,"<$filein[$i]";
52 while (my $aline=<IN>) {
53 my $seq=<IN>;
54 my $n=<IN>;
55 my $qv=<IN>;
56 my $str=substr($seq,0,6);
57 $hash{$str}[$i].=$seq;
58 }
59 close IN;
60 }
61
62 }
63 elsif($format eq "fasta" || $format eq "fa"){
64 for (my $i=0;$i<@filein;$i++) {
65 open IN,"<$filein[$i]";
66 while (my $aline=<IN>) {
67 my $seq=<IN>;
68 my $str=substr($seq,0,6);
69 $hash{$str}[$i].=$seq;
70 }
71 close IN;
72 }
73 }
74
75 open OUT,">$fileout"; #output file
76 my $count=0;
77 foreach my $key (keys %hash) {
78 my %cluster;
79 for (my $i=0;$i<@filein;$i++) {
80 next if(!(defined $hash{$key}[$i]));
81 my @tmp=split/\n/,$hash{$key}[$i];
82 foreach (@tmp) {
83 $cluster{$_}[$i]++;
84 }
85 }
86
87 foreach my $seq (keys %cluster) {
88 my $exp=""; my $ee=0;
89 for (my $i=0;$i<@filein;$i++) {
90 if (defined $cluster{$seq}[$i]) {
91 $exp.="_$cluster{$seq}[$i]";
92 $ee+=$cluster{$seq}[$i];
93 }else{
94 $exp.="_0";
95 }
96 }
97 $count+=$ee;
98 $exp=~s/^_//;
99 print OUT ">$name","_$count:$exp","_x$ee\n$seq\n";
100 }
101 }
102 close OUT;
103
104
105 sub qvcheck{
106 my ($str,$t,$v)=@_;
107 my $qv=0;
108 if($t eq "mean"){
109 $qv=&getMeanQuality($str);
110 }
111 elsif($t eq "min"){
112 $qv=&getMinQuality($str);
113 }
114 if ($qv<$v) {
115 return 0;
116 }
117 return 1;
118 }
119
120 sub getMeanQuality(){
121 chomp $_[0];
122 my @bases = split(//,$_[0]);
123 my $sum = 0;
124 for(my $i = 0; $i <= $#bases; $i++){
125 my $num = ord($bases[$i]) - $pq;
126 $sum += $num;
127 }
128
129 return $sum/($#bases+1);
130
131 }
132
133 ###
134 ### This function gives back the Q-value of the worst base
135 sub getMinQuality(){
136 chomp $_[0];
137 my @bases = split(//,$_[0]);
138 my $worst = 1000;
139 for(my $i = 0; $i <= $#bases; $i++){
140 # printf ("base: $bases[$i] --> %d\n",ord($bases[$i]));
141 my $num = ord($bases[$i]) - $pq;
142 if($num < $worst){
143 $worst = $num;
144 }
145 }
146 return $worst;
147 }
148
149 sub usage{
150 print <<"USAGE";
151 Version $version
152 Usage:
153 $0 -i -format -mark -qual -qv -o
154 options:
155 -i input file#fastq file ##can be multiple -i file1 -i file2 ...
156 -mark string#quary name,default is "seq"
157 -o output file
158 -format string # fastq|fasta|fq|fa
159
160 -qual #reads filter
161 eg:(min:value/mean:value)
162 This parameter just for solexa reads.
163 If the input files are solid and needs filter,please do filter first .
164
165 -qv integer #Phred quality64/33,default 33
166 -h help
167 USAGE
168 exit(1);
169 }
170