Concordo com o que @WeijunZhou disse em seu comentário, você não precisa fazer todos esses arquivos temporários para fazer isso.
O script perl a seguir calculará as contagens para a classificação do seu Método 1 (filiais) e Método 2 (GPS), em uma passagem pelos dois arquivos.
Funciona mantendo uma lista ordenada (array) dos valores de phylop e GPS do arquivo duplicado, e então (para cada linha do arquivo único) calculando onde o valor de phylop e GPS teria ordenado em seus respectivos arrays ordenados .
#!/usr/bin/perl
use strict;
# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';
# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
my @phylop=();
my @GPS=();
open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
chomp;
next if (m/^chromosoom/);
my($chr,$start,$end,$phylop,$GPS) = split;
push @phylop, $phylop + 0; # add 0 to make sure we only ever store a number
push @GPS, $GPS + 0;
};
close(DUPE);
# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$a <=> $b} @phylop;
@GPS = sort {$a <=> $b} @GPS;
print "Method1\tMethod2\n";
# Now find out where the phylop and GPS value from each line of uniqfile
# would have ended up if we had sorted it into dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
next if (m/^chromosoom/);
chomp;
my $phylop_sort_line=1;
my $GPS_sort_line=1;
my($chr,$start,$end,$phylop,$GPS) = split;
for my $i (0..@phylop-1) {
$phylop_sort_line++ if ($phylop < $phylop[$i]);
$GPS_sort_line++ if ($GPS < $GPS[$i]);
};
#printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;
};
close(UNIQ);
quando executado com base nos dados de amostra fornecidos acima, a saída é:
$ ./counts-for-methods.pl
Method1 Method2
2 1
1 1
3 2
4 3
4 5
7 7
8 7
O script ignora completamente a linha de cabeçalho em ambos os arquivos, portanto, esses números de linha podem ser desativados em um caso o seu algoritmo atual os conte.
Ele também pressupõe que os valores do arquivo exclusivo sempre serão classificados imediatamente antes de um valor idêntico no arquivo duplicado. Se isso não for o que você deseja, altere as comparações <
no for my $i (0..@phylop)
loop para <=
.
Se você precisar dos valores do Método 1 e do Método 2 separadamente, poderá extraí-los facilmente com awk
. Ou o script perl
pode ser facilmente modificado para abrir dois arquivos de saída, um para cada método, e imprimir os respectivos valores para cada arquivo.
aqui está uma versão para lidar com 151 campos nas linhas de entrada. Eu não tenho esse arquivo de entrada, então eu testei com a versão '5 fields' que foi comentada no código. A saída foi a mesma que a versão acima.
#!/usr/bin/perl
use strict;
# get uniqfile and dupefile names from cmd line, with defaults
my $uniqfile = shift || 'myUniqueFile';
my $dupefile = shift || 'myDuplicationFile';
my @phylop=();
my @GPS=();
# Read in the dupefile and keep the phylops and GPS values.
# This could take a LOT of memory if dupefile is huge.
# Most modern systems should have no difficulty coping with even
# a multi-gigabyte dupefile.
open(DUPE,"<",$dupefile) || die "couldn't open '$dupefile': $!\n";
while(<DUPE>) {
chomp;
next if (m/^chromosoom/);
my @fields = split;
# 151 fields version:
push @phylop, $fields[42]+0;
push @GPS, $fields[150]+0;
# 5 fields version:
# push @phylop, $fields[3]+0;
# push @GPS, $fields[4]+0;
};
close(DUPE);
# Sort the @phylop and @GPS arrays, numerically descending
@phylop = sort {$b <=> $a} @phylop;
@GPS = sort {$b <=> $a} @GPS;
print "Method1\tMethod2\n";
# Now find out where the phylop and GPS from each line of uniqfile
# would have ended up if we had sorted it into the dupefile
open(UNIQ,"<",$uniqfile) || die "couldn't open '$uniqfile': $!\n";
while (<UNIQ>) {
next if (m/^chromosoom/);
chomp;
my $phylop_sort_line=1;
my $GPS_sort_line=1;
my @fields = split;
for my $i (0..@phylop-1) {
# 151 fields version:
$phylop_sort_line++ if ($fields[42] < $phylop[$i]);
$GPS_sort_line++ if ($fields[150] < $GPS[$i]);
# 5 fields version:
# $phylop_sort_line++ if ($fields[3] < $phylop[$i]);
# $GPS_sort_line++ if ($fields[4] < $GPS[$i]);
};
#printf "%i\t%i\t#%s\n", $phylop_sort_line, $GPS_sort_line, $_;
printf "%i\t%i\n", $phylop_sort_line, $GPS_sort_line;
};
close(UNIQ);