<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The purpose of this notebook is to investigate the Schubert results. Specifically, it looks like 100% of tests are being rejected by qvalue at an alpha cutoff of 0.1. That seems fishy, let's investigate!
Step1:  It looks like OTUs with an uncorrected pvalue of 0.9 get smushed down to 0.08 with qvalue - this seems fishy!
Let's plot a few of these bugs and see if we can understand what's happening with the covariate...
Step2:  So these are, for the most part, singletons, maybe?
Let's check that I calculated the ubiquity correctly here..
Step3:  Hm, ok - so it looks like they're the same. I didn't mess up the ubiquity calculation...
It looks like the least significant bugs (which get moved down to below q = 0.1) are in the lower covariate group - they have ubiquity = 0.02 - 0.03. I want to look back at the covariate boxplots and see if this makes sense, I guess...?
Step4:  Uncorrected pvalue vs. corrected pvalue?
Step5:  Y-axis is the corrected qvalue (specified in subplot title), x-axis is the original pvalue
How does qvalue compare to the other more methods?
Specifically, the methods that have more reasonable results...
Step6:  How does effect size correlate with pvalue and adjusted pvals?
Step7:  And covariate?
<ASSISTANT_TASK:>
Python Code:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set_style('white')
fpvals = 'schubert-sb-table.txt'
fotu = 'data/cdi_schubert_results/RDP/cdi_schubert.otu_table.100.denovo.rdp_assigned'
fmeta = 'data/cdi_schubert_results/cdi_schubert.metadata.txt'
pvals = pd.read_csv(fpvals, sep=' ')
pvals.index = pvals['otu']
df = pd.read_csv(fotu, sep='\t', index_col=0).T
abundf = df.divide(df.sum(axis=1), axis=0)
meta = pd.read_csv(fmeta, sep='\t', index_col=0)
pvals.head()
plt.hist(pvals['unadjusted'])
df.head()
meta.head()
pvals.columns
pvals.sort_values(by='qvalue', ascending=False)[['pval', 'qvalue', 'ubiquity']].head(15)
Explanation: The purpose of this notebook is to investigate the Schubert results. Specifically, it looks like 100% of tests are being rejected by qvalue at an alpha cutoff of 0.1. That seems fishy, let's investigate!
End of explanation
# Tidyfy the OTU table
df.index.name = 'sample'
tidydf = df.reset_index().melt(id_vars='sample', var_name='otu', value_name='abun')
# Add disease state
tidydf = tidydf.join(meta['DiseaseState'], on='sample')
tidydf.head()
otus = pvals.sort_values(by='qvalue', ascending=False).index[0:12].tolist()
fig, ax = plt.subplots(3, 4, figsize=(14,12))
ax = ax.flatten()
for i in range(len(ax)):
    o = otus[i]
    sns.stripplot(
        data=tidydf.query('otu == @o'), 
        x='DiseaseState', y='abun', 
        ax=ax[i],
        jitter=True)
    
# Tidyfy the realtive abundance OTU table
abundf.index.name = 'sample'
tidyabundf = abundf.reset_index().melt(id_vars='sample', var_name='otu', value_name='abun')
# Add disease state
tidyabundf = tidyabundf.join(meta['DiseaseState'], on='sample')
fig, ax = plt.subplots(3, 4, figsize=(14,12))
ax = ax.flatten()
for i in range(len(ax)):
    o = otus[i]
    sns.stripplot(
        data=tidyabundf.query('otu == @o'), 
        x='DiseaseState', y='abun', 
        ax=ax[i],
        jitter=True)
Explanation: It looks like OTUs with an uncorrected pvalue of 0.9 get smushed down to 0.08 with qvalue - this seems fishy!
Let's plot a few of these bugs and see if we can understand what's happening with the covariate...
End of explanation
kept_otus = pvals['otu'].tolist()
df_fromR = df.loc[:, kept_otus]
keep_dis = ['H', 'CDI', 'nonCDI']
df_fromR = df_fromR.loc[meta.query('DiseaseState == @keep_dis').index, :]
df_fromR.shape
# Hm, okay - maybe my problems are coming from
df_fromR.shape, df_fromR.dropna().shape
# Recalculate ubiquity with python
fig, ax = plt.subplots(2, 2, figsize=(10, 8))
ax = ax.flatten()
ax[0].hist((df_fromR.dropna() > 0).sum() / df_fromR.dropna().shape[0])
ax[0].set_title('Ubiquity calculated from python')
ax[1].hist(pvals['ubiquity'])
ax[1].set_title('Ubiquity calculated from R')
ax[2].hist(np.log10((df_fromR.dropna() > 0).sum() / df_fromR.dropna().shape[0]))
ax[2].set_title('Ubiquity calculated from python')
ax[3].hist(np.log10(pvals['ubiquity']))
ax[3].set_title('Ubiquity calculated from R')
Explanation: So these are, for the most part, singletons, maybe?
Let's check that I calculated the ubiquity correctly here..
End of explanation
np.log10(0.02)
Explanation: Hm, ok - so it looks like they're the same. I didn't mess up the ubiquity calculation...
It looks like the least significant bugs (which get moved down to below q = 0.1) are in the lower covariate group - they have ubiquity = 0.02 - 0.03. I want to look back at the covariate boxplots and see if this makes sense, I guess...?
End of explanation
pvals.columns
pcols = [u'bonf', u'bh',
       u'qvalue', u'ihw-a10', 
       u'bl-df03', u'lfdr']
fig, ax = plt.subplots(2, 3, figsize=(14, 8))
ax = ax.flatten()
i = 0
for c in pcols:
    ax[i].scatter(pvals['unadjusted'], pvals[c])
    ax[i].set_title(c)
    if i > 2:
        ax[i].set_xlabel('unadjusted')
    if i in [0, 3]:
        ax[i].set_ylabel('corrected')
        
    i += 1
Explanation: Uncorrected pvalue vs. corrected pvalue?
End of explanation
fig, ax = plt.subplots(1, 4, figsize=(14, 4))
pcols = ['unadjusted', 'ihw-a10', 'bl-df03', 'lfdr']
i = 0
for c in pcols:
    ax[i].scatter(pvals['qvalue'], pvals[c])
    ax[i].set_title(c)
    ax[i].axhline(0.1)
    ax[i].set_xlabel('qvalue')
    i += 1
pvals.columns
Explanation: Y-axis is the corrected qvalue (specified in subplot title), x-axis is the original pvalue
How does qvalue compare to the other more methods?
Specifically, the methods that have more reasonable results...
End of explanation
pcols = ['qvalue', 'unadjusted', 'ihw-a10', 'bl-df03', 'lfdr']
fig, ax = plt.subplots(1, len(pcols), figsize=(15, 4))
i = 0
for c in pcols:
    ax[i].scatter(pvals['effect_size'], pvals[c])
    ax[i].set_title(c)
    ax[i].set_xlabel('effect size')
    i += 1
Explanation: How does effect size correlate with pvalue and adjusted pvals?
End of explanation
fig, ax = plt.subplots(1, len(pcols), figsize=(15, 4))
i = 0
for c in pcols:
    ax[i].scatter(pvals['ubiquity'], -np.log10(pvals[c]))
    ax[i].set_title(c)
    ax[i].set_xlabel('ubiquity')
    #ax[i]
    i += 1
print('Ubiquity vs log10(pvalue)')
import numpy as np
fig, ax = plt.subplots(1, len(pcols), figsize=(15, 4))
i = 0
for c in pcols:
    ax[i].scatter(pvals['ubiquity'].rank(), pvals[c])
    ax[i].set_title(c)
    ax[i].set_xlabel('ubiquity')
    #ax[i]
    i += 1
print('Rank ubiquity vs pvalue')
Explanation: And covariate?
End of explanation
<END_TASK>