Analysis of the complement and molecular evolution of tRNA genes in cow

Step by step analysis of identifying the bovine tRNA repertoire (using the most recent version of the bovine genome) based on the procedure published in this article. This was my first first-authored paper (published as part of the bovine genome companion papers), and in retrospect there could have been many improvements, however the paper was done on what I knew back then.

1. Download bosTau5.fa.gz from http://hgdownload-test.cse.ucsc.edu/goldenPath/bosTau5/bigZips/ and gunzip

2. Download and install tRNAscan-SE from http://lowelab.ucsc.edu/software/tRNAscan-SE.tar.gz

I had problems installing tRNAscan-SE on Ubuntu 11.04 (64 bit) but found the solution here: http://happykhan.com/?p=545

3. Setup RepeatMasker (instructions can be followed here). I downloaded repeatmaskerlibraries-20110419.tar.gz and RepeatMasker-open-3-3-0.tar.gz.

4. ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.2.25+-x64-linux.tar.gz

5. Once tRNAscan-SE is up, perform the search on bosTau5.fa

./tRNAscan-SE -o bos_tau_5_trnascan_se.output -f bos_tau_5_trnascan_se_trna.output -m bos_tau_5_trnascan_se.log ~/genome/bosTau5.fa

tRNAscan-SE v.1.3 (March 2011) scan results (on host tan118-VirtualBox)
Started: Sat Sep 10 14:33:30 JST 2011

————————————————————
Sequence file(s) to search: /home/tan118/genome/bosTau5.fa
Search Mode: Eukaryotic
Results written to: bos_tau_5_trnascan_se.output
Output format: Tabular
Searching with: tRNAscan + EufindtRNA -> Cove
Covariance model: TRNA2-euk.cm
tRNAscan parameters: Strict
EufindtRNA parameters: Relaxed (Int Cutoff= -32.1)
tRNA secondary structure
predictions saved to: bos_tau_5_trnascan_se_trna.output
Search log saved in:
Search statistics saved in: bos_tau_5_trnascan_se.log
————————————————————

First-pass (tRNAscan/EufindtRNA) Stats:
—————
Sequences read: 11894
Seqs w/at least 1 hit: 6798
Bases read: 2918090798 (x2 for both strands)
Bases in tRNAs: 39937512
tRNAs predicted: 519405
Av. tRNA length: 76
Script CPU time: 6153.84 s
Scan CPU time: 1249.34 s
Scan speed: 4671.4 Kbp/sec

First pass search(es) ended: Sat Sep 10 16:41:51 JST 2011

Cove Stats:
———–
Candidate tRNAs read: 519405
Cove-confirmed tRNAs: 225682
Bases scanned by Cove: 48247536
% seq scanned by Cove: 0.8 %
Script CPU time: 1553.15 s
Cove CPU time: 22435.24 s
Scan speed: 2150.5 bp/sec

Cove analysis of tRNAs ended: Sun Sep 11 08:42:34 JST 2011

Summary
——–
Overall scan speed: 185915.6 bp/sec

tRNAs decoding Standard 20 AA: 31889
Selenocysteine tRNAs (TCA): 3157
Possible suppressor tRNAs (CTA,TTA): 243
tRNAs with undetermined/unknown isotypes: 1143
Predicted pseudogenes: 189250
——-
Total tRNAs: 225682

tRNAs with introns: 4279
| Ala-AGC: 6 | Ala-GGC: 33 | Ala-CGC: 2 | Ala-TGC: 29 | Gly-ACC: 122 | Gly-GCC: 1826 | Gly-CCC: 74 | Gly-TCC: 103 | Pro-GGG: 4 | Pro-CGG: 1 | Pro-TGG: 1 | Thr-AGT: 6 | Thr-GGT: 8 | Thr-CGT: 2 | Thr-TGT: 5 | Val-AAC: 3 | Val-GAC: 48 | Val-CAC: 4 | Val-TAC: 8 | Ser-AGA: 8 | Ser-GGA: 13 | Ser-CGA: 2 | Ser-TGA: 3 | Ser-ACT: 22 | Ser-GCT: 148 | Arg-ACG: 8 | Arg-GCG: 34 | Arg-CCG: 1 | Arg-TCG: 7 | Arg-CCT: 12 | Arg-TCT: 15 | Leu-AAG: 6 | Leu-GAG: 5 | Leu-CAG: 8 | Leu-TAG: 2 | Leu-CAA: 9 | Leu-TAA: 4 | Phe-AAA: 7 | Phe-GAA: 46 | Asn-ATT: 7 | Asn-GTT: 16 | Lys-CTT: 3 | Lys-TTT: 8 | Asp-ATC: 14 | Asp-GTC: 174 | Glu-CTC: 12 | Glu-TTC: 47 | His-ATG: 6 | His-GTG: 49 | Gln-CTG: 7 | Gln-TTG: 11 | Ile-AAT: 3 | Ile-GAT: 8 | Ile-TAT: 6 | Met-CAT: 2 | Tyr-ATA: 15 | Tyr-GTA: 94 | Supres-CTA: 6 | Supres-TTA: 9 | Cys-ACA: 136 | Cys-GCA: 918 | Trp-CCA: 55 | SelCys-TCA: 28 |

Isotype / Anticodon Counts:

Ala : 335 AGC: 61 GGC: 54 CGC: 41 TGC: 179
Gly : 8382 ACC: 209 GCC: 2631 CCC: 3431 TCC: 2111
Pro : 37 AGG: 14 GGG: 3 CGG: 7 TGG: 13
Thr : 55 AGT: 25 GGT: 3 CGT: 11 TGT: 16
Val : 258 AAC: 39 GAC: 55 CAC: 99 TAC: 65
Ser : 637 AGA: 41 GGA: 214 CGA: 22 TGA: 29 ACT: 39 GCT: 292
Arg : 932 ACG: 51 GCG: 234 CCG: 57 TCG: 26 CCT: 308 TCT: 256
Leu : 174 AAG: 17 GAG: 25 CAG: 19 TAG: 14 CAA: 63 TAA: 36
Phe : 414 AAA: 48 GAA: 366
Asn : 81 ATT: 13 GTT: 68
Lys : 358 CTT: 110 TTT: 248
Asp : 360 ATC: 37 GTC: 323
Glu : 2493 CTC: 618 TTC: 1875
His : 237 ATG: 36 GTG: 201
Gln : 173 CTG: 104 TTG: 69
Ile : 68 AAT: 28 GAT: 18 TAT: 22
Met : 46 CAT: 46
Tyr : 967 ATA: 179 GTA: 788
Supres: 243 CTA: 162 TTA: 81
Cys : 14133 ACA: 1930 GCA: 12203
Trp : 1749 CCA: 1749
SelCys: 3157 TCA: 3157

6. Make fasta file from tRNAscan-SE output using the script below.

#!/usr/bin/perl

use strict;
use warnings;

my $usage = "Usage: $0 <trnascan.output> <genome.fa> <output.fa>\n";
my $infile = shift or die $usage;
my $genome = shift or die $usage;
my $outfile = shift or die $usage;

my %genome = ();
my $current_acc = '';
open(IN,'<',$genome) || die "Could not open $genome: $!\n";
TEST: while(<IN>){
   chomp;
   #if you are working with bacterial genomes change the line below to
   #if (/^>([A-Za-z0-9\|\.\_]+)\s+/){
   if (/^>(.*)$/){
      $current_acc = $1;
   } else {
      $genome{$current_acc} .= $_;
   }
   #last TEST unless $current_acc eq 'chr1';
}
close(IN);

warn "Stored genome into memory\n";

my %trna = ();
open (IN, '<', $infile) or die "Could not open $infile: $!\n";
open (OUT,'>', $outfile) or die "Could not open $outfile for writing: $!\n";
while (<IN>){
   chomp;
#Sequence                tRNA            Bounds          tRNA    Anti    Intron Bounds   Cove
#Name            tRNA #  Begin           End             Type    Codon   Begin   End     Score
#--------        ------  ----            ------          ----    -----   -----   ----    ------
#chr1            1       32999           33071           Pseudo  GCA     0       0       20.17
#chr1            2       34208           34280           Pseudo  GCA     0       0       31.02
#chr1            3       34729           34800           Pseudo  GAA     0       0       24.24
#chr1            3924    115232036       115232108       Lys     TTT     0       0       79.55
   next if (/^Sequence/);
   next if (/^Name/);
   next if (/^--/);
   next if (/^$/);
   my $strand = '+';
   my ($acc, $trna_number, $trna_start, $trna_end, $trna_type, $anticodon, $intron_start, $intron_end, $cove_score) = split(/\s+/);

   #tRNAscan likes to add a space to the coordinates randomly
   $trna_start =~ s/\s+//g;
   $trna_end =~ s/\s+//g;
   if ($trna_start > $trna_end){
      $strand = '-';
      ($trna_start, $trna_end) = ($trna_end, $trna_start);
   }

   my $length = $trna_end - $trna_start;
   #substr is 0 based and coordinates in tRNAscan-SE are 1 based
   my $seq = substr($genome{$acc},$trna_start-1,$length+1);
   $seq = uc($seq);
   if ($strand eq '-'){
      $seq =~ tr/ACGT/TGCA/;
      $seq = reverse($seq);
   }
   #my $id = "tRNA_${trna_number}_${acc}_${trna_start}_${trna_end}_${strand}_${trna_type}_${anticodon}_${cove_score}";
   if ($acc =~ /chrUn/){
      $acc =~ s/^chr//;
   }
   my $id = "${acc}_${trna_start}_${trna_end}_${strand}_${trna_type}_${anticodon}_${cove_score}";
   print OUT ">$id\n$seq\n";

}
close(IN);
close(OUT);

exit(0);

7. Run RepeatMasker on the fasta file using these parameters:

./RepeatMasker -s -species cow bos_tau_5_trnascan_se_short_id.fa

8. Some stats from RepeatMasker

cat bos_tau_5_trnascan_se_short_id.fa.out | grep -P "\s+\d+" | awk '{print $5}' | sort -u | wc
 225667  225667 9513704
cat bos_tau_5_trnascan_se_short_id.fa | grep ">" | wc
 225682  225682 9740003
cat bos_tau_5_trnascan_se_short_id.fa.out | grep -P "\s+\d+" | awk '{print $11}' | sort | uniq -c | sort -k1rn | less
 149028 SINE/BovA
  74586 SINE/tRNA-Glu
   2264 tRNA
    126 LTR/ERV1
     19 LINE/L1
     13 Simple_repeat
     12 Low_complexity
     10 SINE/MIR
      2 SINE/tRNA
      2 Unknown
      1 SINE/RTE-BovB
      1 snRNA

Perl script called remove_repeat.pl used to remove putative tRNA entries tagged as something other than tRNA by RepeatMasker

#!/usr/bin/perl

use strict;
use warnings;

my $fa = 'bos_tau_5_trnascan_se_short_id.fa';
my $rm = 'bos_tau_5_trnascan_se_short_id.fa.out';

my %repeat = ();
open(IN,'<',$rm) || die "Could not open $rm: $!\n";
while(<IN>){
   chomp;
   next unless /^\s+\d+/;
   my @line_split = split(/\s+/);
   next if $line_split[11] eq 'tRNA';
   $repeat{$line_split[5]} = '1';
}
close(IN);

open(IN,'<',$fa) || die "Could not open $fa: $!\n";
while(<IN>){
   chomp;
   if (/^>(.*)/){
      my $defline = $1;
      my $seq = <IN>;
      chomp($seq);
      next if (exists $repeat{$defline});
      print ">$defline\n$seq\n";
   } else {
      die "Error in line $.: $_\n";
   }
}
close(IN);
remove_repeat.pl | grep ">" | wc
   2239    2239   93827
remove_repeat.pl | grep ">" | grep -v "Pseudo" | wc
   1257    1257   50982
remove_repeat.pl | grep ">" | grep -v "Pseudo" | grep -v ">Un" | wc
   1121    1121   45380
remove_repeat.pl | grep ">" | grep -v "Pseudo" | grep -v ">Un" | cut -f5 -d'_' | sort | uniq -c | sort -k1rn | less
    253 Glu
    105 Lys
     75 Gly
     66 Val
     63 Ala
     59 Leu
     55 Asp
     46 Ser
     43 Gln
     42 Arg
     35 Asn
     34 Thr
     31 Met
     31 Phe
     30 Ile
     29 Cys
     28 Pro
     24 Tyr
     19 SeC
     17 Undet
     16 His
     13 Trp
      6 Sup
      1 SeC(e)
remove_repeat.pl | grep ">" | grep -v "Pseudo" | grep -v ">Un" | cut -f5,6 -d'_' | sort | uniq -c | sort -k1rn | less
    199 Glu_TTC
     64 Lys_TTT
     54 Glu_CTC
     50 Asp_GTC
     41 Lys_CTT
     37 Ala_AGC
     33 Asn_GTT
     31 Gly_TCC
     31 Met_CAT
     29 Cys_GCA
     29 Phe_GAA
     28 Val_CAC
     26 Gly_GCC
     26 Val_AAC
     23 Gln_CTG
     23 Tyr_GTA
     20 Gln_TTG
     19 Ile_AAT
     19 SeC_TCA
     19 Thr_AGT
     17 Gly_CCC
     17 Ser_GCT
     17 Undet_???
     16 Ala_TGC
     16 His_GTG
     16 Leu_TAA
     14 Leu_CAG
     14 Pro_AGG
     13 Arg_CCT
     13 Leu_AAG
     13 Ser_AGA
     13 Trp_CCA
     11 Leu_CAA
     11 Val_TAC
     10 Arg_TCG
     10 Ile_TAT
     10 Thr_TGT
      9 Pro_TGG
      8 Ala_CGC
      8 Arg_ACG
      7 Arg_TCT
      7 Ser_TGA
      5 Asp_ATC
      5 Leu_TAG
      5 Pro_CGG
      5 Ser_CGA
      5 Sup_TTA
      5 Thr_CGT
      4 Arg_CCG
      4 Ser_ACT
      2 Ala_GGC
      2 Asn_ATT
      2 Phe_AAA
      1 Gly_ACC
      1 Ile_GAT
      1 SeC(e)_TCA
      1 Sup_CTA
      1 Tyr_ATA
      1 Val_GAC

8. Comparison with results from the genomic tRNA database.

tRNA Tang et al., 2009 Genomic tRNA database
Glu_TTC 199 450
Lys_TTT 64 81
Glu_CTC 54 171
Asp_GTC 50 48
Lys_CTT 41 40
Ala_AGC 37 30
Asn_GTT 33 40
Gly_TCC 31 398
Met_CAT 31 33
Cys_GCA 29 353
Phe_GAA 29 28
Val_CAC 28 45
Gly_GCC 26 62
Val_AAC 26 23
Gln_CTG 23 42
Tyr_GTA 23 38
Gln_TTG 20 18
Ile_AAT 19 17
SeC_TCA 19 25
Thr_AGT 19 14
Gly_CCC 17 1315
Ser_GCT 17 18
Ala_TGC 16 34
His_GTG 16 22
Leu_TAA 16 9
Leu_CAG 14 9
Pro_AGG 14 12
Arg_CCT 13 111
Leu_AAG 13 11
Ser_AGA 13 14
Trp_CCA 13 182
Leu_CAA 11 15
Val_TAC 11 22
Arg_TCG 10 12
Ile_TAT 10 8
Thr_TGT 10 10
Pro_TGG 9 9
Ala_CGC 8 17
Arg_ACG 8 14
Arg_TCT 7 44
Ser_TGA 7 6
Asp_ATC 5 6
Leu_TAG 5 5
Pro_CGG 5 5
Ser_CGA 5 7
Sup_TTA 5 4
Thr_CGT 5 7
Arg_CCG 4 14
Ser_ACT 4 4
Ala_GGC 2 2
Asn_ATT 2 0
Phe_AAA 2 7
Gly_ACC 1 23
Ile_GAT 1 0
SeC(e)_TCA 1 25
Sup_CTA 1 12
Tyr_ATA 1 10
Val_GAC 1 0

Spearman’s rho of 0.748.

9. Lastly to find evolutionary conserved tRNAs, download the human, mouse, rat, fugu, chicken, horse, chimpanzee and dog genomes. Following the procedures listed above, run tRNAscan-SE on the respective genomes, produce tRNA fasta files for each respective genome using the Perl script in step 6. Remove intronic sequences detected by tRNAscan-SE using this script:

#!/usr/bin/perl

use strict;
use warnings;

my $usage = "Usage: $0 <tRNAscanOutput> <input.fa> <output.fa>\n";
my $infile = shift or die $usage;
my $fasta = shift or die $usage;
my $outfile = shift or die $usage;

my $acc = '';
my %seq = ();
open (IN,'<',$fasta) or die "Can't open $fasta: $!\n";
while (<IN>){
   chomp;
   #>tRNA_9_chr1_164648_164720_+_Pseudo_GCA_25.07
   if (/^>(.*)$/){
      $acc = $1;
   }
   elsif (/[AGTCagtc]+/){
      $seq{$acc} .= $_;
   }
   else {
      die "Can't parse line $.: $_\n";
   }
}
close(IN);

my $counter = '0';
my $intron = '0';
my $no_intron = '0';

open (OUT, '>', $outfile) or die "Can't open $outfile for writing: $!\n";
open (IN, '<', $infile) or die "Can't open $infile\n";
while (<IN>){
   chomp;
   next if (/^Sequence/);
   next if (/^Name/);
   next if (/^--/);
   next if (/^$/);
   next if (/Pseudo/);
   ++$counter;
   #Sequence                tRNA            Bounds          tRNA    Anti    Intron Bounds   Cove
   #Name            tRNA #  Begin           End             Type    Codon   Begin   End     Score
   #--------        ------  ----            ------          ----    -----   -----   ----    ------
   #chr1            1       32999           33071           Pseudo  GCA     0       0       20.17
   my ($acc, $trnaNo, $trnaStart, $trnaEnd, $trnaType, $antiCodon, $intronStart, $intronEnd, $coveScore) = split(/\s+/);
   if ($acc =~ /chrUn/){
      $acc =~ s/chr//;
   }
   #tRNAscan likes to add a space to the coordinates randomly
   $trnaStart =~ s/\s+//g;
   $trnaEnd =~ s/\s+//g;
   my $start = $trnaStart;
   my $end = $trnaEnd;
   my $strand = '+';
   if ($trnaStart > $trnaEnd){
      $strand = '-';
      ($trnaStart, $trnaEnd) = ($trnaEnd, $trnaStart);
   }
   if ($intronStart > $intronEnd){
      ($intronStart,$intronEnd) = ($intronEnd,$intronStart);
   }
   #my $id = "tRNA_${trnaNo}_${acc}_${trnaStart}_${trnaEnd}_${strand}_${trnaType}_${antiCodon}_${coveScore}";
   my $id = "${acc}_${trnaStart}_${trnaEnd}_${strand}_${trnaType}_${antiCodon}_${coveScore}";
   #Check if tRNAs are annotated correctly
   unless (exists $seq{$id}){
      die "tRNA $id does not exist in fasta\n";
   }
   my $sequence = $seq{$id};
   my $length = length($sequence);
   if ($intronStart != 0 && $intronEnd != 0){
      ++$intron;
      #Offset starts at 0 for substring
      #Coordinate for end of first part
      my $firstPartCo = $intronStart - $trnaStart;
      #Coordinate for start of second part
      my $secondPartCo = $intronEnd - $trnaStart;
      #Coordinate for end of second part
      my $thirdPartCo = $trnaEnd - $intronEnd;
      #IF intron is on negative strand need to start counting from the other end
      if ($strand eq '-'){
         $sequence = reverse($sequence);
         #Offset starts at 0 for substring
         #Coordinate for end of first part
         my $firstPartCo = $intronStart - $trnaStart;
         #Coordinate for start of second part
         my $secondPartCo = $intronEnd - $trnaStart;
         #Coordinate for end of second part
         my $thirdPartCo = $trnaEnd - $intronEnd;

         my $first = substr($sequence,0,$firstPartCo);
         $first = reverse($first);
         #0 based so add 1 to move 1 base to the right
         my $second = substr($sequence,$secondPartCo+1,$thirdPartCo);
         $second = reverse($second);
         my $tempSeq = $second . $first;
         $length = length ($tempSeq);
         print OUT ">$id\n";
         print OUT "$second$first\n";
      }
      else {
         my $first = substr($sequence,0,$firstPartCo);
         #0 based so add 1 to move 1 base to the right
         my $second = substr($sequence,$secondPartCo+1,$thirdPartCo);
         my $tempSeq = $first . $second;
         $length = length ($tempSeq);
         print OUT ">$id\n";
         print OUT "$first$second\n";
      }
   }
   else {
      ++$no_intron;
      print OUT ">$id\n";
      print OUT "$sequence\n";
   }
}
close (IN);
close (OUT);
warn "Processed $counter tRNAs and there were $intron tRNAs with intronic sequence and $no_intron without\n";

exit(0);

Using each respective fasta file (where there are no intronic sequences), make BLAST databases and finally BLAST the initial cow tRNA fasta set to each respective tRNA database from the different genomes. We will use a sequence identity of >= 95% as a threshold for orthology.

/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/human/hg19_trnascan_se.output genome/human/hg19.fa hg19_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/mouse/mm9_trnascan_se.output genome/mouse/mm9.fa mm9_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/rat/rn4_trnascan_se.output genome/rat/rn4.fa rn4_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/fugu/fr2_trnascan_se.output genome/fugu/fr2.fa fr2_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/horse/equ_cab_2_trnascan_se.output genome/horse/equ_cab_2.fa equ_cab_2_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/dog/canfam2_trnascan_se.output genome/dog/canfam2.fa canfam2_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/chicken/galgal3_trnascan_se.output genome/chicken/gal_gal_3.fa galgal3_trnascan_se.fa
/home/tan118/work/tRNAscan/trnascan_se_to_fasta.pl genome/chimpanzee/pan_tro_3_trnascan_se.output genome/chimpanzee/panTro3.fa pan_tro_3_trnascan_se.fa

makeblastdb -in hg19_trnascan_se_no_intron.fa -dbtype nucl -title human_trna -out human_trna
makeblastdb -in mm9_trnascan_se_no_intron.fa -dbtype nucl -title mouse_trna -out mouse_trna
makeblastdb -in rn4_trnascan_se_no_intron.fa -dbtype nucl -title rat_trna -out rat_trna
makeblastdb -in fr2_trnascan_se_no_intron.fa -dbtype nucl -title fugu_trna -out fugu_trna
makeblastdb -in equ_cab_2_trnascan_se_no_intron.fa -dbtype nucl -title horse_trna -out horse_trna
makeblastdb -in canfam2_trnascan_se_no_intron.fa -dbtype nucl -title dog_trna -out dog_trna
makeblastdb -in galgal3_trnascan_se_no_intron.fa -dbtype nucl -title chicken_trna -out chicken_trna
makeblastdb -in pan_tro_3_trnascan_se_no_intron.fa -dbtype nucl -title chimpanzee_trna -out chimpanzee_trna

blastn -db chicken_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_chicken_trna_no_intron.blastout
blastn -db mouse_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_mouse_trna_no_intron.blastout
blastn -db horse_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_horse_trna_no_intron.blastout
blastn -db dog_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_dog_trna_no_intron.blastout
blastn -db rat_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_rat_trna_no_intron.blastout
blastn -db human_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_human_trna_no_intron.blastout
blastn -db fugu_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_fugu_trna_no_intron.blastout
blastn -db chimpanzee_trna -outfmt 7 -evalue 1 -query bos_tau_5_trnascan_se_no_intron.fa -out bos_tau_5_trna_to_chimpanzee_trna_no_intron.blastout

To find conserved tRNAs and the corresponding phylogenetic profile, use this script:

#!/usr/bin/perl

use strict;
use warnings;

my $usage = "Usage: $0 <threshold i.e. 1.00 -> 0.50>\n";
my $trna_fasta = 'bos_tau_5_trnascan_se_no_intron.fa';
my $threshold = shift or die $usage;

my @files = qw/
bos_tau_5_trna_to_chicken_trna_no_intron.blastout
bos_tau_5_trna_to_dog_trna_no_intron.blastout
bos_tau_5_trna_to_fugu_trna_no_intron.blastout
bos_tau_5_trna_to_horse_trna_no_intron.blastout
bos_tau_5_trna_to_human_trna_no_intron.blastout
bos_tau_5_trna_to_mouse_trna_no_intron.blastout
bos_tau_5_trna_to_rat_trna_no_intron.blastout
bos_tau_5_trna_to_chimpanzee_trna_no_intron.blastout
/;

my %phyloProfile = ();
my %tRNAlist = ();
open (IN, '<',"$trna_fasta") or die "Can't open $trna_fasta: $!\n";
while (<IN>){
   chomp;
   if (/^>(.*)/){
      my $trna = $1;
      my $seq = <IN>;
      chomp($seq);
      my $length = length($seq);
      $tRNAlist{$trna} = $length;
   } else {
      die "Could not parse line $.: $_\n";
   }
}
close (IN);

foreach my $currentFile (@files){
   open (IN, '<', $currentFile) or die "Can't open $currentFile: $!\n";
   while (<IN>){
      chomp;
      next if /^#/;
      my ($query, $subject, $id, $aliLen, $mismatch, $gapOpen, $qStart, $qEnd, $sStart, $sEnd, $eValue, $bitScore) = split (/\t/);
      if (exists $tRNAlist{$query}){
         my $percentMatch = $id * $aliLen / '100';
         if ($percentMatch =~ /(\d+)\.(\d+)/){
            my $first = $1;
            my $second = $2;
            if ($second =~ /^9/){
               $percentMatch = $first + 1;
            }
            elsif ($second =~ /^0/){
               $percentMatch = $first;
            }
         }
         if (($percentMatch / $tRNAlist{$query}) >= $threshold){
            $phyloProfile{$query}{$currentFile} = '1';
         }
      }
      else {
         die "Query $query is missing\n";
      }
   }
   close (IN);
}

foreach my $trna (keys %tRNAlist){
   print "$trna\t";
   my $line = '';
   foreach my $files (@files){
      if (exists $phyloProfile{$trna}{$files}){
         $line .= '1';
      }
      else {
         $line .= '0';
      }
   }
   my $letter = substr($line,0,1);
   substr($line,0,1,"1$letter");
   print "$line\n";
}

The output from the above script should look something like this:

tRNA_2015_chrX_56490796_56490867_+_Undet_???_33.67 100000000
tRNA_3466_chr23_30023502_30023573_-_Glu_CTC_76.61 111111111
tRNA_1640_chr24_40048717_40048789_+_Gly_GCC_22.38 100000000

To identify the list of perfectly conserved tRNA genes, the above script with a threshold of 1 (i.e. 100%).

phylogenetic_profile_thresholding.pl 1 | grep 111111111 > perfect.list

Then use this script to extract the sequences:

#!/usr/bin/perl

use strict;
use warnings;

#>BTA1.trna335-ValAAC (5010368-5010440) Val (AAC) 73 bp Sc: 47.34

my $current = '';
my %fasta = ();
my %annotation =();
open(IN,'<','bos_tau_5_trnascan_se_no_intron.fa') || die "Could not open bos_tau_5_trnascan_se_no_intron.fa: $!\n";
while(<IN>){
   chomp;
   next if (/^$/);
   if (/^>(.*)$/){
      $current = $1;
      next;
   } else {
      $fasta{$current} .= $_;
   }
}
close(IN);

open(IN,'<','perfect.list') || die "Could not open perfect.list: $!\n";

while(<IN>){
   chomp;
   next if (/^$/);
   #BTA23.trna1281-IleAAT   111111111
   my ($id,$junk) = split(/\t/);
   if (exists $fasta{$id}){
      print ">$id\n$fasta{$id}\n";
   } else {
      die "Missing id $id\n";
   }
}
close(IN);

I came up with 171 tRNA genes (compared with 150 in the paper) using 100% conservation and the unique tRNA gene count was 26 (compared with 22 in the paper). Reasons for the increase are most likely due to improved genome assemblies.

And at the 0.95 threshold, there were 439 cow/bovine tRNA genes.

Print Friendly, PDF & Email



Creative Commons License
This work is licensed under a Creative Commons
Attribution 4.0 International License
.

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.