Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Agat sp compare two annotations 435 #437

Merged
merged 4 commits into from
Mar 7, 2024
Merged

Agat sp compare two annotations 435 #437

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
8c67c49
fix #435 - output list of IDs in files. Result is now a folder not an…
Juke34 Mar 5, 2024
f928bb4
update doc
Juke34 Mar 5, 2024
1259216
update test to reflect it is now an output folder
Juke34 Mar 5, 2024
3f3512c
manage file handler properly for ID output
Juke34 Mar 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 85 additions & 16 deletions bin/agat_sp_compare_two_annotations.pl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,9 +54,22 @@
# --- Manage config ---
$config = get_agat_config({config_file_in => $config});

######################
# Manage output folder #

if (! $opt_output) {
print "Default output name: split_result\n";
$opt_output="comparison_result";
}

if (-d $opt_output){
print "The output directory choosen already exists. Please give me another Name.\n";exit();
}
mkdir $opt_output;

######################
# Manage output file #
my $report = prepare_fileout($opt_output);
my $report = prepare_fileout("$opt_output/report.txt");

# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> MAIN <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

Expand Down Expand Up @@ -323,7 +336,9 @@
my $flat_overlap_1 = $locations1;
my $flat_overlap_2 = $locations2;
my $overlap_A=1;
my @overlap_A_id=($locations1->[0][3]);
my $overlap_B=1;
my @overlap_B_id=($locations2->[0][3]);
my $loop="top";
my $flip=2;
my $current_locs;
Expand Down Expand Up @@ -372,9 +387,11 @@

if($loop eq "top"){
$overlap_A++;
push @overlap_A_id, $locations->[0][3];
}
elsif($loop eq "bot"){
$overlap_B++;
push @overlap_B_id, $locations->[0][3];
}
}
else{
Expand All @@ -398,14 +415,14 @@
}
}
}
$overlap_info{$locations1->[0][2]}{$overlap_A}{$overlap_B}++;
push @{$overlap_info{$locations1->[0][2]}{$overlap_A}{$overlap_B}}, [ [@overlap_A_id], [@overlap_B_id] ];
}
# ----------------------------------- CASE 1 -----------------------------------
# location A ----------------
# location B ---------------
elsif ($locations2->[scalar(@{$locations2})-1][1] < $locations1->[0][0] ){

$overlap_info{$locations1->[0][2]}{0}{1}++;
push @{$overlap_info{$locations1->[0][2]}{0}{1}}, [[undef], [$locations2->[0][3]]];
print "Case1 notoverlap !\n\n" if ($debug);
# throw loc2
remove_loc_by_id($flattened_locations2_clean_sorted, $locusID, $chimere_type, $locations2->[0][3]);
Expand All @@ -420,7 +437,7 @@
my $id1 = $locations1->[0][3];
print "Case2 notoverlap !\n" if ($debug);
if(! exists_keys(\%seen1, ( $id1 ) ) ){ # else it has been dealed by overlap case
$overlap_info{$locations1->[0][2]}{1}{0}++;
push @{$overlap_info{$locations1->[0][2]}{1}{0}}, [ [$locations1->[0][3]], [undef] ];
# throw loc1
remove_loc_by_id($flattened_locations1_clean_sorted, $locusID, $chimere_type, $id1);
$seen1{$id1}++;
Expand All @@ -442,7 +459,7 @@
foreach my $locusID ( keys %{$flattened_locations1_clean_sorted} ){
foreach my $chimere_type ( keys %{$flattened_locations1_clean_sorted->{$locusID}}){
foreach my $locations1 ( @{$flattened_locations1_clean_sorted->{$locusID}{$chimere_type}} ){
$overlap_info{$locations1->[0][2]}{1}{0}++;
push @{$overlap_info{$locations1->[0][2]}{1}{0}}, [ [$locations1->[0][3]], [undef]];
print " Case3 !\n" if ($debug);
}
}
Expand All @@ -454,7 +471,7 @@
foreach my $locusID ( keys %{$flattened_locations2_clean_sorted} ){
foreach my $chimere_type ( keys %{$flattened_locations2_clean_sorted->{$locusID}}){
foreach my $locations2 ( @{$flattened_locations2_clean_sorted->{$locusID}{$chimere_type}} ){
$overlap_info{$locations2->[0][2]}{0}{1}++;
push @{$overlap_info{$locations2->[0][2]}{0}{1}}, [ [undef], [$locations2->[0][3]] ];
print " Case4 !\n" if ($debug);
}
}
Expand All @@ -481,22 +498,71 @@

my $string_to_print = "usage: $0 @copyARGV\nResults of number of genes from file1 that overlap genes from file2:\n\n";


my %file_handler;
foreach my $type_l1 ( sort keys %overlap_info ){

$string_to_print .= "$separator_table|".sizedPrint("$type_l1",92)."|\n";
$string_to_print .= "$separator_table|".sizedPrint($filename1.$ext1,30)."|".sizedPrint($filename2.$ext2,30)."|".sizedPrint("Number of cases",30)."|\n$separator_table";

$total{$type_l1}{'A'}=0;
$total{$type_l1}{'B'}=0;
# Create file handlers

foreach my $value1 ( sort {$a <=> $b} keys %{$overlap_info{$type_l1}} ){
foreach my $value2 ( sort {$a <=> $b} keys %{$overlap_info{$type_l1}{$value1}} ){
# report ids file_handler
my $report_ids = prepare_fileout("$opt_output/$type_l1"."_"."$value1"."_".$value2."_id_list.txt");
$file_handler{$type_l1."_".$value1."_".$value2}=$report_ids;
}
}

$total{$type_l1}{'A'}=0;
$total{$type_l1}{'B'}=0;
foreach my $value1 ( sort {$a <=> $b} keys %{$overlap_info{$type_l1}} ){
foreach my $value2 ( sort {$a <=> $b} keys %{$overlap_info{$type_l1}{$value1}} ){
$string_to_print .= "|".sizedPrint($value1,30)."|".sizedPrint($value2,30)."|".sizedPrint($overlap_info{$type_l1}{$value1}{$value2},30)."|\n";
if ($value1 != 0){
$total{$type_l1}{'A'} += $value1 * $overlap_info{$type_l1}{$value1}{$value2};
$string_to_print .= "|".sizedPrint($value1,30)."|".sizedPrint($value2,30)."|".sizedPrint(scalar(@{$overlap_info{$type_l1}{$value1}{$value2}}),30)."|\n";

# Use proper file handler for outputing IDs
my $report_ids = $file_handler{$type_l1."_".$value1."_".$value2};

# file1
foreach my $array ( @{$overlap_info{$type_l1}{$value1}{$value2}}) {
my $cpt=0;
my $last = scalar(@{$array->[0]});
foreach my $value ( @{$array->[0]} ) { # array0 is id overlarpA
$cpt++;
if(! $value){
print $report_ids "-";
} else {
if ($last == $cpt){
print $report_ids $value;
} else {
print $report_ids $value.", ";
}
}
}
print $report_ids " | ";
my $cpt2=0;
my $last2 = scalar(@{$array->[1]});
foreach my $value ( @{$array->[1]} ) { # array1 is id overlarpB
$cpt2++;
if(! $value){
print $report_ids "-\n";
} else {
if ($last2 == $cpt2){
print $report_ids "$value"
} else {
print $report_ids "$value, "
}
}
}
print $report_ids "\n";
}


if ($value1 != 0){
$total{$type_l1}{'A'} += $value1 * scalar(@{$overlap_info{$type_l1}{$value1}{$value2}});
}
if ($value2 != 0){
$total{$type_l1}{'B'} += $value2 * $overlap_info{$type_l1}{$value1}{$value2};
$total{$type_l1}{'B'} += $value2 * scalar(@{$overlap_info{$type_l1}{$value1}{$value2}});
}
}
}
Expand Down Expand Up @@ -648,7 +714,7 @@ =head1 DESCRIPTION

The script aims to compare two annotation of the same assembly. It provided
information about split/fusion of genes between the two annotations.
The most common case are:
The most common cases are:
1 => 0 ( gene uniq to file1)
0 => 1 ( gene uniq to file2)
1 => 1 ( 1 gene from file 1 overlaps only 1 gene from file2)
Expand All @@ -658,6 +724,9 @@ =head1 DESCRIPTION
Then you can get more complex cases:
<many> => <many> (<many> genes from file 1 overlap <many> genes from file2)

The script output a folder containing a report of number of different cases as well as a file
per case type listing per line the gene feature's ID involved in each case.

=head1 SYNOPSIS

agat_sp_compare_two_annotations.pl -gff1 infile.gff [ -o outfile ]
Expand All @@ -677,8 +746,8 @@ =head1 OPTIONS

=item B<-o> , B<--output> or B<--out>

Output GFF file. If no output file is specified, the output will be
written to STDOUT.
Output folder. It contains a report that resume the type and number of cases, as well as a file per case type
containing one case per line with the list of gene feature's ID (or other type of feature level1) from file1 then file2 separated by a |.

=item B<--debug> or B<-d>

Expand Down
7 changes: 5 additions & 2 deletions docs/tools/agat_sp_compare_two_annotations.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,9 @@ The most common case are:
Then you can get more complex cases:
&lt;many> => &lt;many> (&lt;many> genes from file 1 overlap &lt;many> genes from file2)

The script output a folder containing a report of number of different cases as well as a file
per case type listing per line the gene feature's ID involved in each case.

## SYNOPSIS

```
Expand All @@ -33,8 +36,8 @@ agat_sp_compare_two_annotations.pl --help

- **-o** , **--output** , **--out** or **--outfile**

Output GFF file. If no output file is specified, the output will be
written to STDOUT.
Output folder. It contains a report that resume the type and number of cases, as well as a file per case type
containing one case per line with the list of gene feature's ID (or other type of feature level1) from file1 then file2 separated by a |.

- **-v**

Expand Down
12 changes: 6 additions & 6 deletions t/scripts_output.t
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -238,22 +238,22 @@ $script = $script_prefix."bin/agat_sp_compare_two_annotations.pl";
$result = "$output_folder/agat_sp_compare_two_annotations_1.txt";
system(" $script --gff1 $input_folder/1.gff --gff2 $input_folder/1.gff -o $outtmp 2>&1 1>/dev/null");
#run test
ok( system("diff $result $outtmp") == 0, "output $script");
unlink $outtmp;
ok( system("diff $result $outtmp/report.txt") == 0, "output $script");
rmtree $outtmp;

$script = $script_prefix."bin/agat_sp_compare_two_annotations.pl";
$result = "$output_folder/agat_sp_compare_two_annotations_2.txt";
system(" $script --gff1 $input_folder/agat_sp_compare_two_annotations/file1.gff --gff2 $input_folder/agat_sp_compare_two_annotations/file2.gff -o $outtmp 2>&1 1>/dev/null");
#run test
ok( system("diff $result $outtmp") == 0, "output $script");
unlink $outtmp;
ok( system("diff $result $outtmp/report.txt") == 0, "output $script");
rmtree $outtmp;

$script = $script_prefix."bin/agat_sp_compare_two_annotations.pl";
$result = "$output_folder/agat_sp_compare_two_annotations_3.txt";
system(" $script --gff1 $input_folder/agat_sp_compare_two_annotations/file2.gff --gff2 $input_folder/agat_sp_compare_two_annotations/file1.gff -o $outtmp 2>&1 1>/dev/null");
#run test
ok( system("diff $result $outtmp") == 0, "output $script");
unlink $outtmp;
ok( system("diff $result $outtmp/report.txt") == 0, "output $script");
rmtree $outtmp;

# --------check agat_sp_compare_two_BUSCOs.pl -------------

Expand Down
Loading