annotate collapseReads2Tags.pl @ 31:7321a6f82492 draft

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