【正文】 With score: Query: hsa:51144 hsa:2052 hsa:2053 hsa:51004 hsa:9420 %_id: || || || |0| || Sbjct: gi4098957 gi207689 gi55930 gi1226240 Hit 4: Ortholog pathway for: Caenorhabditis elegans. With score: Query: hsa:51144 hsa:2052 hsa:2053 hsa:51004 hsa:9420 %_id: || || || || || Sbjct: gi726418 gi1465805 gi3876864 gi2088820 gi13775482 !/usr/bin/perl program: pm_ortholog purpose: finds an orthlogous pathway for a query pathway in a given species. Prints the output in alignment format. author: Grace Yang Pathmetrics, Inc. 10/14/2022 usage: pm_ortholog query_pw query_aa protein_db were contains the pathway information contains all the proteins in query use strict。 Part 1. Parse input, check files my ($usage, $q_id, $q_aa, $q_pnu, $q_pw, $aa_db)。 my (%gn2spec, %score, %total_score, $file)。 my (@q, @arr, %qu2spec, $spec, @time_st)。 $usage = \n $0 query_pw query_aa protein_db\n query_pw: query pathway file query_aa: query aa file protein_db: protein db to search\n\n。 if (@ARGV1) { die $usage。} ($q_pw, $q_aa, $aa_db)=@ARGV。 for $file ($q_pw, $q_aa, $aa_db) { if (!(e $file)) { die Did not find $file file\n。} } open (QSEQ, $q_pw)。 while (QSEQ) { $file=$_。 chomp ($file)。 if ($file=~/^(\S+)\s/) { $q_id=$1。 next。} push(@q,split(/\s+/, $file))。 $q_pnu=@q。 } close (QSEQ)。 @time_st=localtime。 amp。print_header。 amp。big_matrix_sort($aa_db, $q_aa)。 open (AA, /usr/local/biobin/im_retrieve $aa_db /tmp/$$. |)。 while (AA) { if ($_=~/^(\S+)\s+.*\[([\w\s]+)\]/) { $gn2spec{$1}=$2。}} close (AA)。 get the best hit for each query id and each spec open (MAT, /tmp/$$.)。 while(MAT) { chomp。 @arr = split(/\t/)。 if($qu2spec{$arr[0]}{$gn2spec{$arr[1]}}) {next。} $qu2spec{$arr[0]}{$gn2spec{$arr[1]}} = $arr[1]。 $score{$arr[0]}{$arr[1]} = $arr[2]。 if($total_score{$gn2spec{$arr[1]}}){ $total_score{$gn2spec{$arr[1]}} += $arr[2]*20。 }else{ $total_score{$gn2spec{$arr[1]}} = $arr[2]*20。} } close(MAT)。 my ($qid, $i, $j, $ln)。 $ii=0。 foreach $spec (sort by_score keys (%total_score)) { $ii++。 printf Hit%3d: Ortholog pathway for: %20s. With score: %\n\n, $ii,$spec, $total_score{$spec}。 for ($i=0。 $i(@q/6)。 $i++) { my (@ln1, @ln2, @ln3, $sc, $hid, $k)。 for ($j=0。 $j6。 $j++) { $k = $i*6+$j。 if ($k @q){ $sc = $score{$q[$k]}{$qu2spec{$q[$k]}{$spec}}。 if ($qu2spec{$q[$k]}{$spec}) {$hid=$qu2spec{$q[$k]}{$spec}。 } else {$hid =。} if (!defined($sc)) {$sc=。} push (@ln1,$q[$k])。push (@ln2, \|$sc\|)。 push (@ln3, $hid)。} } format STDOUT= Query: @|||||||||| @|||||||||| @|||||||||| @|||||||||| @|||||||||| @|||||||||| $ln1[0], $ln1[1], $ln1[2],$ln1[3],$ln1[4],$ln1[5] %_id: @||||| @||||| @||||| @||||| @||||| @||||| $ln2[0], $ln2[1], $ln2[2],$ln2[3],$ln2[4],$ln2[5] Sbjct: @|||||||||| @|||||||||| @|||||||||| @|||||||||| @|||||||||| @|||||||||| $ln3[0], $ln3[1], $ln3[2],$ln3[3],$ln3[4],$ln3[5] . write STDOUT。 } } amp。print_end。 sub by_score { return $total_score{$b}=$total_score{$a}。} sub big_matrix_sort { my (@arr, $q_len, $m_len, $pct_id, $pct_pos, $l, $tp)。 my ($bg, $end,$hsp_len,$pm_score)。 my ($aa_db, $qu_aa)=@_。 open (IN, /usr/local/biobin/im_cycle blastp $aa_db $q_aa S=100 | /usr/local/biobin/pm_pblast |)。 open(HIT, /tmp/$$.matrix)。 while(IN){ chomp。 @arr = split(/\t/)。 ($q_len, $m_len) = split(/:/,$arr[2])。 ($pct_id, $pct_pos) = split(/:/, $arr[5])。 ($l, $tp) = split(/:/, $arr[6])。 ($bg, $end) = split(//, $l)。 $hsp_len = abs($end$bg)+1。 $pm_score = get_pm_score($pct_id, $pct_pos, $hsp_len, $q_len, $m_len)。 if($pm_score = 0) { next。 } printf HIT %s\t%s\t%\n, $arr[0],$arr[1],$pm_score。 } close(IN)。close(HIT)。 system (sort k 3rn /tmp/$$.matrix /tmp/$$.)。 system (cut f2 /tmp/$$.matrix |sort u /tmp/$$.)。 } sub get_pm_score { my ($pct_id, $pct_pos, $hsp_len, $q_len, $m_len) = @_。 my $len = ($q_len$m_len) ? $q_len : $m_len。 if($len = 0) { print STDERR warn : length of sequence is calculated to = 0\n。 return 1。 }else{ return * ($pct_id + $pct_pos) * $hsp_len / $len。} } sub print_header { my ($aa_nu)。 print \n。 print PM_ORTHOLOG , Pathmetrics,[Oct202022] [Build linuxx86]\n\n。 print Ref.: US . \Methods for Establishing Pathway Database\n。 print and Perform Pathway Searches\. XXX Feb. 20, 2022.\n\n。 print Query pathway= $q_id\n。 print ($q_pnu proteins)\n\n。 print Database: $aa_db\n。 open (DB, $)。 while (DB) {if ($_=~/Total keys\s(\d+)/) {$aa_nu=$1。 last。}} close (DB)。 print $aa_nu proteins.\n。 } Pathway Prediction Engines ? They are the crown jewels of Pathmetrics software tools ? Can predict many novel interactions ? Use diverse input data, including sequence data, expression data, and known interaction data ? Employ plex numerical algorithms such as dynamical programming and clustering