diff bin/piPipe.pl @ 0:198009598544 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/bin/piPipe.pl	Wed Oct 11 09:57:58 2017 -0400
@@ -0,0 +1,269 @@
+#!/usr/bin/perl
+use strict;
+use warnings;
+use Getopt::Long;
+use Parallel::ForkManager;
+use File::Basename;
+use File::Copy::Recursive qw( dircopy );
+use POSIX;
+use FindBin;
+use lib $FindBin::Bin;
+use resize qw ( size_distribution );
+use subgroups qw (subgroups );
+use ppp qw ( ping_pong_partners );
+use Rcall qw (pie_chart bg_to_png );
+use align qw ( to_build get_unique sam_count sam_count_mis sam_sorted_bam rpms_rpkm BWA_call get_fastq_seq extract_sam sam_to_bam_bg );
+use html qw ( main_page details_pages menu_page ppp_page );
+use File::Copy;
+
+my ( @fastq, @fastq_n, $dir, $min, $max, $mis, $misTE, $help, $Pcheck, $Dcheck, $mapnumf, $html_out);
+my ( $ref, $tRNAs, $rRNAs, $snRNAs, $miRNAs, $exons, $TE );
+my ( $si_min, $si_max, $pi_min, $pi_max );
+my ( $build_index, $build_tRNAs, $build_rRNAs, $build_snRNAs, $build_miRNAs, $build_exons, $build_TE );
+my $max_procs = 4;
+
+( $build_index, $build_tRNAs, $build_rRNAs, $build_snRNAs, $build_miRNAs, $build_exons, $build_TE ) = (0,0,0,0,0,0,0);
+( $min, $max, $mis, $misTE, $si_min, $si_max, $pi_min, $pi_max, $dir ) = ( 18, 29, 0, 3, 21, 21, 23, 29 );
+$Pcheck ='true';
+
+GetOptions (
+	"fastq=s" => \@fastq,
+	"fastq_n=s" => \@fastq_n,
+	"dir=s" => \$dir,
+	"min:i" => \$min,
+	"max:i" => \$max,
+	"si_min:i" => \$si_min,
+	"si_max:i" => \$si_max,
+	"pi_min:i" => \$pi_min,
+	"pi_max:i" => \$pi_max,
+	"mis:i" => \$mis,
+	"misTE:i" => \$misTE,
+	"html:s" => \$html_out,
+	"PPPon:s" => \$Pcheck,
+	"Dison:s" => \$Dcheck,
+	"help"   =>  \$help,
+	"ref:s" => \$ref,
+	"tRNAs:s" => \$tRNAs,
+	"rRNAs:s" => \$rRNAs,
+	"snRNAs:s" => \$snRNAs,
+	"miRNAs:s" => \$miRNAs,
+	"exons:s" => \$exons,
+	"TE:s" => \$TE,
+	"build_index" => \$build_index,
+	"build_tRNAs" => \$build_tRNAs,
+	"build_snRNAs" => \$build_snRNAs,
+	"build_miRNAs" => \$build_miRNAs,
+	"build_exons" => \$build_exons,
+	"build_rRNAs" => \$build_rRNAs,
+	"build_TE" => \$build_TE
+);
+
+my $fq_collection = 'fastq_dir/';
+mkdir $dir; mkdir $fq_collection;
+$dir = $dir.'/' unless $dir =~ /\/$/;
+mkdir $dir.'/css';mkdir $dir.'/js';
+dircopy( $FindBin::Bin.'/css', $dir.'/css' );
+dircopy( $FindBin::Bin.'/js', $dir.'/js' );
+
+my $file = $dir.'report.txt';
+open my $report, '>', $file or die "Cannot open $file $!\n";
+
+my %toBuild;
+@toBuild{ ( $ref, $tRNAs, $rRNAs, $snRNAs, $miRNAs, $exons, $TE ) } =  ( $build_index, $build_tRNAs, $build_rRNAs, $build_snRNAs, $build_miRNAs, $build_exons, $build_TE ) ;
+to_build ( \%toBuild, $report );
+
+my $proc_child = ceil($max_procs / scalar(@fastq));
+my $proc_grand_child = ceil($proc_child/5);
+my $pm = Parallel::ForkManager->new($max_procs);
+my $pm3 = Parallel::ForkManager->new($proc_grand_child);
+
+$pm->run_on_finish( sub {
+	my ($pid, $exit_code, $ident) = @_;
+	print $report "Fastq fork $ident just finished ".
+	"with PID $pid and exit code: $exit_code\n";
+	die "Something went wrong!\n" if $exit_code != 0;
+	});
+$pm->run_on_start( sub {
+	my ($pid,$ident)=@_;
+	print $report "Fastq fork : $ident started, pid: $pid\n";
+	});
+$pm3->run_on_finish( sub {
+	my ($pid, $exit_code, $ident) = @_;
+	print $report "** Subgroup fork $ident just finished ".
+	"with PID $pid and exit code: $exit_code\n";
+	die "Something went wrong!\n" if $exit_code != 0;
+	});
+$pm3->run_on_start( sub {
+	my ($pid,$ident)=@_;
+	print $report "** Subgroup fork $ident started, pid: $pid\n";
+	});
+
+
+foreach my $child ( 0 .. $#fastq )
+{
+	my @suffix = ('.fastq', '.fastq.gz,', '.fq', '.fq.gz', 'ref', '.dat', '.fa','.fas','.fasta', '.txt');
+	my ( $name, $path, $suffix ) = fileparse( $fastq[$child], @suffix );
+	my ( $ref_name, $ref_path, $ref_suffix ) = fileparse( $ref, @suffix );
+	my ( $TE_name, $TE_path, $TE_suffix ) = fileparse( $TE, @suffix );
+	my ( $ex_name, $ex_path, $ex_suffix ) = fileparse( $exons,@suffix );
+
+	$pm->start($fastq[$child]) and next;
+
+	my $dir_fq = $dir.$name.'/';
+	mkdir $dir_fq;
+
+	my $gen_dir = $dir_fq.'genome/';
+	mkdir $gen_dir;
+
+	my $size_dir = $dir_fq.'size/';
+	mkdir $size_dir;
+
+	my $fastq_resized = $dir_fq.$name.'_'.$min.'-'.$max.'.fastq';
+	size_distribution (  $fastq[$child], $fastq_resized, $size_dir, $min, $max );
+
+	my $sam_genome = $gen_dir.$name.'_'.$min.'-'.$max.'_'.$ref_name.'.sam';
+	my $sam_genome_unique = $gen_dir.$name.'_'.$min.'-'.$max.'_'.$ref_name.'_unique.sam';
+	my $fastq_prefix = $gen_dir.$name.'_'.$min.'-'.$max.'_'.$ref_name;
+
+	BWA_call ( $ref, $fastq_resized, $sam_genome, $mis, $proc_child, $report );
+	my ( $fai_ref_hashP, $ma, $ma_uni ) = get_unique ( $sam_genome, $sam_genome_unique, $gen_dir, 1, $report );
+
+	my $scale = 1000000 / $ma;
+	sam_to_bam_bg ( $sam_genome_unique, $scale, $proc_child );
+	sam_to_bam_bg ( $sam_genome, $scale, $proc_child );
+
+	my $Gviz_dir = $gen_dir.'Gviz/';
+	my $fai_file =  $gen_dir.'fai';
+	mkdir $Gviz_dir;
+	my $Gviz_dir_rand = $Gviz_dir.'rand/';
+	mkdir $Gviz_dir_rand;
+	my $Gviz_dir_uni = $Gviz_dir.'unique/';
+	mkdir $Gviz_dir_uni;
+
+	open my $gfai, '>', $fai_file;
+	foreach my $k  ( sort keys %{$fai_ref_hashP} )
+	{
+		print $gfai "$k\t$fai_ref_hashP->{$k}\n";
+	}
+	close $gfai;
+	bg_to_png ( $fai_file, $fastq_prefix.'_unique_plus.bedgraph', $fastq_prefix.'_unique_minus.bedgraph', $Gviz_dir_uni, 'Mb' );
+	bg_to_png ( $fai_file, $fastq_prefix.'_plus.bedgraph', $fastq_prefix.'_minus.bedgraph', $Gviz_dir_rand, 'Mb' );
+
+	my $group_dir = $dir_fq.'subgroups/';
+	my $fastq_uni = $gen_dir.'unique.fastq';
+	my $fastq_all = $gen_dir.'all.fastq';
+	my ($bo, $mi, $pi) = subgroups ( $fastq_all, $group_dir, $mis, $misTE, $proc_child, $tRNAs, $rRNAs, $snRNAs, $miRNAs, $exons, $TE, $si_min, $si_max, $pi_min, $pi_max, $report);
+
+	pie_chart($group_dir);
+
+	open (my $dupnum, $gen_dir.'dup_mapnum.txt') || die "cannot open dup_mapnum.txt $!";
+	my %dupnum_genome;
+	my $header = <$dupnum>;
+	while (<$dupnum>)
+	{
+		chomp $_;
+		my @dupline = split /\t/, $_;
+		$dupnum_genome{$dupline[0]} = [$dupline[1], $dupline[2]];
+	}
+	close $dupnum;
+
+	my $mi_sam = $group_dir.'miRNAs.sam';
+	mkdir $group_dir.'miRNAs/';
+	my $mi_count_file =  $group_dir.'miRNAs/miRNAs_reads_counts.txt';
+	my ( $mi_count, $mi_ref_size ) = sam_count ( $mi_sam );
+
+	rpms_rpkm( $mi_count, $mi_ref_size, $ma, $mi_count_file, $pi, $mi, $bo );
+
+	my (  $sam_exons, $sam_TEs ) = ( $group_dir.'exons.sam', $group_dir.'TEs.sam' );
+	my @types = ($group_dir.'bonafide_reads.fastq', $group_dir.'miRNAs.fastq', $group_dir.'siRNAs.fastq', $group_dir.'piRNAs.fastq' );
+	my @types_names = ('bonafide_reads', 'miRNAs', 'siRNAs', 'piRNAs');
+	foreach my $grand_child ( 0 .. $#types )
+	{
+		my $type_dir = $group_dir.$types_names[$grand_child].'/';
+		my $type_prefix = $types_names[$grand_child].'-';
+		mkdir  $type_dir;
+
+		my ( $type_sam_genome, $type_sam_TEs, $type_sam_exons ) = ( $type_dir.$type_prefix.'genome.sam', $type_dir.$type_prefix.'TEs.sam', $type_dir.$type_prefix.'exons.sam' );
+		my ( $type_sam_uni_genome, $type_sam_uni_TEs,  $type_sam_uni_exons ) = ( $type_dir.$type_prefix.'genome_unique.sam', $type_dir.$type_prefix.'TEs_unique.sam', $type_dir.$type_prefix.'exons_unique.sam' );
+		my ( $type_uni_genome_fastq, $type_uni_TEs_fastq,  $type_uni_exons_fastq ) = ( $fq_collection.$type_prefix.'genome_uni.fastq', $fq_collection.$type_prefix.'TEs_uni.fastq', $fq_collection.$type_prefix.'exons_uni.fastq');
+		my ( $type_genome_fastq, $type_TEs_fastq,  $type_exons_fastq ) = ( $fq_collection.$type_prefix.'genome.fastq', $fq_collection.$type_prefix.'TEs.fastq', $fq_collection.$type_prefix.'exons.fastq');
+		my $type_sequence_hashP = get_fastq_seq ( $types[$grand_child] );
+
+		if ( $grand_child == 1 )
+		{
+			BWA_call ( $TE, $types[$grand_child],  $type_sam_TEs, $misTE, $proc_child, $report );
+			BWA_call ( $exons, $types[$grand_child], $type_sam_exons, $mis, $proc_child, $report );
+		 	BWA_call ( $ref, $types[$grand_child], $type_sam_genome, $mis, $proc_child, $report );
+			extract_sam ( undef, $type_sam_TEs, $type_sam_TEs, $type_sam_uni_TEs, $type_uni_TEs_fastq, $type_uni_TEs_fastq );
+			extract_sam ( undef, $type_sam_exons, $type_sam_exons, $type_sam_uni_exons, $type_exons_fastq, $type_uni_exons_fastq );
+			extract_sam ( undef, $type_sam_genome, $type_sam_genome, $type_sam_uni_genome, $type_genome_fastq, $type_uni_genome_fastq );
+		}
+		else
+		{
+			extract_sam ( $type_sequence_hashP, $sam_TEs, $type_sam_TEs, $type_sam_uni_TEs, $type_TEs_fastq, $type_uni_TEs_fastq );
+			extract_sam ( $type_sequence_hashP, $sam_exons, $type_sam_exons, $type_sam_uni_exons, $type_exons_fastq, $type_uni_exons_fastq );
+			extract_sam ( $type_sequence_hashP, $sam_genome, $type_sam_genome, $type_sam_uni_genome, $type_genome_fastq, $type_uni_genome_fastq );
+		}
+
+		my $ex_count_file =  $type_dir.'exons_reads_counts.txt';
+		my ( $ex_count, $ex_ref_size ) =  sam_count ( $type_sam_exons );
+		rpms_rpkm( $ex_count, $ex_ref_size, $ma, $ex_count_file, $pi, $mi, $bo );
+
+		my ( $TEs_count, $TEs_ref_size, $TEs_count_NoM, $TEs_count_M ) = sam_count_mis ( $type_sam_TEs );
+		my $TEs_count_file = $type_dir.$type_prefix.'TEs_reads_counts.txt';
+		my $TEs_count_file_M = $type_dir.$type_prefix.'TEs_reads_counts_mismatches.txt';
+		my $TEs_count_file_noM = $type_dir.$type_prefix.'TEs_reads_counts_nomismatches.txt';
+		rpms_rpkm( $TEs_count, $TEs_ref_size, $ma, $TEs_count_file, $pi, $mi, $bo );
+		rpms_rpkm( $TEs_count_NoM, $TEs_ref_size, $ma, $TEs_count_file_M, $pi, $mi, $bo );
+		rpms_rpkm( $TEs_count_M, $TEs_ref_size, $ma, $TEs_count_file_noM, $pi, $mi, $bo );
+
+		sam_to_bam_bg ( $type_sam_TEs, $scale, $grand_child );
+		sam_sorted_bam ( $type_sam_exons, $grand_child ); sam_sorted_bam ( $type_sam_uni_exons, $grand_child ); 
+		sam_sorted_bam ( $type_sam_uni_TEs, $grand_child ); 
+
+		my $Gviz_TEs =  $type_dir.'Gviz_TEs/';
+		mkdir $Gviz_TEs;
+		bg_to_png ( $group_dir.'TEs.fai', $type_dir.$type_prefix.'TEs_plus.bedgraph', $type_dir.$type_prefix.'TEs_minus.bedgraph', $Gviz_TEs, 'Kb' );
+
+		my $Gviz_genome=  $type_dir.'Gviz_genome/';
+		my $Gviz_genome_rand = $Gviz_genome.'rand/';
+		my $Gviz_genome_uni = $Gviz_genome.'unique/';
+		mkdir $Gviz_genome; mkdir $Gviz_genome_uni; mkdir $Gviz_genome_rand;
+
+		sam_to_bam_bg ( $type_sam_genome, $scale, $grand_child );
+		sam_to_bam_bg ( $type_sam_uni_genome, $scale, $grand_child );
+
+		bg_to_png ( $fai_file, $type_dir.$type_prefix.'genome_unique_plus.bedgraph', $type_dir.$type_prefix.'genome_unique_minus.bedgraph', $Gviz_genome_uni, 'Mb' );
+		bg_to_png ( $fai_file, $type_dir.$type_prefix.'genome_plus.bedgraph', $type_dir.$type_prefix.'genome_minus.bedgraph', $Gviz_genome_rand, 'Mb' );
+
+		#HTML Details
+		my $prefix_details_pages = $dir.$fastq_n[$child].'-'.$types_names[$grand_child];
+		details_pages ( $type_dir, $prefix_details_pages, \@fastq_n, $fastq_n[$child], $misTE, $dir );
+
+		$pm3->finish($grand_child);
+	}
+	$pm3->wait_all_children;
+
+	if ( $Pcheck eq 'true' )
+	{
+		my $ppp = $group_dir.'PPPartners/'; mkdir $ppp;
+		print $report "ping_pong_partners $group_dir/bonafide_reads/TEs.sam $ppp\n";
+		ping_pong_partners ( $group_dir.'TEs.fai', $group_dir.'bonafide_reads/bonafide_reads-TEs_sorted.bam', $ppp, $min );
+		my $ppp_page = $dir.$fastq_n[$child].'-bonafide_reads-PPP.html';
+		ppp_page ( $group_dir, $ppp_page, \@fastq_n, $fastq_n[$child], $ppp, $dir );
+	}
+
+	#HTML Main Webpage
+	my $index_page = $dir.$fastq_n[$child].'.html';
+  main_page ( $gen_dir, $index_page, \@fastq_n, $fastq_n[$child], $ma, $ma_uni, $dir );
+  copy ($index_page, $html_out) if $child == 0;
+	#HTML Menu
+	my $menu_page = $dir.$fastq_n[$child].'-sub.html';
+	menu_page ( $group_dir, $menu_page, \@fastq_n, $fastq_n[$child], $min, $max, $si_min, $si_max, $pi_min, $pi_max, $dir );
+
+	$pm->finish($child); # pass an exit code to finish
+}
+$pm->wait_all_children;
+
+print $report "Job done!\n";
+close $report;