annotate preProcess.pl @ 53:f5a2e8308836 draft default tip

Uploaded
author big-tiandm
date Mon, 08 Dec 2014 01:51:16 -0500
parents c75593f79aa9
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
47
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
1 #!/usr/bin/perl -w
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
2 #Filename:
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
3 #Author: Tian Dongmei
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
4 #Email: tiandm@big.ac.cn
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
5 #Date: 2014-12-2
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
6 #Modified:
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
7 #Description: RNA-seq data pre-process
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
8 my $version=1.00;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
9
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
10 use strict;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
11 use Getopt::Long;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
12 use threads;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
13 #use threads::shared;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
14 use File::Path;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
15 use File::Basename;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
16 #use RNA;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
17 #use Term::ANSIColor;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
18
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
19 my %opts;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
20 GetOptions(\%opts,"i:s@","tag:s@","format=s","phred:i","gfa=s","rfam:s","idx:s","idx2:s","mis:i","v:i","a:s","M:i","t:i","min:i","max:i","o:s","path:s","h");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
21 if (!(defined $opts{i} and defined $opts{format} and defined $opts{gfa} ) || defined $opts{h}) { #necessary arguments
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
22 &usage;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
23 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
24
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
25 my $time=&Time();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
26 print "miPlant program start:\n The time is $time!\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
27 print "Command line:\n $0 @ARGV\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
28
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
29 my $format=$opts{'format'};
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
30 if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
31 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
32 die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
33 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
34
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
35 my $phred_qv=64;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
36 if (defined $opts{'phred'}) {$phred_qv=$opts{'phred'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
37
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
38 my @inputfiles=@{$opts{'i'}};
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
39 my @inputtags=@{$opts{'tag'}};
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
40
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
41 my $mypath=`pwd`;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
42 chomp $mypath;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
43
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
44 my $dir=defined $opts{'o'} ? $opts{'o'} : "$mypath/preProcess/";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
45
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
46
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
47 unless ($dir=~/\/$/) {$dir.="/";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
48 if (not -d $dir) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
49 mkdir $dir;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
50 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
51 my $config=$dir."/input_config";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
52 open CONFIG,">$config";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
53 for (my $i=0;$i<@inputfiles;$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
54 print CONFIG $inputfiles[$i],"\t",$inputtags[$i],"\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
55 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
56 close CONFIG;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
57
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
58 my $scipt_path=defined $opts{'path'} ? $opts{'path'} : "/Users/big/galaxy-dist/tools/myTools/";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
59
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
60 my $a="ATCTCGTATG"; #adapter
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
61 if (defined $opts{'a'}) {$a=$opts{'a'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
62
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
63 my $m=6; #adapter minimum mapped nt
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
64 if (defined $opts{'M'}) {$m=$opts{'M'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
65
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
66 my $t=1; #threads number
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
67 if (defined $opts{'t'}) {$t=$opts{'t'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
68
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
69 my $min_nt=19; # minimum reads length
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
70 if (defined $opts{'min'}) {$min_nt=$opts{'min'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
71
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
72 my $max_nt=28; #maximum reads length
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
73 if (defined $opts{'max'}) {$max_nt=$opts{'max'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
74
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
75 my $mis=0; #mismatch number for microRNA
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
76 if (defined $opts{'mis'}) {$mis=$opts{'mis'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
77
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
78 my $mis_rfam=0;# mismatch number for rfam
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
79 if (defined $opts{'v'}) {$mis_rfam=$opts{'v'};}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
80
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
81 my (@filein,@mark,@clean);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
82 #&read_config();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
83 @filein=@inputfiles;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
84 @mark=@inputtags;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
85
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
86 &checkfa($opts{gfa});
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
87
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
88
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
89 ##### clip adpter --> clean data start
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
90 my $preprocess=$dir."preProcess_clean/";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
91 mkdir $preprocess;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
92 my $can_use_threads = eval 'use threads; 1';
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
93 if ($can_use_threads) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
94 # Do processing using threads
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
95 print "Do processing using threads\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
96 my @filein1=@filein; my @mark1=@mark;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
97 while (@filein1>0) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
98 my @thrs; my @res;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
99 for (my $i=0;$i<$t ;$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
100 last if(@filein1==0);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
101 my $in=shift @filein1;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
102 my $out=shift @mark1;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
103 push @clean,$preprocess.$out."_clips_adapter.fq";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
104 $thrs[$i]=threads->create(\&clips,$in,$out);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
105 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
106 for (my $i=0;$i<@thrs;$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
107 $res[$i]=$thrs[$i]->join();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
108 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
109 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
110 } else {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
111 # Do not processing using threads
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
112 print "Do not processing using threads\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
113 for (my $i=0;$i<@filein ;$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
114 my $in=$filein[$i];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
115 my $out=$mark[$i];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
116 push @clean,$preprocess.$out."_clips_adapter.fq";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
117 &clips($in,$out);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
118 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
119 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
120
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
121 ##### clip adpter --> clean data end
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
122
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
123 my $collapsed=$preprocess."collapse_reads.fa";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
124 my $data=$preprocess."collapse_reads_${min_nt}_${max_nt}.fa"; ## raw clean data
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
125 &collapse(\@clean,$collapsed); #collapse reads to tags
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
126
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
127 &filterbylength(); # filter <$min_nt && >$max_nt
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
128
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
129 print "The final clean data file is $data, only contains reads which length is among $min_nt\~$max_nt\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
130
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
131
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
132 $time=Time();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
133 print "$time: Genome alignment!\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
134 my $genome_map=$dir."genome_match";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
135 &genome($data);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
136 #my $genome_map=&search($dir,"genome_match_");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
137 my $mapfile=$genome_map."/genome_mapped.bwt";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
138 my $mapfa=$genome_map."/genome_mapped.fa";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
139 my $unmap=$genome_map."/genome_not_mapped.fa";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
140
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
141 chdir $dir;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
142 my $pathfile="$dir/path.txt";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
143 open PA,">$pathfile";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
144 print PA "$config\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
145 print PA "$preprocess\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
146 print PA "$genome_map\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
147
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
148 if (defined $opts{'rfam'}) { #rfam mapping and analysis
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
149 $time=Time();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
150 print "$time: RNA annotate!\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
151 $time=~s/:/-/g;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
152 $time=~s/ /-/g;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
153 my $rfam_exp_dir=$dir."rfam_match";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
154 &rfam();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
155 #my $rfam_exp_dir=&search($dir,"rfam_match_");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
156 print PA "$rfam_exp_dir\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
157
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
158 my $tag=join "\\;" ,@mark;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
159 system("perl $scipt_path/count_rfam_express.pl -i $rfam_exp_dir/rfam_mapped.bwt -tag $tag -o rfam_non-miRNA_annotation.txt");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
160 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
161
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
162
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
163 close PA;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
164 system("perl $scipt_path/html_preprocess.pl -i $pathfile -format $format -min $min_nt -max $max_nt -o $dir/preprocessResult.html");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
165
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
166 $time=Time();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
167 print "$time: Program end!!\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
168
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
169 ############################## sub programs ###################################
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
170 sub genome{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
171 my ($file)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
172 if(defined $opts{'idx'}){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
173 system("perl $scipt_path/matching.pl -i $file -g $opts{gfa} -r 1000 -v $mis -p $t -o $dir -index $opts{idx}") ;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
174 # print "\nmatching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -index $opts{idx} -time $time\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
175 }else{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
176 system("perl $scipt_path/matching.pl -i $file -g $opts{gfa} -r 1000 -v $mis -p $t -o $dir") ;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
177 # print "\nmatching.pl -i $file -g $opts{gfa} -v $mis -p $t -r $hit -o $dir -time $time\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
178 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
179 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
180 sub rfam{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
181 if (defined $opts{'idx2'}) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
182 system("perl $scipt_path/rfam.pl -i $mapfa -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -index $opts{idx2} ");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
183 # print "\nrfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -index $opts{idx2} -time $time\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
184 }else{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
185 system("perl $scipt_path/rfam.pl -i $mapfa -ref $opts{rfam} -v $mis_rfam -p $t -o $dir ");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
186 # print "\nrfam.pl -i $data2 -ref $opts{rfam} -v $mis_rfam -p $t -o $dir -time $time\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
187 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
188 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
189 sub filterbylength{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
190 my $tmpmark=join ",", @mark;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
191 system("perl $scipt_path/filterReadsByLength.pl -i $collapsed -o $data -min $min_nt -max $max_nt -mark $tmpmark");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
192 system("perl $scipt_path/Length_Distibution.pl -i $preprocess/reads_length_distribution.txt -o $preprocess/length.html");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
193 # print "\nfilterReadsByLength.pl -i $collapsed -o $data -min $min_nt -max $max_nt -mark $tmpmark\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
194
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
195 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
196 sub collapse{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
197 my ($ins,$data)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
198 my $str="";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
199 for (my $i=0;$i<@{$ins};$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
200 $str .="-i $$ins[$i] ";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
201 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
202 system ("perl $scipt_path/collapseReads2Tags.pl $str -mark seq -o $data -format $format");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
203 # print "\ncollapseReads2Tags.pl $str -mark seq -o $data -format $format\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
204 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
205
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
206 sub clips{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
207 my ($in,$out)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
208 my $adapter=$preprocess.$out."_clips_adapter.fq";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
209 if($format eq "fq" || $format eq "fastq"){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
210 system("fastx_clipper -a $a -M $m -Q $phred_qv -i $in -o $adapter") ;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
211 # print "\nfastx_clipper -a $a -M $m -Q $phred_qv -i $in -o $adapter\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
212 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
213 if($format eq "fa" || $format eq "fasta"){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
214 system("fastx_clipper -a $a -M $m -i $in -o $adapter") ;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
215 # print "\nfastx_clipper -a $a -M $m -i $in -o $adapter\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
216 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
217 #my $clean=$preprocess.$out."_clean.fq";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
218 #system("filterReadsByLength.pl -i $adapter -o $clean -min $min_nt -max $max_nt ");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
219
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
220 return;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
221 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
222
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
223 sub read_config{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
224 open CON,"<$config";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
225 while (my $aline=<CON>) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
226 chomp $aline;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
227 my @tmp=split/\t/,$aline;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
228 push @filein,$tmp[0];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
229 push @mark,$tmp[1];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
230 &check_rawdata($tmp[0]);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
231 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
232 close CON;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
233 if (@filein != @mark) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
234 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
235 die "Maybe config file have some wrong!!!\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
236 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
237 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
238 sub check_rawdata{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
239 my ($fileforcheck)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
240 if (!(-s $fileforcheck)) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
241 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
242 die "Can not find $fileforcheck, or file is empty!!!\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
243 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
244 if ($format eq "fasta" || $format eq "fa") {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
245 &checkfa($fileforcheck);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
246 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
247 if ($format eq "fastq" || $format eq "fq") {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
248 &checkfq($fileforcheck);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
249 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
250 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
251 sub checkfa{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
252 my ($file_reads)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
253 open N,"<$file_reads";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
254 my $line=<N>;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
255 chomp $line;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
256 if($line !~ /^>\S+/){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
257 #printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
258 die "The first line of file $file_reads does not start with '>identifier'
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
259 Reads file $file_reads is not a valid fasta file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
260 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
261 if(<N> !~ /^[ACGTNacgtn]*$/){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
262 #printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
263 die "File $file_reads contains not allowed characters in sequences
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
264 Allowed characters are ACGTN
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
265 Reads file $file_reads is not a fasta file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
266 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
267 close N;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
268 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
269 sub checkfq{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
270 my ($file_reads)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
271
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
272 open N,"<$file_reads";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
273 for (my $i=0;$i<10;$i++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
274 my $a=<N>;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
275 my $b=<N>;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
276 my $c=<N>;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
277 my $d=<N>;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
278 chomp $a;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
279 chomp $b;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
280 chomp $c;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
281 chomp $d;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
282 if($a!~/^\@/){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
283 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
284 die "$file_reads is not a fastq file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
285 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
286 if($b!~ /^[ACGTNacgtn]*$/){
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
287 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
288 die "File $file_reads contains not allowed characters in sequences
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
289 Allowed characters are ACGTN
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
290 Reads file $file_reads is not a fasta file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
291 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
292 if ($c!~/^\@/ && $c!~/^\+/) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
293 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
294 die "$file_reads is not a fastq file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
295 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
296 if ((length $b) != (length $d)) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
297 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
298 die "$file_reads is not a fastq file\n\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
299 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
300 my @qv=split //,$d;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
301 for (my $j=0;$j<@qv ;$j++) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
302 my $q=ord($qv[$j])-64;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
303 if($q<0){$phred_qv=33;}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
304 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
305 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
306 close N;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
307 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
308
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
309 sub search{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
310 my ($dir,$str)=@_;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
311 opendir I,$dir;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
312 my @ret;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
313 while (my $file=readdir I) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
314 if ($file=~/$str/) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
315 push @ret, $file;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
316 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
317 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
318 closedir I;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
319 if (@ret != 1) {
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
320 #&printErr();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
321
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
322 die "Can not find directory or file which name has string: $str !!!\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
323 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
324 return $ret[0];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
325 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
326
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
327 sub Time{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
328 my $time=time();
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
329 my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
330 $month++;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
331 $year+=1900;
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
332 if (length($sec) == 1) {$sec = "0"."$sec";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
333 if (length($min) == 1) {$min = "0"."$min";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
334 if (length($hour) == 1) {$hour = "0"."$hour";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
335 if (length($day) == 1) {$day = "0"."$day";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
336 if (length($month) == 1) {$month = "0"."$month";}
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
337 #print "$year-$month-$day $hour:$min:$sec\n";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
338 return("$year-$month-$day $hour:$min:$sec");
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
339 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
340
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
341
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
342 sub usage{
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
343 print <<"USAGE";
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
344 Version $version
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
345 Usage:
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
346 $0 -i -format -gfa -index -rfam -a -M -min -max -mis -v -t -o -path
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
347 options:
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
348 -i input files, # raw data file, can be multipe eg. -i xxx.fq -i xxx .fq ...
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
349 -tag string # raw data file names, -tag xxx -tag xxx
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
350
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
351 -format string,#specific input rawdata file format : fastq|fq|fasta|fa
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
352 -phred int # phred quality number, default is 64
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
353
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
354 -path scirpt path
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
355
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
356 -gfa string, input file # genome fasta. sequence file
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
357 -idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
358 string must be the prefix of the bowtie index. For instance, if
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
359 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
360 the prefix is 'h_sapiens_37_asm'.##can be null
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
361
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
362 -rfam string, input file# rfam database file, microRNAs must not be contained in this file## if not define, rfam small RNA will not be count.
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
363 -idx2 string, rfam file index, file-prefix #(must be indexed by bowtie-build) The parameter
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
364 string must be the prefix of the bowtie index. For instance, if
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
365 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
366 the prefix is 'h_sapiens_37_asm'.##can be null
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
367
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
368 -a string, ADAPTER string. default is ATCTCGTATG.
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
369 -M int, require minimum adapter alignment length of N. If less than N nucleotides aligned with the adapter - don't clip it.
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
370 -min int, reads min length,default is 19.
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
371 -max int, reads max length,default is 28.
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
372
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
373 -mis [int] number of allowed mismatches when mapping reads to genome, default 0
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
374 -v <int> report end-to-end hits w/ <=v mismatches; ignore qualities,default 0; used in rfam alignment
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
375
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
376 -t int, number of threads [1]
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
377
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
378 -o output directory# absolute path
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
379 -h help
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
380 USAGE
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
381 exit(1);
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
382 }
c75593f79aa9 Uploaded
big-tiandm
parents:
diff changeset
383