# stats from blast hits table (btall1b = best 1 target source per ref prot) head arath15ap-pine9set.aa.btall1b Query Source Bits Ident Align Qlen Slen arath15:AT1G01010.1 pitaeg8:Tv1R002904t1 109 79 158 429 892 arath15:AT1G01010.1 pita2trpacbio:c768_f2p1_3288 109 79 158 429 866 arath15:AT1G01010.1 tsaGECO:01009440_1 109 79 158 429 869 arath15:AT1G01010.1 tsaGDQR:01010264_1 107 78 158 429 887 arath15:AT1G01010.1 pita1gmodl:PITA_000051105_RA 107 78 156 429 361 arath15:AT1G01010.1 pilav1gm:PILAhq_040488_RA 107 66 152 429 280 arath15:AT1G01010.1 tsaGBLJ:01035593_1 105 69 137 429 723 arath15:AT1G01010.1 pitav1tr:1A_I9_VO_L_1_T_114208_116761 105 69 137 429 741 arath15:AT1G01010.1 tsaGEUZ:01002827_1 105 70 137 429 707 --- head vitvi15nc-pine9set.aa.btall1b Query Source Bits Ident Align Qlen Slen vitvi15nc:NP_001267815.1 tsaGBLJ:01023584_1 808 399 535 537 552 vitvi15nc:NP_001267815.1 tsaGDQR:01020649_1 801 395 530 537 551 vitvi15nc:NP_001267815.1 tsaGECO:01026399_1 799 396 535 537 552 vitvi15nc:NP_001267815.1 pita2trpacbio:c7523_f1p2_2153 798 395 535 537 552 vitvi15nc:NP_001267815.1 pitav1tr:3A_I18_VO_L_2_T_124691_136250 798 396 535 537 572 vitvi15nc:NP_001267815.1 pita1gmodl:PITAhm_002825_RA 797 395 535 537 572 vitvi15nc:NP_001267815.1 pitaeg8:Tv1R010364t1 796 395 535 537 579 vitvi15nc:NP_001267815.1 tsaGEUZ:01005322_1 795 393 535 537 552 vitvi15nc:NP_001267815.1 pilav1gm:PILAhm_001377_RA 652 340 535 537 559 --- sort -k6,6nr -k1,1 -k5,5nr -k3,3nr arath*-pine9set.aa.btall1b | env mineq=0.3 nref=1 mainref=1 perl -ne 'BEGIN{ $MINEQ=$ENV{mineq}||0.5; }($rd,$td,$bs,$idn,$al,$rw,$tw)=split; next unless($rw>0); next if($ENV{mainref} and not($rd=~/\.1$/)); ($ts)= $td=~m/(\w+):/; if($rd ne $lrd) { $topbs=$bs; $topal=$al; $ok=1; } else { $ok=($al>=$MINEQ*$topal); } if($ok){ $sum{$ts}{n}++; $sum{$ts}{al}+=$al; $sum{$ts}{idn}+=$idn; $pi=int(100*$idn/$al); $sum{$ts}{pi}+=$pi; $pa=int(100*$al/$rw); $pa=100 if($pa>100); $sum{$ts}{pa}+=$pa; } $rid{$rd}++; $lrd=$rd; END{ @ts=sort keys %sum; $nrd=scalar(keys %rid); @st=qw(al idn pa pi); @hd=qw(Geneset nHit pHit Algn Ident pAln pIden); map{ $_=sprintf"%-12s",$_; } @hd[0,1]; print join("\t",@hd)." for nref=$nrd\n"; for $ts (@ts) { $n=$sum{$ts}{n}; $nd= ($ENV{nref})?$nrd:$n; @av=map{ $v=$sum{$ts}{$_}; $nnd=($_ eq "pi")?$n:$nd; $av=int(10*$v/$nnd)/10; $av; } @st; $pnd=int(1000*$n/$nrd)/10; print join("\t",sprintf("%-12s",$ts),"$n/$nd",$pnd,@av)."\n"; } print "\n"; } ' # ------------ vitvi15nc, mainref=0, nref=1, mineq=0.3 Geneset nHit pHit Algn Ident pAln pIden for nref=35618 pilav1gm 30981/35618 86.9 356.4 204.3 70 57.9 pita1gmodl 31837/35618 89.3 361.6 205.1 71 57.4 pita2trpacbio 30121/35618 84.5 314.2 182.1 66.7 57.6 pitaeg8 35094/35618 98.5 428.2 250 82.7 58.6 pitav1tr 31260/35618 87.7 380 214.5 72.7 56.1 tsaGBLJ 33750/35618 94.7 393.5 227.8 77.1 57.8 tsaGDQR 30200/35618 84.7 397.5 228.8 70.8 56.9 tsaGECO 34530/35618 96.9 426.8 249.5 82.4 58.3 tsaGEUZ 32152/35618 90.2 396.1 232.3 76.7 58.6 #------------------ arath15ap, mainref=1 nref=1 mieq=0.3 Geneset nHit pHit Algn Ident pAln pIden for nref=23161 pilav1gm 20275/23161 87.5 317.4 183.6 70.1 57.6 pita1gmodl 20741/23161 89.5 318 181.6 70.5 57 pita2trpacbio 19792/23161 85.4 289.3 167.1 67.1 57 pitaeg8 22652/23161 97.8 374.1 219.3 80.9 58.1 pitav1tr 19987/23161 86.2 332.8 188.8 70.8 55.9 tsaGBLJ 21322/23161 92 344.6 200 74.7 57.4 tsaGDQR 18830/23161 81.3 344.5 198.7 67.6 56.3 tsaGECO 22122/23161 95.5 375.6 220.3 80.4 57.8 tsaGEUZ 20553/23161 88.7 351.7 207.1 75.3 58.3 arath15ap, mainref=1 nref=0 mieq=0.3 Geneset nHit pHit Algn Ident pAln pIden for nref=23161 pilav1gm 20275/20275 87.5 362.6 209.7 80.1 57.6 pita1gmodl 20741/20741 89.5 355.2 202.8 78.8 57 pita2trpacbio 19792/19792 85.4 338.5 195.6 78.5 57 pitaeg8 22652/22652 97.8 382.5 224.2 82.7 58.1 pitav1tr 19987/19987 86.2 385.7 218.8 82.1 55.9 tsaGBLJ 21322/21322 92 374.3 217.2 81.2 57.4 tsaGDQR 18830/18830 81.3 423.7 244.4 83.1 56.3 tsaGECO 22122/22122 95.5 393.2 230.6 84.2 57.8 tsaGEUZ 20553/20553 88.7 396.3 233.4 84.9 58.3 arath15ap, mainref=1 nref=1 mieq=0.3, longest 1000 REF Geneset nHit pHit Algn Ident pAln pIden for nref=999 pilav1gm 803/999 80.3 809.9 473.6 60 56.8 pita1gmodl 830/999 83 860.5 493.5 63.2 56 pita2trpacbio 742/999 74.2 591.1 346.3 46.4 57 pitaeg8 994/999 99.4 1182.9 696.7 83.5 57.7 pitav1tr 866/999 86.6 971.5 556.8 70.3 55.6 tsaGBLJ 971/999 97.1 997 585.7 72.5 57.3 tsaGDQR 963/999 96.3 1132.6 663 80.4 57 tsaGECO 980/999 98 1131.2 671.8 81.8 57.7 tsaGEUZ 909/999 90.9 1009.6 602.8 74 57.8 #------------------ ## ordered by pHit arath15ap, mainref=1 nref=1 mieq=0.3 Geneset nHit pHit Algn Ident pAln pIden for nref=23161 pitaeg8 22652/23161 97.8 374.1 219.3 80.9 58.1 tsaGECO 22122/23161 95.5 375.6 220.3 80.4 57.8 tsaGBLJ 21322/23161 92 344.6 200 74.7 57.4 pita1gmodl 20741/23161 89.5 318 181.6 70.5 57 tsaGEUZ 20553/23161 88.7 351.7 207.1 75.3 58.3 pilav1gm 20275/23161 87.5 317.4 183.6 70.1 57.6 pitav1tr 19987/23161 86.2 332.8 188.8 70.8 55.9 pita2trpacbio 19792/23161 85.4 289.3 167.1 67.1 57 tsaGDQR 18830/23161 81.3 344.5 198.7 67.6 56.3 vitvi, mainref=0 nref=1 mieq=0.3 Geneset nHit pHit Algn Ident pAln pIden for nref=35618 pitaeg8 35094/35618 98.5 428.2 250 82.7 58.6 tsaGECO 34530/35618 96.9 426.8 249.5 82.4 58.3 tsaGBLJ 33750/35618 94.7 393.5 227.8 77.1 57.8 tsaGEUZ 32152/35618 90.2 396.1 232.3 76.7 58.6 pita1gmodl 31837/35618 89.3 361.6 205.1 71 57.4 pitav1tr 31260/35618 87.7 380 214.5 72.7 56.1 pilav1gm 30981/35618 86.9 356.4 204.3 70 57.9 tsaGDQR 30200/35618 84.7 397.5 228.8 70.8 56.9 pita2trpacbio 30121/35618 84.5 314.2 182.1 66.7 57.6 #------------------