EndHiC/hcluster_suitable_stop.pl at master · whc2/EndHiC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#!/usr/bin/perl

=head1 Name

hcluster_suitable_stop.pl  -- stop clustering at suitable loop

=head1 Description


=head1 Version

  Author: Fan Wei, fanw@genomics.org.cn
  Version: 1.0,  Date: 2006-12-6
  Note:

=head1 Usage

  hcluster_suitable_stop.pl <hcluster_result_file>
  --chr_num <int>    expected chromosome number
  --chr_len <int>    minimum chromosome length
  --verbose   output running progress information to screen
  --help      output help information to screen

=head1 Exmple


=cut

use strict;
use Getopt::Long;
use FindBin qw($Bin $Script);
use File::Basename qw(basename dirname);
use Data::Dumper;


##get options from command line into variables and set default values

my $chromosome_number;
my $chromosome_length;
my ($Verbose,$Help);
GetOptions(
	"chr_num:i"=>\$chromosome_number,
	"chr_len:i"=>\$chromosome_length,
	"verbose"=>\$Verbose,
	"help"=>\$Help
);
$chromosome_number ||= 18;
$chromosome_length ||= 10000000;
die `pod2text $0` if (@ARGV == 0 || $Help);

#print STDERR "Parameters: $chromosome_number\t$chromosome_length\t$anchor_rate\n";

my $hcluster_log_file = shift;
my @DATA;

$/ = "Clustering loop";
open IN, $hcluster_log_file || die "fail $hcluster_log_file";
<IN>;
while (<IN>) {
	chomp;
	s/Clustering finished\s+//;
	#print "\n--------------------------------------------\n";
	my @lines = split /\n+/;
	my $loop_line = shift @lines;
	my $loop_id = $1 if($loop_line =~ /^\s+(\d+)/);
	my $group_num_line = shift @lines;
	#print Dumper \@lines;
	#print $loop_id."\n";


	my @data;
	foreach  (@lines) {
		my @t = split /\s+/;
		my ($cluster_id, $contig_num,$cluster_len, $robustness, $contigs_str) = @t;
		push @data, [$cluster_id, $contig_num,$cluster_len, $robustness, $contigs_str];
	}
	push @DATA, \@data;

}
close IN;

@DATA = reverse @DATA;


#print Dumper \@DATA;


for (my $i = $chromosome_number - 1; $i < @DATA; $i++) {

	my $cluster_number = 0;

	my $p = $DATA[$i];
	foreach my $pp (@$p) {
		if($pp->[2] > $chromosome_length){
			$cluster_number ++;
		}
	}

	#print STDERR "$i\t$cluster_number\n";
	if ($cluster_number == $chromosome_number) {
		my $loop_id = $i + 1;
		print "#Cluster_id\tcontig_count\tcluster_length\trobustness\tincluded_contigs(no_order_and_orient) [Loop_id: $loop_id]\n";
		foreach my $pp (@$p) {
			my $line = join("\t",@$pp);
			print $line."\n";
		}
		last;
	}

}


####################################################
################### Sub Routines ###################
####################################################