forked from schraderL/GAGA
-
Notifications
You must be signed in to change notification settings - Fork 0
/
06c_analyze_absrel.pl
59 lines (39 loc) · 1.75 KB
/
06c_analyze_absrel.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
#!/usr/bin/perl
use strict;
use warnings;
use File::Basename;
my $dirname = dirname(__FILE__);
# Quick way: grep ' p-value =' absrel/*output | less -S
# usage: perl analyze_absrel.pl
# INPUT
my $cdsdir = "/home/projects/ku_00039/people/joeviz/orthology_trees/all_orthologs/Cleaned_alignments/run_all_withgard/hmmcleaner_50_aln/";
my $absreldir = "/home/projects/ku_00039/people/joeviz/orthology_trees/all_orthologs/Hyphy_absrel/absrel_withgardpartitions_hmmclean50aln/absrel_allhog/";
my $outdir = "absrel_output_tables";
system ("mkdir -p $outdir");
my $pvalue = "0.001"; # p-value to filter significant hits
my $fdr = "0.01"; # FDR to filter significant hits
my ($line, $name);
# Reading each alignment file for which absrel was run
#my %protfasta;
system ("ls $cdsdir\/\*\.nonstop.aln > tmp_protlist.txt");
open(File, "<", "tmp_protlist.txt");
while(<File>){
chomp;
$line = $_;
next if ($line !~ /\S+/);
my $id = "";
if ($line =~ /.*\/(\S+)\.cds/){
$id = $1;
} else {
die "Can't find OG id in $line\n";
}
my $absrelout = "$absreldir\/$id\_ABSREL.json";
my $tableout = "$outdir\/$id\_absrel_wdnds.tsv";
my $tablesignout = "$outdir\/$id\_absrel_significant\_pval$pvalue\.tsv";
my $tablesignfdrout = "$outdir\/$id\_absrel_significant\_fdr$fdr\.tsv";
system ("python3 /home/projects/ku_00039/people/joeviz/orthology_trees/all_orthologs/significant_aBSREL_2_table.py $absrelout $tablesignout 0 $pvalue \n");
system ("python3 /home/projects/ku_00039/people/joeviz/orthology_trees/all_orthologs/significant_aBSREL_2_table.py $absrelout $tablesignfdrout $fdr 0\n");
system ("python3 /home/projects/ku_00039/people/joeviz/orthology_trees/all_orthologs/aBSREL_2_table_withdnds.py $absrelout $tableout \n");
}
close File;
system ("rm tmp_protlist.txt");