Mercurial > repos > bigrna > gpsrna

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