-
Notifications
You must be signed in to change notification settings - Fork 0
/
39.homoeologTriadsWithKOs.pl
75 lines (57 loc) · 2.77 KB
/
39.homoeologTriadsWithKOs.pl
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
#!/usr/bin/perl
##########################################
my $script='39.homoeologTriadsWithKOs.pl';
# Paul Bailey 15.12.2015
##########################################
use strict;
use warnings;
use Getopt::Long;
my ($triadsFile, $geneList, $usage);
$usage="
-----------------------------------------------------
Usage: perl $script -f <triads file> -g <gene list>
Required options:
-f or --infile: homoeolog triads table file
-g or --genelist: gene list
-----------------------------------------------------
";
GetOptions(
'f|triadsfile:s' => \$triadsFile,
'g|genelist:s' => \$geneList,
'h|help:s' => die $usage
);
my $geneListCountr = 0;
my $triadCountr = 0;
my $threeGenesCountr = 0;
my $twoGenesCountr = 0;
my $oneGeneCountr = 0;
my $noGenesCountr = 0;
open GENE_LIST, $geneList or die "Error opening file $geneList: $!\n\n$usage";
my %geneListHash;
while(my $line = <GENE_LIST>) {
chomp $line;
$geneListHash{$line} = 0;
$geneListCountr++;
}
open TRIADS_FILE, $triadsFile or die "Error opening $triadsFile: $!\n\n$usage";
my $headr = <TRIADS_FILE>;
while(my $line = <TRIADS_FILE>) {
chomp $line;
my($tribeId, $n_total_members, $n_conf_members, $tribe_genome_configuration, $tribe_chromosome_configuration, $members_a, $members_b, $members_c) = split '\t', $line;
if( exists $geneListHash{$members_a} && exists $geneListHash{$members_b} && exists $geneListHash{$members_c} ) {$threeGenesCountr++}
elsif( (exists $geneListHash{$members_a} && !exists $geneListHash{$members_b} && exists $geneListHash{$members_c})) {$twoGenesCountr++}
elsif( (exists $geneListHash{$members_a} && exists $geneListHash{$members_b} && !exists $geneListHash{$members_c})) {$twoGenesCountr++}
elsif( (!exists $geneListHash{$members_a} && exists $geneListHash{$members_b} && exists $geneListHash{$members_c})) {$twoGenesCountr++}
elsif(( exists $geneListHash{$members_a} && !exists $geneListHash{$members_b} && !exists $geneListHash{$members_c}) ) {$oneGeneCountr++}
elsif(( !exists $geneListHash{$members_a} && exists $geneListHash{$members_b} && !exists $geneListHash{$members_c}) ) {$oneGeneCountr++}
elsif(( !exists $geneListHash{$members_a} && !exists $geneListHash{$members_b} && exists $geneListHash{$members_c}) ) {$oneGeneCountr++}
elsif(( !exists $geneListHash{$members_a} && !exists $geneListHash{$members_b} && !exists $geneListHash{$members_c}) ) {$noGenesCountr++}
$triadCountr++;
}
print "\n# genes in list: $geneListCountr";
print "\n# gene triads in set: $triadCountr";
print "\n# triads with all genes present: $threeGenesCountr";
print "\n# triads with two genes present: $twoGenesCountr";
print "\n# triads with one gene present: $oneGeneCountr";
print "\n# triads with no gene present: $noGenesCountr";
print "\n";