Mercurial > repos > big-tiandm > sirna_plant

comparison siRNA.pl @ 19:e0884a4b996b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author big-tiandm
date Wed, 05 Nov 2014 01:17:26 -0500
parents 22d79320085c
children
comparison
equal deleted inserted replaced
18:22d79320085c 19:e0884a4b996b
137 my $rpkm=$dir."cluster\/"."sample_rpkm\.cluster"; 137 my $rpkm=$dir."cluster\/"."sample_rpkm\.cluster";
138 my $preprocess; 138 my $preprocess;
139 my $cluster_file; 139 my $cluster_file;
140 my $annotate_dir; 140 my $annotate_dir;
141 my $deg_dir; 141 my $deg_dir;
142 my $plot_dir;
142 my %id; 143 my %id;
143 for (my $i=0;$i<@mark ;$i++) { 144 for (my $i=0;$i<@mark ;$i++) {
144 $id{$mark[$i]}=$i+4; 145 $id{$mark[$i]}=$i+4;
145 } 146 }
146 147
147 148 print "\n######## tiandm test start ###########\n";
149 print "\@mark: @mark\n\%id keys number:";
150 print scalar keys %id;
151 print "\n";
152 foreach my $kyess (keys %id){
153 print $kyess," --> $id{$kyess}\n";
154 }
155 print "\n######## tiandm test end ############\n";
148 group_and_filter(); #collapse reads to tags 156 group_and_filter(); #collapse reads to tags
149 157
150 rfam(); 158 rfam();
151 159
152 my @map_read; 160 my @map_read;
278 my $l_out=$dir."preProcess\/"."collapse_reads_18-40.fa"; 286 my $l_out=$dir."preProcess\/"."collapse_reads_18-40.fa";
279 my $length_f=`perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark`; 287 my $length_f=`perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark`;
280 print "perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark\n\n"; 288 print "perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark\n\n";
281 my $cout_f=`perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark`; 289 my $cout_f=`perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark`;
282 print "perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark\n\n"; 290 print "perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark\n\n";
291 my $plot_l_D=`perl $path/Length_Distibution.pl -i $dir/preProcess/reads_length_distribution_after_count_filter.txt -o $dir/preProcess/length.html `;
292 print "perl $path\/Length_Distibution.pl -i $dir\/preProcess\/reads_length_distribution_after_count_filter.txt -o $dir\/preProcess\/length\.html\n\n";
283 return 0; 293 return 0;
284 } 294 }
285 295
286 sub rfam{ 296 sub rfam{
287 if (defined $options{'idx2'}) { 297 if (defined $options{'idx2'}) {
339 if ($cluster_mothod eq "F") { 349 if ($cluster_mothod eq "F") {
340 my $cluster=`perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt`; 350 my $cluster=`perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt`;
341 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"; 351 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";
342 } 352 }
343 elsif($cluster_mothod eq "T"){ 353 elsif($cluster_mothod eq "T"){
344 my $cluster=`perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $read_txt -mark $sample_mark`; 354 my $cluster=`perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $read_txt -k $sample_mark`;
345 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"; 355 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";
346 } 356 }
347 else{print "\-p is wrong!\n\n";} 357 else{print "\-p is wrong!\n\n";}
348 return 0; 358 return 0;
349 } 359 }
350 360
464 print "perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length\n\n" 474 print "perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length\n\n"
465 475
466 } 476 }
467 477
468 sub plot{ 478 sub plot{
469 my $plot_file="$dir\/plot\/"; 479 $plot_dir="$dir\/plot\/";
470 mkdir ("$plot_file"); 480 mkdir ("$plot_dir");
471 my $genome_plot="$dir\/plot\/genome\/";
472 mkdir ("$genome_plot");
473 #genome cluster
474 my $span=defined($options{span})?$options{span}:50000; 481 my $span=defined($options{span})?$options{span}:50000;
475 foreach (1..$sample_number) { 482 my $cen="";
476 my $mark=$mark[$_-1]; 483 if (defined $options{cen}) {
477 my $cen=""; 484 $cen="-cen $options{cen}";
478 if (defined $options{cen}) { 485 }
479 $cen="-cen $options{cen}"; 486 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 `;
480 } 487 "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";
481 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`; 488
482 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";
483 }
484
485 my $chr_plot_dir="$dir\/plot\/chr\/";
486 mkdir("$chr_plot_dir");
487 my %chr;
488 open LEN,"<$dir\/ref\/genome\.length";
489 while (my $aline=<LEN>) {
490 next if($aline=~/^\#/);
491 chomp $aline;
492 my @temp=split/\t/,$aline;
493 $chr{$temp[0]}=$temp[1];
494 }
495 close LEN;
496 foreach my $chr (sort keys %chr) {
497 my $cen="";
498 if (defined $options{cen}) {
499 $cen="-cen $options{cen}";
500 }
501 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`;
502 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";
503 }
504 } 489 }
505 490
506 sub html{ 491 sub html{
507 my $pathfile="$dir/path.txt"; 492 my $pathfile="$dir/path.txt";
508 open PA,">$pathfile"; 493 open PA,">$pathfile";
510 print PA "$preprocess\n"; 495 print PA "$preprocess\n";
511 print PA "$dir"."rfam_match\n"; 496 print PA "$dir"."rfam_match\n";
512 print PA "$dir"."genome_match\n"; 497 print PA "$dir"."genome_match\n";
513 print PA "$cluster_file\n"; 498 print PA "$cluster_file\n";
514 print PA "$annotate_dir\n"; 499 print PA "$annotate_dir\n";
500 print PA "$plot_dir\n";
515 if (defined($deg_dir)) { 501 if (defined($deg_dir)) {
516 print PA "$deg_dir\n"; 502 print PA "$deg_dir\n";
517 } 503 }
518 close PA; 504 close PA;
519 my $html=`perl $path\/html.pl -i $pathfile -format $format -o $dir/result.html`; 505 my $html=`perl $path\/html.pl -i $pathfile -format $format -o $dir/result.html`;