Mercurial > repos > big-tiandm > sirna_plant

comparison siRNA.pl @ 0:07745c0958dd draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Thu, 18 Sep 2014 21:40:25 -0400
parents
children f466394ee1fd
comparison
equal deleted inserted replaced
-1:000000000000 0:07745c0958dd
1 #!/usr/bin/perl -w
2 my $version=1.00;
3 use strict;
4 use warnings;
5 use Getopt::Long;
6 use Getopt::Std;
7 use threads;
8 use threads::shared;
9 use Parallel::ForkManager;
10 use lib '/leofs/biotrans/chentt/perl_module/';
11 #perl ../siRNA.pl -i config -g /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/genome.fa -f /share_bio/hs4/disk3-4/Reference/Plants/Rice_TIGR/Reference/TIGR/version_6.1/all.dir/all.gff3 -path /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/ -o /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test -t 3 -rfam /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/Rfam.fasta -idx /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/genome -idx2 /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/rfam -deg deg -n 25 -nat class/nat_1 -repeat class/repeat_1 -cen centromere_TIGR.txt -format fastq
12 print "
13 #####################################
14 # #
15 # sRNA cluster #
16 # #
17 #####################################
18 ";
19 ###########################################################################################
20 my $usage="$0
21 Options:
22 -i input file# raw data file
23 -tag string #raw data sample name
24 -g genome file
25 -f gff file
26
27 -o workdir file
28 -path script path
29 -t int, number of threads [1]
30 -format fastq, fq, fasta or fa
31 -idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter
32 string must be the prefix of the bowtie index. For instance, if
33 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
34 the prefix is 'h_sapiens_37_asm'.##can be null
35 -mis int number of allowed mismatches when mapping reads to genome, default 0
36 -rfam string, input file# rfam database file.
37 -idx2 string, rfam file index, file-prefix #(must be indexed by bowtie-build) The parameter
38 string must be the prefix of the bowtie index. For instance, if
39 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
40 the prefix is 'h_sapiens_37_asm'.##can be null
41
42 -v int report end-to-end hits w/ <=v mismatches; ignore qualities,default 0; used in rfam alignment
43
44 -a string, ADAPTER string. default is ATCTCGTATG.
45 -n int max hits number,default 25; used in genome alignment
46 -d int distance of tag to merged a cluster; default 100
47 -p cluster method F :conventional default is F
48 T :NIBLES
49 -l int the length of the upstream and downstream,default 1000;used in position annotate
50
51 -nat natural antisense transcripts file
52 -repeat repeat information file out of Repeatmasker
53 -deg file config of de sample
54 -cen centromere file input
55 -span plot span, default 50000
56 ";
57
58 my %options;
59 GetOptions(\%options,"i:s@","tag:s@","g=s","f=s","o=s","a:s","path:s","p=s","format=s","nat:s","repeat:s","deg:s","n:i","mis:i","rfam:s","t:i","v:i","d:i","l:i","idx:s","idx2:s","cen:s","span:s","h");
60
61 my @inputfiles=@{$opts{'i'}};
62 my @inputtags=@{$opts{'tag'}};
63
64 #my $config=$options{'i'};
65 my $genome_fa=$options{'g'};
66 my $gff=$options{'f'};
67 ##########################################################################################
68 my $predir=`pwd`;
69 chomp $predir;
70 my $workdir=defined($options{'o'}) ? $options{'o'}:$predir;
71
72 my $path=$options{'path'};
73
74 my $t=defined($options{'t'})? $options{'t'}:1; #threads number
75
76 my $mis=defined $options{'mis'} ? $options{'mis'}:0;
77
78 my $mis_rfam=defined $options{'v'} ? $options{'v'}:0;
79
80 my $hit=defined $options{'n'}?$options{'n'}:25;
81
82 my $distance_of_merged_tag=defined $options{'d'} ? $options{'d'}:100;
83
84 my $up_down_dis=defined $options{'l'} ?$options{'l'}:1000;
85
86 my $cluster_mothod=defined $options{'p'}?$options{'p'}:"F";
87
88 my $format=$options{'format'};
89 #if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") {
90 # die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
91 #}
92
93 my $adpter="ATCTCGTATG"; #adapter
94 if (defined $opts{'a'}) {$a=$opts{'a'};}
95
96 #print help if that option is used
97 if($options{h}){die $usage;}
98
99 my $phred_qv=64;
100 my $sample_number;
101 my ($dir,$dir_tmp);
102 ################################ MAIN ##################################################
103 print "\ncluster program start:";
104 my $time=Time();
105 make_dir_tmp();
106
107 my (@filein,@mark,@clip);
108 my $mark;
109 my $sample_mark;
110
111 my $config=$workdir."/input_config";
112 open CONFIG,">$config";
113 for (my $i=0;$i<@inputfiles;$i++) {
114 print CONFIG $inputfiles[$i],"\t",$inputtags[$i],"\n";
115 }
116 close CONFIG;
117
118 read_config();
119
120 trim_adapter_and_filter();
121
122 my $filter_out=$dir."preProcess\/"."collapse_reads_out.fa";## raw clean data
123 my $data2=$filter_out; ### mirbase not mapped reads
124 my $data3=$dir."\/rfam_match\/rfam_not_mapped\.fa"; ### rfam not mapped reads
125 my $bed=$dir."cluster\/"."sample\.bed";
126 my $read=$dir."cluster\/"."sample_reads\.cluster";
127 my $read_txt=$dir."cluster\/"."cluster\.txt";
128 my $rpkm=$dir."cluster\/"."sample_rpkm\.cluster";
129 my $preprocess;
130 my $cluster_file;
131 my $annotate_dir;
132 my $deg_dir;
133 my %id;
134 for (my $i=0;$i<@mark ;$i++) {
135 $id{$mark[$i]}=$i+4;
136 }
137 group_and_filter(); #collapse reads to tags
138
139 rfam();
140
141 my @map_read;
142 my $map_tag=0;
143 genome();
144
145 bwt2bed();
146
147 cluster();
148
149 quantify();
150
151 phase();
152
153 class();
154
155 annotate();
156
157 genome_length();
158
159 plot();
160
161 my @pairdir;
162 if (defined($options{'deg'})) {
163 dec();
164 infor_merge();
165 }
166 html();
167 print "\ncluster program end:";
168 Time();
169 ############################sub program###################################################
170 sub make_dir_tmp{
171
172 #make temporary directory
173 if(not -d "$workdir\/cluster_runs_$time"){
174 mkdir("$workdir\/cluster_runs_$time");
175 mkdir("$workdir\/cluster_runs_$time\/ref\/");
176 }
177
178 $dir="$workdir\/cluster_runs_$time\/";
179 print STDERR "mkdir $dir\n\n";
180 return;
181 }
182
183 sub read_config{
184 open IN,"<$config";
185 while (my $aline=<IN>) {
186 chomp $aline;
187 my @tmp=split/\t/,$aline;
188 push @filein,$tmp[0];
189 push @mark,$tmp[1];
190 }
191 close IN;
192 if (@filein != @mark) {
193 die "Maybe config file have some wrong!!!\n";
194 }
195 $sample_number=@mark;
196 $mark=join "\t",@mark;
197 $sample_mark=join "\#",@mark;
198 }
199
200
201 sub trim_adapter_and_filter{
202 my $time=time();
203 $preprocess=$dir."preProcess/";
204 mkdir $preprocess;
205 my $can_use_threads = eval 'use threads; 1';
206 if ($can_use_threads) {
207 # Do processing using threads
208 my @filein1=@filein; my @mark1=@mark;
209 while (@filein1>0) {
210 my @thrs; my @res;
211 for (my $i=0;$i<$t ;$i++) {
212 last if(@filein1==0);
213 my $in=shift @filein1;
214 my $out=shift @mark1;
215 push @clip,$dir."preProcess\/$out\_clip\.fq";
216 $thrs[$i]=threads->create(\&clips,$in,$out);
217 }
218 for (my $i=0;$i<@thrs;$i++) {
219 $res[$i]=$thrs[$i]->join();
220 }
221 }
222 }
223 else {
224 # Do not processing using threads
225 for (my $i=0;$i<@filein ;$i++) {
226 my $in=$filein[$i];
227 my $out=$mark[$i];
228 push @clip,$dir."preProcess\/$out\_clip\.fq";
229 &clips($in,$out);
230 }
231 }
232 }
233
234 sub clips{
235 my ($filein,$fileout)=@_;
236 my $adapter=$dir."preProcess\/$fileout\_clip\.fq";
237 if($format eq "fq" || $format eq "fastq"){
238 my $clip=`$path\/fastx_clipper -a $adpter -M 6 -Q $phred_qv -i $filein -o $adapter`;
239 }
240 if($format eq "fa" || $format eq "fasta"){
241 my $clip=`$path\/fastx_clipper -a $adpter -M 6 -i $filein -o $adapter`;
242 }
243 #my $clean=$dir."preProcess\/$fileout\_clean.fq";
244 #my $filter=`filterReadsByLength.pl -i $adapter -o $clean -min 18 -max 40 `;
245 return $fileout;
246 }
247
248 sub group_and_filter{
249 #my ($ins,$data)=@_;
250 my @ins=@clip;
251 my $str="";
252 my $group_out_file=$dir."preProcess\/"."collapse_reads.fa";
253 #print "$$ins[0]\t$$ins[0]\n";
254 for (my $i=0;$i<@clip;$i++) {
255 $str .="-i $clip[$i] ";
256 #print "$$ins[$i]\n";
257 }
258 my $group=`perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format`;
259 print "perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format\n\n";
260
261 my $l_out=$dir."preProcess\/"."collapse_reads_18-40.fa";
262 my $length_f=`perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark`;
263 print "perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark\n\n";
264 my $cout_f=`perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark`;
265 print "perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark\n\n";
266 return 0;
267 }
268
269 sub rfam{
270 if (defined $options{'idx2'}) {
271 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir -index $options{idx2}");
272 }else{
273 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir");
274 }
275 my $tag=join "\\;" ,@mark;
276 my $rfam_count=`perl $path\/count_rfam_express.pl -i $dir\/rfam_match\/rfam_mapped.bwt -tag $tag -o $dir\/rfam_match\/rfam_non-miRNA_annotation.txt`;
277 return 0;
278 }
279 sub genome{
280 if(defined $options{'idx'}){
281 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir -index $options{idx}") ;
282 }else{
283 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir ") ;
284 }
285 #=================== mapping sta ===================================================
286 my $map_file=$dir."genome_match\/genome_mapped\.fa";
287 open (MAP,"<$map_file")||die"$!";
288 print "\n#each sample mapping reads sta:\n\n";
289 print "#$mark\ttotal\n";
290 while (my $ID=<MAP>) {
291 chomp $ID;
292 my @tmp=split/\:/,$ID;
293 my @exp=split/\_/,$tmp[1];
294 $exp[-1] =~ s/^x//;
295 for (my $i=0;$i<@exp ;$i++) {
296 $map_read[$i]+=$exp[$i];
297 }
298 $map_tag++;
299 my $seq=<MAP>;
300 }
301 my $map_read=join"\t",@map_read;
302 print "$map_read\n\n";
303 print "#total mapped tags:$map_read\n\n";
304 close MAP;
305 return 0;
306 }
307
308 sub bwt2bed{
309 $cluster_file=$dir."cluster\/";
310 mkdir ("$cluster_file");
311 print "sam file changed to bed file\n";
312 my ($file) = $dir."genome_match\/genome_mapped\.bwt";
313
314 my $sam2bed=`perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed `;
315 print "perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed\n\n";
316 return 0;
317 }
318
319 sub cluster{
320 print "tags is ready to merged clusters\n\n";
321 my ($file) =$bed;
322 if ($cluster_mothod eq "F") {
323 my $cluster=`perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt`;
324 print "Using converntional method\n perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt\n\n";
325 }
326 elsif($cluster_mothod eq "T"){
327 my $cluster=`perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $read_txt -mark $sample_mark`;
328 print "Using nibls method\n perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $dir\/cluster.txt -mark $sample_mark\n\n";
329 }
330 else{print "\-p is wrong!\n\n";}
331 return 0;
332 }
333
334
335 sub quantify{
336 print "clusters is ready to quantified\n\n";
337 my @depth=@map_read;
338 pop @depth;
339 my $depth=join ",",@depth;
340 my $quantify=`perl $path\/quantify.pl -i $read -d $depth -o $rpkm`;
341 print "perl $path\/quantify.pl -i $read -d $depth -o $rpkm\n\n\n";
342 return 0;
343 }
344
345 sub phase{
346 $annotate_dir=$dir."annotate\/";
347 mkdir ("$annotate_dir");
348 print "clusters is to predict phase siRNA\n";
349 my $phase=`perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out`;
350 print "perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out\n\n\n";
351 return 0;
352 }
353
354 sub class{
355 print "clusters is ready to annotate by source\n\n";
356 my $nat=$options{'nat'};
357 my $repeat=$options{'repeat'};
358 my $class=`perl $path\/ClassAnnotate.pl -i $rpkm -g $gff -n $nat -r $repeat -p $annotate_dir\/phase.out -o $annotate_dir\/sample_class.anno -t $annotate_dir\/nat.out -l $dir\/ref\/genelist.txt`;
359 print "perl $path\/ClassAnnotate.pl -i $rpkm -g $gff -n $nat -r $repeat -p $annotate_dir\/phase.out -o $annotate_dir\/sample_class.anno -t $annotate_dir\/nat.out -l $dir\/ref\/genelist.txt\n\n";
360 }
361
362 sub annotate{
363 print "clusters is ready to annotate by gff file\n\n";
364 my $file="$annotate_dir\/sample_class.anno";
365 my $annotate=`perl $path\/Annotate.pl -i $file -g $dir\/ref\/genelist.txt -d $up_down_dis -o $annotate_dir\/sample_c_p.anno`;
366 print "perl $path\/Annotate.pl -i $file -g $dir\/ref\/genelist.txt -d $up_down_dis -o $annotate_dir\/sample_c_p.anno\n\n";
367 return 0;
368 }
369
370 sub dec{
371 print "deg reading\n\n";
372 my $deg_file=$options{'deg'};
373 open IN,"<$deg_file";
374 my @deg;
375 my $s=0;
376 while (my $aline=<IN>) {
377 chomp $aline;
378 next if($aline=~/^\#/);
379 $deg[$s]=$aline;
380 my @ea=split/\s+/,$aline;
381 push @pairdir,"$ea[0]_VS_$ea[1]\/";
382 #print "$deg[$s]\n";
383 $s++;
384 }
385 close IN;
386 $deg_dir=$dir."deg\/";
387 mkdir ("$deg_dir");
388 my $max_process = 10;
389 my $pm = new Parallel::ForkManager( $max_process );
390 my $number=@deg-1;
391 foreach(0..$number){
392 $pm->start and next;
393 &dec_pel($deg[$_]);
394 $pm->finish;
395 }
396 $pm->wait_all_children;
397 }
398
399 sub dec_pel{
400 print "start:\n";
401 Time();
402 my $sample=shift(@_);
403 my @each=split/\s+/,$sample;
404 print "$each[0]\t$each[1]\n";
405 my $deg_sample_dir=$deg_dir."$each[0]_VS_$each[1]\/";
406 mkdir ("$deg_sample_dir");
407 my $deg=`perl $path\/DEGseq_2.pl -i $read -outdir $deg_sample_dir -column1 $id{$each[0]} -mark1 $each[0] -column2 $id{$each[1]} -mark2 $each[1]`; #-depth1 -depth2
408 my $time2=time();
409 print "end:\n";
410 Time();
411 sleep 1;
412 }
413
414 sub infor_merge{
415 my ($input,$mark);
416 foreach (@pairdir) {
417 print "@pairdir\n";
418 $mark.=" -mark $_ ";
419 $input.=" -i $dir/deg\/$_\/output_score\.txt ";
420 print "$input\n$mark\n";
421 }
422 my $infor_merge=`perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result `;
423 print "perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result\n\n";
424 }
425
426 sub genome_length{
427 my $length=`perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length`;
428 print "perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length\n\n"
429
430 }
431
432 sub plot{
433 my $plot_file="$dir\/plot\/";
434 mkdir ("$plot_file");
435 my $genome_plot="$dir\/plot\/genome\/";
436 mkdir ("$genome_plot");
437 #genome cluster
438 my $span=defined($options{span})?$options{span}:50000;
439 foreach (1..$sample_number) {
440 my $mark=$mark[$_-1];
441 my $cen="";
442 if (defined $options{cen}) {
443 $cen="-cen $options{cen}";
444 }
445 my $plot=`perl $path\/sRNA_rpkm_distribution_along_genome.pl -c $rpkm -n $_ -mark $mark -span $span -l $dir\/ref\/genome\.length $cen -o $genome_plot\/$mark\.html -out $genome_plot\/$mark\.txt`;
446 print "perl $path\/sRNA_rpkm_distribution_along_genome.pl -c $rpkm -n $_ -mark $mark -span $span -l $dir\/ref\/genome\.length $cen -o $genome_plot\/$mark\.html -out $genome_plot\/$mark\.txt\n\n";
447 }
448
449 my $chr_plot_dir="$dir\/plot\/chr\/";
450 mkdir("$chr_plot_dir");
451 my %chr;
452 open LEN,"<$dir\/ref\/genome\.length";
453 while (my $aline=<LEN>) {
454 next if($aline=~/^\#/);
455 chomp $aline;
456 my @temp=split/\t/,$aline;
457 $chr{$temp[0]}=$temp[1];
458 }
459 close LEN;
460 foreach my $chr (sort keys %chr) {
461 my $cen="";
462 if (defined $options{cen}) {
463 $cen="-cen $options{cen}";
464 }
465 my $chr_plot=`perl $path\/chr_plot.pl -l $chr{$chr} -chro $chr -g $dir\/ref\/genelist.txt -span $span -c $rpkm -mark $sample_mark -o $chr_plot_dir\/$chr\.html`;
466 print "perl $path\/chr_plot.pl -l $chr{$chr} -chro $chr -g $dir\/ref\/genelist.txt -span $span -c $rpkm -mark $sample_mark -o $chr_plot_dir\/$chr\.html\n";
467 }
468 }
469
470 sub html{
471 my $pathfile="$dir/path.txt";
472 open PA,">$pathfile";
473 print PA "$config\n";
474 print PA "$preprocess\n";
475 print PA "$dir"."rfam_match\n";
476 print PA "$dir"."genome_match\n";
477 print PA "$cluster_file\n";
478 print PA "$annotate_dir\n";
479 print PA "$deg_dir\n";
480 close PA;
481 my $html=`perl $path\/html.pl -i $pathfile -format $format -o $dir/result.html`;
482 }
483
484 sub Time{
485 my $time=time();
486 my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
487 $month++;
488 $year+=1900;
489 if (length($sec) == 1) {$sec = "0"."$sec";}
490 if (length($min) == 1) {$min = "0"."$min";}
491 if (length($hour) == 1) {$hour = "0"."$hour";}
492 if (length($day) == 1) {$day = "0"."$day";}
493 if (length($month) == 1) {$month = "0"."$month";}
494 print "$year-$month-$day $hour:$min:$sec\n";
495 return("$year-$month-$day-$hour-$min-$sec");
496 }
497 #################################################################################