物种丰度计算

由于Qiime出了点问题，ITS项目先缓几天，这两天先忙着做meta的内容。

物种丰度计算准备工作：

1　　使用SOAPAligner对过滤好的数据进行比对，得到相应的.soap文件，里面记录匹配到的reads的情况；

2　　还需要将所有用到的reference做一个TAX，tax文件记录reference的物种信息，每一行分别是一个物种的gi号、界、门、纲、目、科、属、种、亚种的名称，分别用制表符隔开；

3　　对于每个亚种genome，需要计算其长度，做成SIZE文件，每一行分别是一个亚种的名称和genome长度，用制表符隔开；

计算过程：

以亚种为基础，一个亚种的丰度Ab(G) = Ab(U) + Ab(M)，其中

G：物种

U：unique数目，一条reads只比对上单一物种称为unique reads，一个物种的所有unique总和为U

M：multiple数目，一条reads比对上多个物种成为multiple reads，一个物种的所有multiple总和为M

Ab(U) = U/L(g)； L(g)为相应物种的genome长度；

Ni 表示 multiple reads 比对上的所有物种；

  1 #!/usr/bin/perl
  2 use strict;
  3 
  4 my $usage = "usage:get_profiling.pl <.soap file> <outprefix> <TAX file> <SIZE file> depth list\n";
  5    $usage .= "depth 1..8 stand for Kindom..SubSpecies\n";
  6 
  7 die $usage unless @ARGV>=5;
  8 
  9 my $soapfile = shift @ARGV;
 10 my $outprefix = shift @ARGV;
 11 my $taxfile = shift @ARGV;
 12 my $sizefile = shift @ARGV;
 13 my @depth = @ARGV;
 14 
 15 
 16 #########################################################################################################################
 17 #                                                            #
 18 #    %tax{gi}[Kindom][Phylum][Class][Order][Family][Genus][Species][SubSpecies]                    #
 19 #        1    2    3    4    5    6    7    8                        #
 20 #    %size{strname} = length of genome                                        #
 21 #    %soap{reads id}[0] = unique or multiple;                                    #
 22 #    %soap{reads id}[1..n] stand for matched strname;                                #
 23 #                                                            #
 24 #                    set tables                                    #    
 25 #                                                            #
 26 #########################################################################################################################
 27 open TAX,$taxfile or die $!;
 28 my %tax;
 29 while(<TAX>){
 30     chomp;
 31     my @a = split /\s+/;
 32     for(my $i=1;$i<=8;$i++){
 33         $tax{
   $a[0]}[$i]=$a[$i];