GINI Index的计算方法

基尼系数是1943年美国经济学家阿尔伯特·赫希曼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。在RNA结构中使用GINI系数首先由Weissman在文献[1]中提出，他在使用这个系数时引用了文献[2]。而文献[2]既没有其他引用其他的文献，也没有给出GINI的计算方法，只是做了如下描述：

张老师的脚本中给出的计算方法如下所示：

 
1
sub calcGini 
2
{  
3
    my $ref_array = shift;
4
    my %parameters = @_;
5
    
6
    my $len = scalar ( @{$ref_array} );
7
    my $start = defined ( $parameters{start} ) ? $parameters{start} : 0;
8
    my $end = defined ( $parameters{end} ) ? $parameters{end} : $len - 1;
9
 
10
    my $gini = "NULL";
11
    my @num = (); 
12
    for ( my $idx = $start; $idx < $end; $idx++ ) { 
13
        if ( $ref_array->[$idx] ne "NULL" ) {
14
            push @num, $ref_array->[$idx];
15
        }
16
    }
17
    @num = sort { $a <=> $b } @num;
18
​
19
    my $count = scalar ( @num );
20
    if ( $count > 10 ) {
21
        my $total = eval ( join ( "+", @num ) );
22
        if ( $total > 0 )  {
23
            my $giniB = 0;  my $accum = 0;
24
            for ( my $idx = 0; $idx < scalar(@num); $idx++ ) { 
25
                $accum += $num[$idx];
26
                #$giniB += ( ( $accum - $num[$idx]/2 ) / $total / $count ) ; 
27
                $giniB += ( $accum / $total / $count ) ; 
28
            }
29
            if ( $parameters{format} ) {
30
                $gini = sprintf ( "%.3f", 1 - 2*$giniB );
31
            }   
32
            else {
33
                $gini = 1 - 2*$giniB;
34
            }
35
        }
36
    }       
37
 
38
    return $gini;
39
};

​ 上面有个小bug：my $end = defined ( $parameters{end} ) ? $parameters{end} : $len - 1;这一行使最后一个值不会被加入到数组中，在计算5'UTR、3'UTR、CDS时$parameters{end}这个参数通常是指定的，所以这个bug不会有任何影响，在计算整条RNA时，由于最后一个位点往往是NULL，这也没有影响。

​ 我们可以把上面的代码整理成公式：

​ 把上面的Perl代码修改为Python:

 
1
def GINI_Zhang(list_of_values):
2
    "We suppose all input values are valid float value"
3
    length = len(list_of_values)
4
    total = sum(list_of_values)
5
    if total == 0 or length <= 10: return -1
6
    Sorted_Array = sorted(list_of_values)
7
    accum = 0
8
    giniB = 0
9
    for i in Sorted_Array:
10
        accum += i
11
        giniB += float(accum)/total/length
12
    return 1 - 2*giniB

​ 我们在网上搜索到另外一个计算GINI系数的公式（目前使用的方法）：

 
1
def GINI_Web(list_of_values):
2
    "We suppose all input values are valid float value"
3
    length = len(list_of_values)
4
    total = sum(list_of_values)
5
    if total == 0 or length <= 10: return -1
6
    Sorted_Array = sorted(list_of_values)
7
    accum, giniB = 0, 0
8
    for i in Sorted_Array:
9
        accum += i
10
        giniB += accum - i / 2.0
11
    fair_area = accum * length / 2.0
12
    return (fair_area - giniB) / fair_area

我们把上面的两段代码整理成公式：

这两种计算方式的区别是：

从上面的公式可以看出张老师的计算方法得到的值要小，差1/n。下面我们给出几组测试值来验证两种方法得到的结果的区别：

>>> test_Array = [0.4,]*15 		 #15个0.4
>>> GINI_Zhang(test_Array)
-0.06666666666666665			#小于0，接近于0
>>> GINI_Web(test_Array)
1.578983857244667e-16			#更加接近于0

>>> test_Array = [1,1,] + [0,]*15
>>> GINI_Zhang(test_Array)
0.8235294117647058
>>> GINI_Web(test_Array)
0.8823529411764706				#更加接近于1

注意：序列越长，GINI Index越大

随机数生成短序列长度长序列长度图片
uniform(0.1, 0.9)201000
uniform(0.1, 0.9)1001000
uniform(0.1, 0.9)2001000
uniform(0.1, 0.9)5001000

结论

​ 当n的值比较大时，两种计算方法不会相差很大。所以，就算用张老师的计算方法和我们目前用的方法会得到非常相似的结果。

[1] Genome-wide probing of RNA structure reveals active unfolding of mRNA structures in vivo. Nature. 2014

[2] Wittebolle, L. et al. Initial community evenness favours functionality under selective stress. Nature. 2009