Mercurial > repos > big-tiandm > mirplant2

comparison siRNA_pipeline.pl @ 50:7b5a48b972e9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Fri, 05 Dec 2014 00:11:02 -0500
parents
children
comparison
equal deleted inserted replaced
49:f008ab2cadc6 50:7b5a48b972e9
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","phed:i","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 #print help if that option is used
61 if($options{h}){die $usage;}
62
63 my @filein=@{$options{'i'}};
64 my @mark=@{$options{'tag'}};
65
66 #my $config=$options{'i'};
67 my $genome_fa=$options{'g'};
68 my $gff=$options{'f'};
69
70
71 ##########################################################################################
72 my $predir=`pwd`;
73 chomp $predir;
74 my $workdir=defined($options{'o'}) ? $options{'o'}:$predir;
75
76 my $path=$options{'path'};
77
78 my $t=defined($options{'t'})? $options{'t'}:1; #threads number
79
80 my $mis=defined $options{'mis'} ? $options{'mis'}:0;
81
82 my $mis_rfam=defined $options{'v'} ? $options{'v'}:0;
83
84 my $hit=defined $options{'n'}?$options{'n'}:25;
85
86 my $distance_of_merged_tag=defined $options{'d'} ? $options{'d'}:100;
87
88 my $up_down_dis=defined $options{'l'} ?$options{'l'}:1000;
89
90 my $cluster_mothod=defined $options{'p'}?$options{'p'}:"F";
91
92 my $format=$options{'format'};
93 #if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") {
94 # die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
95 #}
96
97 my $adpter="ATCTCGTATG"; #adapter
98 if (defined $options{'a'}) {$adpter=$options{'a'};}
99
100
101 my $phred_qv=64;
102 if(defined $opts{'phred'}){$phred_qv=$opts{'phred'};}
103 my $sample_number;
104 my ($dir,$dir_tmp);
105 ################################ MAIN ##################################################
106 print "\ncluster program start:";
107 my $time=Time();
108 make_dir_tmp();
109
110 my @clip;
111 my $mark;
112 my $sample_mark;
113
114 my $config=$dir."/input_config";
115 open CONFIG,">$config";
116 for (my $i=0;$i<@filein;$i++) {
117 print CONFIG $filein[$i],"\t",$mark[$i],"\n";
118 }
119 close CONFIG;
120 if (@filein != @mark) {
121 die "Maybe config file have some wrong!!!\n";
122 }
123 $sample_number=@mark;
124 $mark=join "\t",@mark;
125 $sample_mark=join "\#",@mark;
126
127
128 #read_config();
129
130 trim_adapter_and_filter();
131
132 my $filter_out=$dir."preProcess\/"."collapse_reads_out.fa";## raw clean data
133 my $data2=$filter_out; ### mirbase not mapped reads
134 my $data3=$dir."\/rfam_match\/rfam_not_mapped\.fa"; ### rfam not mapped reads
135 my $bed=$dir."cluster\/"."sample\.bed";
136 my $read=$dir."cluster\/"."sample_reads\.cluster";
137 my $read_txt=$dir."cluster\/"."cluster\.txt";
138 my $rpkm=$dir."cluster\/"."sample_rpkm\.cluster";
139 my $preprocess;
140 my $cluster_file;
141 my $annotate_dir;
142 my $deg_dir;
143 my $plot_dir;
144 my %id;
145 for (my $i=0;$i<@mark ;$i++) {
146 $id{$mark[$i]}=$i+4;
147 }
148
149 print "\n######## tiandm test start ###########\n";
150 print "\@mark: @mark\n\%id keys number:";
151 print scalar keys %id;
152 print "\n";
153 foreach my $kyess (keys %id){
154 print $kyess," --> $id{$kyess}\n";
155 }
156 print "\n######## tiandm test end ############\n";
157 group_and_filter(); #collapse reads to tags
158
159 rfam();
160
161 my @map_read;
162 my $map_tag=0;
163 genome();
164
165 bwt2bed();
166
167 cluster();
168
169 quantify();
170
171 phase();
172
173 if (defined($options{'nat'})&&defined($options{'repeat'})) {
174 class();
175 }
176 else{
177 get_genelist();
178 }
179
180 annotate();
181
182 genome_length();
183
184 plot();
185
186 my @pairdir;
187 if (defined($options{'deg'})) {
188 dec();
189 infor_merge();
190 }
191 else{infor_merge_no_dec()}
192 html();
193 print "\ncluster program end:";
194 Time();
195 ############################sub program###################################################
196 sub make_dir_tmp{
197
198 #make temporary directory
199 if(not -d "$workdir\/cluster_runs"){
200 mkdir("$workdir\/cluster_runs");
201 mkdir("$workdir\/cluster_runs\/ref\/");
202 }
203
204 $dir="$workdir\/cluster_runs\/";
205 #print STDERR "mkdir $dir\n\n";
206 return;
207 }
208
209 #sub read_config{
210 # open IN,"<$config";
211 # while (my $aline=<IN>) {
212 # chomp $aline;
213 # my @tmp=split/\t/,$aline;
214 # push @filein,$tmp[0];
215 # push @mark,$tmp[1];
216 # }
217 # close IN;
218 # if (@filein != @mark) {
219 # die "Maybe config file have some wrong!!!\n";
220 # }
221 # $sample_number=@mark;
222 # $mark=join "\t",@mark;
223 # $sample_mark=join "\#",@mark;
224 #}
225
226
227 sub trim_adapter_and_filter{
228 my $time=time();
229 $preprocess=$dir."preProcess/";
230 mkdir $preprocess;
231 my $can_use_threads = eval 'use threads; 1';
232 if ($can_use_threads) {
233 # Do processing using threads
234 my @filein1=@filein; my @mark1=@mark;
235 while (@filein1>0) {
236 my @thrs; my @res;
237 for (my $i=0;$i<$t ;$i++) {
238 last if(@filein1==0);
239 my $in=shift @filein1;
240 my $out=shift @mark1;
241 push @clip,$dir."preProcess\/$out\_clip\.fq";
242 $thrs[$i]=threads->create(\&clips,$in,$out);
243 }
244 for (my $i=0;$i<@thrs;$i++) {
245 $res[$i]=$thrs[$i]->join();
246 }
247 }
248 }
249 else {
250 # Do not processing using threads
251 for (my $i=0;$i<@filein ;$i++) {
252 my $in=$filein[$i];
253 my $out=$mark[$i];
254 push @clip,$dir."preProcess\/$out\_clip\.fq";
255 &clips($in,$out);
256 }
257 }
258 }
259
260 sub clips{
261 my ($filein,$fileout)=@_;
262 my $adapter=$dir."preProcess\/$fileout\_clip\.fq";
263 if($format eq "fq" || $format eq "fastq"){
264 my $clip=`fastx_clipper -a $adpter -M 6 -Q $phred_qv -i $filein -o $adapter`;
265 }
266 if($format eq "fa" || $format eq "fasta"){
267 my $clip=`fastx_clipper -a $adpter -M 6 -i $filein -o $adapter`;
268 }
269 #my $clean=$dir."preProcess\/$fileout\_clean.fq";
270 #my $filter=`filterReadsByLength.pl -i $adapter -o $clean -min 18 -max 40 `;
271 return $fileout;
272 }
273
274 sub group_and_filter{
275 #my ($ins,$data)=@_;
276 my @ins=@clip;
277 my $str="";
278 my $group_out_file=$dir."preProcess\/"."collapse_reads.fa";
279 #print "$$ins[0]\t$$ins[0]\n";
280 for (my $i=0;$i<@clip;$i++) {
281 $str .="-i $clip[$i] ";
282 #print "$$ins[$i]\n";
283 }
284 my $group=`perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format`;
285 print "perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format\n\n";
286
287 my $l_out=$dir."preProcess\/"."collapse_reads_18-40.fa";
288 my $tmpmark=join ",", @mark;
289
290 my $length_f=`perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $tmpmark`;
291 print "perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $tmpmark\n\n";
292 my $cout_f=`perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark`;
293 print "perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark\n\n";
294 my $plot_l_D=`perl $path/Length_Distibution.pl -i $dir/preProcess/reads_length_distribution_after_count_filter.txt -o $dir/preProcess/length.html `;
295 print "perl $path\/Length_Distibution.pl -i $dir\/preProcess\/reads_length_distribution_after_count_filter.txt -o $dir\/preProcess\/length\.html\n\n";
296 return 0;
297 }
298
299 sub rfam{
300 if (defined $options{'idx2'}) {
301 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir -index $options{idx2}");
302 }else{
303 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir");
304 }
305 my $tag=join "\\;" ,@mark;
306 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`;
307 return 0;
308 }
309 sub genome{
310 if(defined $options{'idx'}){
311 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir -index $options{idx}") ;
312 }else{
313 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir ") ;
314 }
315 #=================== mapping sta ===================================================
316 my $map_file=$dir."genome_match\/genome_mapped\.fa";
317 open (MAP,"<$map_file")||die"$!";
318 print "\n#each sample mapping reads sta:\n\n";
319 print "#$mark\ttotal\n";
320 while (my $ID=<MAP>) {
321 chomp $ID;
322 my @tmp=split/\:/,$ID;
323 my @exp=split/\_/,$tmp[1];
324 $exp[-1] =~ s/^x//;
325 for (my $i=0;$i<@exp ;$i++) {
326 $map_read[$i]+=$exp[$i];
327 }
328 $map_tag++;
329 my $seq=<MAP>;
330 }
331 my $map_read=join"\t",@map_read;
332 print "$map_read\n\n";
333 print "#total mapped tags:$map_read\n\n";
334 close MAP;
335 return 0;
336 }
337
338 sub bwt2bed{
339 $cluster_file=$dir."cluster\/";
340 mkdir ("$cluster_file");
341 print "sam file changed to bed file\n";
342 my ($file) = $dir."genome_match\/genome_mapped\.bwt";
343
344 my $sam2bed=`perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed `;
345 print "perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed\n\n";
346 return 0;
347 }
348
349 sub cluster{
350 print "tags is ready to merged clusters\n\n";
351 my ($file) =$bed;
352 if ($cluster_mothod eq "F") {
353 my $cluster=`perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt`;
354 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";
355 }
356 elsif($cluster_mothod eq "T"){
357 my $cluster=`perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $read_txt -k $sample_mark`;
358 print "Using nibls method\n perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $dir\/cluster.txt -k $sample_mark\n\n";
359 }
360 else{print "\-p is wrong!\n\n";}
361 return 0;
362 }
363
364
365 sub quantify{
366 print "clusters is ready to quantified\n\n";
367 my @depth=@map_read;
368 pop @depth;
369 my $depth=join ",",@depth;
370 my $quantify=`perl $path\/quantify.pl -i $read -d $depth -o $rpkm`;
371 print "perl $path\/quantify_siRNA.pl -i $read -d $depth -o $rpkm\n\n\n";
372 return 0;
373 }
374
375 sub phase{
376 $annotate_dir=$dir."annotate\/";
377 mkdir ("$annotate_dir");
378 print "clusters is to predict phase siRNA\n";
379 my $phase=`perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out`;
380 print "perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out\n\n\n";
381 return 0;
382 }
383
384 sub class{
385 print "clusters is ready to annotate by sources\n\n";
386 my $nat=$options{'nat'};
387 my $repeat=$options{'repeat'};
388 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`;
389 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";
390 }
391
392 sub annotate{
393 print "clusters is ready to annotate by gff file\n\n";
394 my $file;
395 if (defined($options{'nat'})&&defined($options{'repeat'})) {
396 $file="$annotate_dir\/sample_class.anno";
397 }
398 else{
399 $file=$rpkm;
400 }
401 my $annotate=`perl $path\/Annotate.pl -i $file -g $dir\/ref\/genelist.txt -d $up_down_dis -o $annotate_dir\/sample_c_p.anno`;
402 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";
403 return 0;
404 }
405 sub get_genelist{
406
407 my $get_genelist=`perl $path\/get_genelist.pl -i $gff -o $dir\/ref\/genelist.txt`;
408 print "perl $path\/get_genelist.pl -i $gff -o $dir\/ref\/genelist.txt";
409 }
410
411 sub dec{
412 print "deg reading\n\n";
413 my $deg_file=$options{'deg'};
414 open IN,"<$deg_file";
415 my @deg;
416 my $s=0;
417 while (my $aline=<IN>) {
418 chomp $aline;
419 next if($aline=~/^\#/);
420 $deg[$s]=$aline;
421 my @ea=split/\s+/,$aline;
422 push @pairdir,"$ea[0]_VS_$ea[1]\/";
423 #print "$deg[$s]\n";
424 $s++;
425 }
426 close IN;
427 $deg_dir=$dir."deg\/";
428 mkdir ("$deg_dir");
429 my $max_process = 10;
430 my $pm = new Parallel::ForkManager( $max_process );
431 my $number=@deg-1;
432 foreach(0..$number){
433 $pm->start and next;
434 &dec_pel($deg[$_]);
435 $pm->finish;
436 }
437 $pm->wait_all_children;
438 }
439
440 sub dec_pel{
441 print "\n******************\nstart:\n";
442 Time();
443 my $sample=shift(@_);
444 my @each=split/\s+/,$sample;
445 print "$each[0]\t$each[1]\n";
446 my $deg_sample_dir=$deg_dir."$each[0]_VS_$each[1]\/";
447 mkdir ("$deg_sample_dir");
448 print "read: $read\n";
449 print "deg_sample_dir: $deg_sample_dir\n";
450 print "$id{$each[0]}\t$each[0]\n";
451 print "$id{$each[1]}\t$each[1]\n";
452 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
453 my $time2=time();
454 print "end:\n*************************\n";
455 Time();
456 sleep 1;
457 }
458
459 sub infor_merge{
460 my ($input,$mark);
461 foreach (@pairdir) {
462 print "@pairdir\n";
463 $mark.=" -mark $_ ";
464 $input.=" -i $dir/deg\/$_\/output_score\.txt ";
465 print "$input\n$mark\n";
466 }
467 my $infor_merge=`perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result `;
468 print "perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result\n\n";
469 }
470
471 sub infor_merge_no_dec{
472 my $infor_merge_no_dec=`cp $annotate_dir\/sample_c_p.anno $dir\/total.result`;
473 }
474
475 sub genome_length{
476 my $length=`perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length`;
477 print "perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length\n\n"
478
479 }
480
481 sub plot{
482 $plot_dir="$dir\/plot\/";
483 mkdir ("$plot_dir");
484 my $span=defined($options{span})?$options{span}:50000;
485 my $cen="";
486 if (defined $options{cen}) {
487 $cen="-cen $options{cen}";
488 }
489 my $plot=`perl $path/sRNA_plot.pl -c $rpkm -g $dir/ref/genelist.txt -span 50000 -mark $sample_mark -l $dir/ref/genome\.length $cen -o $plot_dir/cluster.html -out $plot_dir/cluster.txt `;
490 "print perl $path/sRNA_plot.pl -c $rpkm -g $dir/ref/genelist.txt -span 50000 -mark $sample_mark -l $dir/ref/genome.length $cen -o $plot_dir/cluster.html -out $plot_dir/cluster.txt \n";
491
492 }
493
494 sub html{
495 my $pathfile="$dir/path.txt";
496 open PA,">$pathfile";
497 print PA "$config\n";
498 print PA "$preprocess\n";
499 print PA "$dir"."rfam_match\n";
500 print PA "$dir"."genome_match\n";
501 print PA "$cluster_file\n";
502 print PA "$annotate_dir\n";
503 print PA "$plot_dir\n";
504 if (defined($deg_dir)) {
505 print PA "$deg_dir\n";
506 }
507 close PA;
508 my $html=`perl $path\/html.pl -i $pathfile -format $format -o $dir/result.html`;
509 }
510
511 sub Time{
512 my $time=time();
513 my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
514 $month++;
515 $year+=1900;
516 if (length($sec) == 1) {$sec = "0"."$sec";}
517 if (length($min) == 1) {$min = "0"."$min";}
518 if (length($hour) == 1) {$hour = "0"."$hour";}
519 if (length($day) == 1) {$day = "0"."$day";}
520 if (length($month) == 1) {$month = "0"."$month";}
521 print "$year-$month-$day $hour:$min:$sec\n";
522 return("$year-$month-$day-$hour-$min-$sec");
523 }
524 #################################################################################