Mercurial > repos > dereeper > pangenome_explorer

comparison Perl/ConvertPanacotaMatrix.pl @ 3:e42d30da7a74 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Thu, 30 May 2024 11:52:25 +0000
parents
children
comparison
equal deleted inserted replaced
2:97e4e3e818b6 3:e42d30da7a74
1 #!/usr/bin/perl
2
3 use strict;
4
5 my $indir = $ARGV[0];
6 my $matrix = $ARGV[1];
7 my $out = $ARGV[2];
8 my $strain_names = $ARGV[3];
9
10 my %strains_of_gb;
11 open(F,$strain_names);
12 while(<F>){
13 my $line = $_;
14 $line =~s/\n//g;$line =~s/\r//g;
15 my ($gb,$strain) = split(/\t/,$line);
16 $strains_of_gb{$gb} = $strain;
17 }
18 close(F);
19
20 my @strains;
21 my %corr;
22 open(D,"ls $indir/*pep |");
23 while(<D>){
24 my $file = $_;
25 my $prot_num = 0;
26 my $strain;
27 if ($file =~/\/([^\/]*).pep/){
28 $strain = $1;
29 $strain = $strains_of_gb{$strain};
30 push(@strains,$strain);
31 }
32 open(F,"$file");
33 while(<F>){
34 if (/>(.*)/){
35 my $prot_id = $1;
36 $prot_num++;
37 my $new_id = "$strain"."_".$prot_num;
38 $corr{$new_id} = $prot_id;
39 }
40 }
41 close(F);
42 }
43 close(D);
44
45 my $cl_num = 0;
46 my $nb_strains = 1;
47 open(O,">$out");
48 open(U,">$out.upsetr.txt");
49 open(M,">$out.accessory_01matrix.txt");
50 open(F,$matrix);
51 print O "ClutserID";
52 print U "ClutserID";
53 print M "Gene";
54 my %hash_place_strains;
55 my $num_cell = 0;
56 foreach my $strain(@strains){
57 $num_cell++;
58 my @words = split(/_/,$strain);
59 my $genus = $words[0];
60 my $species = $words[1];
61 my $shortname = substr($genus,0,3) . "_". substr($species,0,2);
62 for (my $j = 2; $j <= $#words; $j++){
63 $shortname.="_".$words[$j];
64 }
65 $shortname = substr($shortname,0,25);
66 print O "\t".$strain;
67 print U "\t".$shortname;
68 print M "\t".$shortname;
69 $hash_place_strains{$strain} = $num_cell;
70 $nb_strains++;
71 }
72 print O "\n";
73 print U "\n";
74 print M "\n";
75 while(<F>){
76 $cl_num++;
77 my $line = $_;
78 $line =~s/\n//g;$line =~s/\r//g;
79 my @infos = split(/ /,$line);
80 my %cells;
81 for (my $i = 1; $i <= $#infos; $i++){
82 my $new_id = $infos[$i];
83 my $prot_id = $corr{$new_id};
84 my $strain;
85 if ($new_id =~/^(.*)_\d+$/){$strain=$1;}
86 my $num_cell = $hash_place_strains{$strain};
87 #print "$strain $num_cell $prot_id $new_id\n";
88 $cells{$strain}.= $prot_id.",";
89 }
90 print O $cl_num;
91 print U $cl_num;
92 my $concat_accessory = "";
93 foreach my $strain(@strains){
94 my $val;
95 if ($cells{$strain}){
96 $val = $cells{$strain};
97 chop($val);
98 }
99 else{
100 $val = "-";
101 }
102 if ($val =~/\w+/){
103 print U "\t1";
104 $concat_accessory .= "\t1";
105 }
106 else{
107 print U "\t0";
108 $concat_accessory .= "\t0";
109 }
110 my $concat = $val;
111 print O "\t".$concat;
112 }
113 if ($concat_accessory =~/0/){
114 print M $cl_num.$concat_accessory."\n";
115 }
116 print O "\n";
117 print U "\n";
118 }
119 close(F);
120 close(O);