物种丰度计算
由于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]; 34 } 35 } 36 close TAX; 37 open SIZE,$sizefile or die $!; 38 my %size; 39 while(<SIZE>){ 40 chomp; 41 my @a = split /\s+/; 42 if (exists $tax{$a[0]}) 43 { 44 $size{$tax{$a[0]}[8]} = int($a[1]); 45 } 46 } 47 close SIZE; 48 open SOAP,$soapfile or die $!; 49 my %soap; 50 while(<SOAP>){ 51 chomp; 52 my @a = split /\s+/; 53 my $id = $a[0]; 54 $id = substr($id,0,$id.length()-2); 55 my $strname = $tax{$a[7]}[8]; 56 57 if ( exists $soap{$id} ){ 58 $soap{$id}[0] = "multiple" unless ($soap{$id}[1] eq $strname); 59 }else{ 60 $soap{$id}[0] = "unique"; 61 } 62 push(@{$soap{$id}},$strname); 63 64 my $reads2 = <SOAP>; 65 } 66 close SOAP; 67 ######################################################################################################################### 68 # # 69 # table setted successful # 70 # # 71 ######################################################################################################################### 72 my %uniquenum; 73 my %multiplestr; 74 my %str_reads; 75 my %reads_str; 76 my %abu; 77 my %abm; 78 my %sum; 79 my %ab_str; #the result hash; 80 foreach my $id (sort keys %soap){ 81 my $type = shift @{$soap{$id}}; 82 if ( $type eq "unique" ){ 83 $uniquenum{$soap{$id}[0]}++; 84 }elsif ($type eq "multiple"){ 85 foreach my $strname($soap{$id}){ 86 $multiplestr{$strname}++; 87 $reads_str{$id}{$strname}++; 88 $str_reads{$strname}{$id}++; 89 } 90 } 91 } 92 foreach my $strname(sort keys %uniquenum){ 93 $abu{$strname} = $uniquenum{$strname} / $size{$strname} if(exists $size{$strname}); 94 } 95 foreach my $id(sort keys %soap){ 96 foreach my $strname (sort keys %{$reads_str{$id}}){ 97 $sum{$id} += $abu{$strname} if (exists $abu{$strname}); 98 } 99 } 100 foreach my $strname(sort keys %multiplestr){ 101 foreach my $id(sort keys %{$str_reads{$strname}}){ 102 $abm{$strname} += $abu{$strname}/$sum{$id} if($sum{$id}); 103 } 104 } 105 foreach my $strname(sort keys %abu){ 106 if(exists $abm{$strname}) 107 { 108 $ab_str{$strname} = $abu{$strname} + $abm{$strname}; 109 }else{ 110 $ab_str{$strname} = $abu{$strname}; 111 } 112 } 113 undef %tax; 114 undef %size; 115 undef %soap; 116 undef %uniquenum; 117 undef %multiplestr; 118 undef %str_reads; 119 undef %reads_str; 120 undef %abu; 121 undef %abm; 122 undef %sum; 123 ######################################################################################################################### 124 # # 125 # abundance of subSpecies calculated successful # 126 # # 127 ######################################################################################################################### 128 129 open OUT,">test" or die $!; 130 foreach my $keys(sort keys %ab_str){ 131 print OUT "$keys\t$ab_str{$keys}\n"; 132 } 133 close OUT; 134 foreach my $d (@depth){ 135 die "Please input correct depth !\n" unless ($d>=1 && $d<=8); 136 &Abundance($d); 137 } 138 sub Abundance{ 139 my $depth = shift @_; 140 my @text=("num","Kindom","Phylum","Class","Order","Family","Genus","Species","SubSpecies"); 141 open TAX,$taxfile or die $!; 142 my %strtable; 143 my %ab; 144 while(<TAX>){ 145 chomp; 146 my @a = split /\s+/; 147 $strtable{$a[$depth]}{$a[8]}++; 148 } 149 foreach my $name (sort keys %strtable){ 150 foreach my $strname (sort keys %{$strtable{$name}}){ 151 if(exists $ab_str{$strname}){ 152 $ab{$name} += $ab_str{$strname}; 153 } 154 } 155 } 156 open AB,">$outprefix"."_"."$text[$depth].abundance" or die $!; 157 print AB "$text[$depth]\tabundance\n"; 158 foreach my $name(sort keys %ab){ 159 print AB "$name\t$ab{$name}\n" 160 } 161 close AB; 162 }