6.6 perl--根据rdp置信度修改命名

将Otu中的代表序列上传到网页版的RDP进行注释后，会返回一个文件，该文件为每条序列所能匹配上的物种的分类地位，同时在各个分类地位上都有置信度。如果是代表片段小于250bp，我们一般选择的是置信度大于50为可信，这样我们需要将那些置信度小于50的命名为Unclassfied。因此我写了这样一个脚本rdpformat.pl。

脚本难点与思路：

这个脚本难点在于你不知道每一行在什么分类水平上的置信度会小于50，所以你需要识别小于50的那个数。所以需要设置多重If循环，如果在门水平上置信度小于50，后面的都命名为Unclassfied,如果是在纲水平上，那么后面的分类水平为Unclassfied,

脚本编写：rdpformat.pl

#!/usr/bin/perl
use warnings;
use strict;
my $usage="
#########################################################Promgram#########################################

rdpformat.pl
This script is used to translate the the output of rdp website to unclassfied if the confidence <50%

History:

   2014/08/10    Sam    First release

example:

   perl /path/way/to/rdpformat.pl < input > < output >
###########################################################################################################
";
die $usage if (@ARGV!=2);
my ($input,$output)=@ARGV;

open IN,"$input" or die;
open OUT,'>',"$output" or die;

my @p=();
while (< IN >){
  chomp;
  my @p=split /;/;
  if ($p[5]<50){
    print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\tunclassfied_$p[2]\t$p[5]\tunclassfied_$p[2]\tp[7]\tunclassfied_$p[2]\t$p[9]\tunclassfied_$p[2]\t$p[11]\tunclassfied_$p[2]\t$p[13]\tunclassfied_$p[2]\t$p[15]\n";

  }else {

        if ($p[7]<50) {
        print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\tunclassfied_$p[4]\t$p[7]\tunclassfied_$p[4]\t$p[9]\tunclassfied_$p[4]\t$p[11]\tunclassfied_$p[4]\t$p[13]\tunclassfied_$p[4]\t$p[15]\n";

         }else {
             if ($p[9]<50) {
             print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\tunclassfied_$p[6]\t$p[9]\tunclassfied_$p[6]\t$p[11]\tunclassfied_$p[6]\t$p[13]\tunclassfied_$p[6]\t$p[15]\n";

             }else {
                  if ($p[11]<50) {
                  print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\t$p[8]\t$p[9]\tunclassfied_$p[8]\t$p[11]\tunclassfied_$p[8]\t$p[13]\tunclassfied_$p[8]\t$p[13]\n";

                   }else {

                       if (defined($p[15])){ 

                              if ($p[13]<50) {

                         print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\t$p[8]\t$p[9]\t$p[10]\t$p[11]\tunclassfied_$p[10]\t$p[13]\tunclassfied_$p[10]\t$p[15]\n";

                               }else {

                                       if ($p[15]<50) {

                               print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\t$p[8]\t$p[9]\t$p[10]\t$p[11]\t$p[12]\t$p[13]\tunclassfied_$p[12]\t$p[15]\n";

                                      }else {

                                print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\t$p[8]\t$p[9]\t$p[10]\t$p[11]\t$p[12]\t$p[13]\t$p[14]\t$p[15]\n";
                                      }

                                                   }
                                      }else{
                                  print OUT "$p[0]\t$p[1]\t$p[2]\t$p[3]\t$p[4]\t$p[5]\t$p[6]\t$p[7]\t$p[8]\t$p[9]\tunclassfied_$p[8]\t$p[11]\tunclassfied_$p[8]\t$p[11]\tunclassfied_$p[8]\t$p[11]\n";
                                          }
                                            }

                               }
                                                                  }
                   }
            }
close IN;
close OUT;

脚本存在的问题：

use of uninitialized values in concatenation(.) or string at rdpformat,pl at line line 36, < IN > line 3.

所以只好把 if ($p[11]<50) {

print …… $p[15]换成$p[13],问题的原因是不能识别输入的数据的第15行。好奇怪的东西额？？？

对rdp的数据预处理：

sed ‘s/%//’ all_**_classfied.txt |sed ‘s/”//g’ >abc:

#将数据中的%和””去掉。

ps:这应该是我的第4个脚本了吧，呵呵，越写月有感觉了。