logo标签列表

从pandas DataFrame计算p值

pythonpandasdataframescipystatistics
6

我有一个名为stats的DataFrame，它具有Multindex和8个样本（这里只显示了两个）以及每个样本的8个基因。

 In[13]:stats
    Out[13]: 
                       ARG/16S                                            \
                         count          mean           std           min   
    sample      gene                                                       
    Arnhem      IC        11.0  2.319050e-03  7.396130e-04  1.503150e-03   
                Int1      11.0  7.243040e+00  6.848327e+00  1.364879e+00   
                Sul1      11.0  3.968956e-03  9.186019e-04  2.499074e-03   
                TetB       2.0  1.154748e-01  1.627663e-01  3.816936e-04   
                TetM       4.0  1.083125e-04  5.185259e-05  5.189226e-05   
                blaOXA     4.0  4.210963e-06  3.783235e-07  3.843571e-06   
                ermB       4.0  4.111081e-05  7.894879e-06  3.288865e-05   
                ermF       4.0  2.335210e-05  4.519758e-06  1.832037e-05   
    Basel       Aph3a      4.0  7.815592e-06  1.757242e-06  5.539389e-06   
                IC        11.0  5.095161e-03  5.639278e-03  1.302205e-03   
                Int1      12.0  1.333068e+01  1.872207e+01  4.988048e-02   
                Sul1      11.0  1.618617e-02  1.988817e-02  2.970397e-03

我正在尝试计算p值（学生t检验），比较每个基因之间的这些样本。
我已经使用了scipy.stats.ttest_ind_from_stats，但我只能得到一个基因的不同样本的p值，而且只有相邻的样本的p值。
Experiments = list(values1_16S['sample'].unique())
for exp in Experiments:
    if Experiments.index(exp)<len(Experiments)-1:
        second = Experiments[Experiments.index(exp)+1]
    else:
        second = Experiments[0]
    tstat, pvalue = scipy.stats.ttest_ind_from_stats(stats.loc[(exp,'Sul1')]['ARG/16S','mean'],
                                    stats.loc[(exp,'Sul1')]['ARG/16S','std'],
                                    stats.loc[(exp,'Sul1')]['ARG/16S','count'],
                                    stats.loc[(second,'Sul1')]['ARG/16S','mean'],
                                    stats.loc[(second,'Sul1')]['ARG/16S','std'],
                                    stats.loc[(second,'Sul1')]['ARG/16S','count'])
    d.append({'loc1':exp, 'loc2':second, 'pvalue':pvalue})


stats_Sul1 = pd.DataFrame(d)
stats_Sul1

如何获取所有样本之间的p值？是否有一种方法可以同时对所有基因进行此操作，而不是逐个运行每个基因的代码？

每个样本的基因都一样吗？我可以看到第二个样本有 Aph3a 基因，而第一个样本中没有。 - Ben.T
是的，每个样本都有相同的基因。第一个样本中Aph3为0，因此它在这里不出现。 - Gabriela Catalina
你可以查看我的编辑，我猜这个解决方案将适用于任何数量的样本和基因，即使它们不相同。 - Ben.T
1个回答
10

假设您对Y个样本有相同的X基因。我使用X = 3和Y = 2来尝试我的方法，但我想您可以推广。我从以下内容开始：

df1 = 
             count       mean        std       min
sample gene                                       
Arnhem IC       11   0.002319   0.000740  0.001503
       Int1     11   7.243040   6.848327  1.364879
       Sul1     11   0.003969   0.000919  0.002499
Basel  IC       11   0.005095   0.005639  0.001302
       Int1     12  13.330680  18.722070  0.049880
       Sul1     11   0.016186   0.019888  0.002970

请注意，基因需要按照相同的顺序排列。

首先使用 reset_index() 和 df_reindex = df1.reset_index()，我不确定在多索引情况下是否可行：

df_reindex =
   sample  gene  count       mean        std       min
0  Arnhem    IC     11   0.002319   0.000740  0.001503
1  Arnhem  Int1     11   7.243040   6.848327  1.364879
2  Arnhem  Sul1     11   0.003969   0.000919  0.002499
3   Basel    IC     11   0.005095   0.005639  0.001302
4   Basel  Int1     12  13.330680  18.722070  0.049880
5   Basel  Sul1     11   0.016186   0.019888  0.002970

我创建了一个滚动的DF并将其与df_reindex连接：

nb_genes = 3
df_rolled = pd.DataFrame(pd.np.roll(df_reindex,nb_genes,0), columns = df_reindex.columns)
df_joined = df_reindex.join(df_rolled, rsuffix='_')
# rsuffix='_' is to be able to perform the join

现在，在同一行上，我拥有你需要计算pvalue并使用apply创建列的所有数据:

df_joined['pvalue'] = df_joined.apply(lambda x: stats.ttest_ind_from_stats(x['mean'],x['std'],x['count'], x['mean_'],x['std_'],x['count_'])[1],axis=1)

最后，我使用你想要的数据创建了一个DF并重新命名列：

df_output = df_joined[['sample','sample_','gene','pvalue']].rename(columns = {'sample':'loc1', 'sample_':'loc2'})

您最终得到的数据为：

df_output = 
     loc1    loc2  gene    pvalue
0  Arnhem   Basel    IC  0.121142
1  Arnhem   Basel  Int1  0.321072
2  Arnhem   Basel  Sul1  0.055298
3   Basel  Arnhem    IC  0.121142
4   Basel  Arnhem  Int1  0.321072
5   Basel  Arnhem  Sul1  0.055298

您可以根据需要重新索引。

如果您希望对每个样本进行比较，我认为可以使用循环for来完成。

编辑: 使用pivot_table，我认为有更简单的方法。

以您的输入stats作为仅包含ARG / 16S（不确定如何处理此级别）的多层索引表格，因此我从以下内容开始（这可能是您的stats ['ARG / 16S']）：

df=
               count       mean           std       min
sample gene                                            
Arnhem IC         11   0.002319  7.396130e-04  0.001503
       Int1       11   7.243040  6.848327e+00  1.364879
       Sul1       11   0.003969  9.186019e-04  0.002499
       TetB        2   0.115475  1.627663e-01  0.000382
       TetM        4   0.000108  5.185259e-05  0.000052
       blaOXA      4   0.000004  3.783235e-07  0.000004
       ermB        4   0.000041  7.894879e-06  0.000033
       ermF        4   0.000023  4.519758e-06  0.000018
Basel  Aph3a       4   0.000008  1.757242e-06  0.000006
       IC         11   0.005095  5.639278e-03  0.001302
       Int1       12  13.330680  1.872207e+01  0.049880
       Sul1       11   0.016186  1.988817e-02  0.002970

使用函数 pivot_table，您可以重新排列数据，例如：

df_pivot = df.pivot_table(values = ['count','mean','std'], index = 'gene', 
                               columns = 'sample', fill_value = 0)
在这个 df_pivot 中（我为了易读性没有在此处打印，但在新列中会展示），您可以使用 itertools 和 apply 为每个 (sample1, sample2) 创造一列。
import itertools
for sample1, sample2 in itertools.combinations(df.index.levels[0],2):
    # itertools.combinations create all combinations between your samples
    df_pivot[sample1+ '_' + sample2 ] = df_pivot.apply(lambda x: stats.ttest_ind_from_stats(x['mean'][sample1],x['std'][sample1],x['count'][sample1], 
                                                                                        x['mean'][sample2 ],x['std'][sample2 ],x['count'][sample2 ],)[1],axis=1).fillna(1)

我认为这种方法与样本数量、基因数量无关，即使基因不完全相同，最终得到的 df_pivot 结果如下：

        count            mean                      std            Arnhem_Basel
sample Arnhem Basel    Arnhem      Basel        Arnhem      Basel             
gene                                                                          
Aph3a       0     4  0.000000   0.000008  0.000000e+00   0.000002     1.000000
IC         11    11  0.002319   0.005095  7.396130e-04   0.005639     0.121142
Int1       11    12  7.243040  13.330680  6.848327e+00  18.722070     0.321072
Sul1       11    11  0.003969   0.016186  9.186019e-04   0.019888     0.055298
TetB        2     0  0.115475   0.000000  1.627663e-01   0.000000     1.000000
TetM        4     0  0.000108   0.000000  5.185259e-05   0.000000     1.000000
blaOXA      4     0  0.000004   0.000000  3.783235e-07   0.000000     1.000000
ermB        4     0  0.000041   0.000000  7.894879e-06   0.000000     1.000000
ermF        4     0  0.000023   0.000000  4.519758e-06   0.000000     1.000000

如果它有效，请告诉我。

编辑2：为了回复这条评论，我认为你可以这样做：

df_pivot 不需要更改，然后创建一个多索引 DF df_multi 来写入你的结果：

df_multi = pd.DataFrame(index = df.index.levels[1], 
                        columns = pd.MultiIndex.from_tuples([p for p in itertools.combinations(df.index.levels[0],2)])).fillna(0)

然后您使用循环for来实现这个df_multi中的数据：

for sample1, sample2 in itertools.combinations(df.index.levels[0],2):
    # itertools.combinations create all combinations between your samples
    df_multi.loc[:,(sample1,sample2)] = df_pivot.apply(lambda x: stats.ttest_ind_from_stats(x['mean'][sample1],x['std'][sample1],x['count'][sample1], 
                                                                                        x['mean'][sample2 ],x['std'][sample2 ],x['count'][sample2 ],)[1],axis=1).fillna(1)

最后，您可以在第1级上使用 transpose 和 unstack 来获得您所要求的方式（如果我误解了，请关闭）

df_output = df_multi.transpose().unstack(level=[1]).fillna(1)

如果你想看到索引中的最后一个样本和列中的第一个样本（因为它们在我构建的过程中并不存在），你需要将 itertools.combinations替换为itertools.combinations_with_replacement，分别在创建df_multi和循环for中进行替换（我没有尝试过，但应该可以正常工作）。

谢谢你的编辑！第一个解决方案对我没有起作用。现在尝试您的编辑，但我必须先删除第一行MultiIndex ('ARG/16S')。如果我使用drop()删除它，它也会删除下面的所有内容。如果我让它起作用，我会告诉你的。 - Gabriela Catalina
如果你执行df = stats['ARG/16S']，可以吗？ - Ben.T
成功了！非常感谢！您知道是否还有一种方法可以将比较的样本作为列和行索引，并使用基因作为列Multiindex吗？（只是为了更容易找到p值） - Gabriela Catalina
@GabrielaCatalina 是的，有方法，我在考虑创建输出数据的最佳方式。我会在我的EDIT2中尝试一下。 - Ben.T
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接