由于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 }
时间： 2024-11-13 02:36:47