scripts/coding_potential_track_faster.pl at main · cnag-aat/scripts · GitHub

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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#!/usr/bin/env perl
use strict;
use Getopt::Long;
# conda activate /software/assembly/conda/perl-bioperl-1.7.8
#use lib "/project/devel/aateam/perlmods";
use lib "/scratch_isilon/groups/assembly/talioto/myperlmods";
use Bio::Seq;
use Bio::SeqIO;
use Bio::DB::Fasta;
use File::Basename qw( fileparse );
use SeqOp;

# Need to add GFF2 FITS (geneid-like) output
my $fa= 0;
my $verbose = 0;
my $geneid_param = "~/repos/geneid/param/human.101007.scoring.param";
my $frame = 1;
my $strand = "+";
my $window = 60;
GetOptions(
           'seq:s'          => \$fa,
	   'param:s'        => \$geneid_param,
	   'frame|f:s'       => \$frame,
	   'strand:s'       => \$strand,
	   'w:s'            => \$window
	   );
if ($strand =~/(-|r|R|minus)/){$strand = "-";}else{$strand = "+";}
# open the first file: table of codon usage frequencies

open(PARAM,"<$geneid_param") or die "$0: the file $geneid_param can not be opened: $!\n";
if($frame!~/[123]/){die "-frame must be 1, 2, or 3";}
# load the frequencies of codon usage into the hash table %pcodons
# this hash is indexed by a triplet of nucleotides or codon

my %llhex;
while (my $line = <PARAM>) {
  next if $line!~/Markov_Transition_probability_matrix/;
  while(my $entry = <PARAM>){
    last if $entry!~/\w/;
    last if $entry=~/^#/;
    chomp $entry;
    my($hex,$i,$frame,$lls)=split /\s+/,$entry;
    $llhex{$hex}{$frame}=$lls;
  }
}
close(PARAM);

my $db = Bio::DB::Fasta->new($fa);

my @seqids = sort $db->ids;
my $start=int($window/2)+($frame-1)+1;

foreach my $id (@seqids){
  my $chrlen =$db->length($id);
  my $rem = $chrlen % 3;
  my $sobj = $db->get_Seq_by_id($id);
  my $seq = $sobj->seq;
  print "fixedStep  chrom=$id  start=$start  step=3\n";
  if($strand eq '+'){
    for (my $base = $frame;$base<($chrlen-$window+1);$base+=3){
      #my $dna = SeqOp::get_seq_BioDBFasta($db,$id,$base,$base+$window-1,'+');
      #my $dna = $sobj->subseq($base,$base+$window-1);
      my $dna = substr($seq,$base-1,$window);
      print compute_cp($dna,0),"\n";
    }
  }else{
    for (my $base = $frame;$base<($chrlen-$window+1);$base+=3){
      #my $dna = SeqOp::get_seq_BioDBFasta($db,$id,$base,$base+$window-1,'+');
      #my $dna = $sobj->subseq($base,$base+$window-1);
      #my $rdna = SeqOp::get_seq_BioDBFasta($db,$id,$base,$base+$window-1,'-');
      my $dna = substr($seq,$base-1,$window);
      my $rdna = reverse($dna);
      $rdna =~ tr/ACGTacgt/TGCATGCA/;
      print compute_cp($rdna,0),"\n";
    }
  }

}
sub compute_cp {
  my $seq   = shift;
  my $ucseq = uc($seq);
  #print STDERR $ucseq,"\n";
  my $frame = shift;
  my $f = $frame;
  my $score = 0;
  for(my $p=0;$p<(length($ucseq)-5);$p++){
    my $h=substr($ucseq,$p,6);
    #print STDERR "$h\t$f\t$llhex{$h}{$f}\n";
    $score+=$llhex{$h}{$f};
    $f++;
    if($f==3){$f=0;}
  }
  return $score;
}