<jupyter_start><jupyter_text>MaAsLin2 Prepare data<jupyter_code>!mkdir /tmp/maaslin2
matrix = sdf.groupby(['genus', 'sample'])["abundance"].sum().rename('abundance').reset_index()
matrix = pd.pivot_table(sdf, index='sample', columns='genus', values='abundance', fill_value=0)

outfile = "../tmp/maaslin2/bracken_self_combine_result.genus_sum.matrix.tsv"
matrix.to_csv(outfile, sep="\t")
matrix<jupyter_output><empty_output><jupyter_text>**MaAsLin2 is run in `maaslin3_run.ipynb`** Supplementary figure 4: Treat MaAsLin2 results<jupyter_code>def read_maaslin_eth(fname):
    ref = os.path.basename(os.path.dirname(fname).split('_')[4]).title()
    df = pd.read_csv(fname, sep='\t')
    df['ref'] = df['metadata'].apply(lambda metadata: ref if metadata == 'Ethnicity' else 'None')
    return df

fnames = [
    '../tmp/maaslin2/abundance_genus_mean_cplm_chinese/significant_results.tsv',
    '../tmp/maaslin2/abundance_genus_mean_cplm_indian/significant_results.tsv',
    '../tmp/maaslin2/abundance_genus_mean_cplm_malay/significant_results.tsv'
]

sigres = pd.concat([
    read_maaslin_eth(fname) for fname in fnames
])

sigres['id'] = sigres[['feature', 'value', 'ref']].apply(lambda row: hash(tuple(sorted(row.values))), axis=1)
sigres = sigres.drop_duplicates('id')

sigres = sigres[['feature', 'metadata', 'value', 'ref', 'coef', 'stderr', 'pval', 'qval']].sort_values('pval')
sigres.columns = ['Genus', 'Feature', 'Value', 'Reference', 'Coefficient', 'StandardError', 'pValue', 'qValue']

sigres.head()
abu = sdf.groupby(['genus', 'sample', 'ethnicity'])["abundance"].sum().rename('abundance').reset_index()
mabu = abu.groupby('genus')['abundance'].sum().rename('abundance').sort_values(ascending=False)
top10 = list(mabu.head(10).index)
mabu.head(10)
saf = sigres[(sigres['Genus'].str.strip().isin(set(top10))) & (sigres['Feature'].str.strip() == 'Ethnicity')]
saf = saf[['Genus', 'Value', 'Reference', 'pValue']]
saf = saf.to_records(index=False)
saf
def get_panno_star(pvalue):
    annos = [[1e-4, "****"], [1e-3, "***"], [1e-2, "**"], [0.05, "*"]]
    for anno, name in annos:
        if pvalue <= anno:
            return name
    return 'ns'

def add_pvalue_hue(ax, pdata, odata, x, y, hue):
    # https://stackoverflow.com/questions/36578458/how-does-one-insert-statistical-annotations-stars-or-p-values-into-matplotlib
    
    sumx = odata.groupby(x)[y].max().to_dict()
    xorder = list(odata.sort_values(x)[x].unique())
    
    horder = list(odata.sort_values(hue)[hue].unique())
    # need to change that, don't have time right now
    horder = {horder[0]: -.3, horder[1]: 0, horder[2]: .3}
    
    height = (ax.get_ylim()[1] - ax.get_ylim()[0]) / 20
    found = defaultdict(int)
    
    for x, h1, h2, pvalue in pdata:
        xidx = xorder.index(x)
        h1idx = horder[h1]
        h2idx = horder[h2]       
        
        x1 = xidx + h1idx
        x2 = xidx + h2idx
        
        xheight = height * found[x] * 2
        found[x] += 1
        
        y1 = sumx[x] + xheight + height
        y2 = sumx[x] + xheight + height * 1.5
        y3 = sumx[x] + xheight + height * 1.5
        
        ax.plot([x1, x1, x2, x2], [y1, y2, y2, y1], lw=1.5, color='black')
        
        star = get_panno_star(pvalue)
        ax.text((x1 + x2) * .5, y3, star, ha='center', va='bottom')
abu_top10 = abu[abu['genus'].isin(set(top10))]

abu_top10['genus'] = pd.Categorical(abu_top10['genus'], categories=top10)
abu_top10['ethnicity'] = pd.Categorical(abu_top10['ethnicity'], categories=['Chinese', 'Indian', 'Malay'])

graph = Graph(abu_top10)
graph.sns.boxplot(x='genus', y='abundance', hue='ethnicity', palette=cpal[3:], fliersize=1)
graph.change_boxplot_width(.8)

add_pvalue_hue(graph.ax, saf, abu_top10,  x="genus", y="abundance", hue="ethnicity")

graph.ax.set_ylim(-2, 60)

graph.ax.set_xlabel('')
graph.ax.set_ylabel('Relative abundance', size=12)
graph.ax.legend(title='')

graph.set_size_inches(14, 5)
graph.tight_layout()<jupyter_output><ipython-input-12-4e5be3a7423b>:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  abu_top10['genus'] = pd.Categorical(abu_top10['genus'], categories=top10)
<ipython-input-12-4e5be3a7423b>:4: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  abu_top10['ethnicity'] = pd.Categorical(abu_top10['ethnicity'], categories=['Chinese', 'Indian', 'Malay'])<jupyter_text>Supplementary figure 8.A<jupyter_code>def set_label(row):
    if row['label'] == 'null':
        return 'SLC ' + str(row['SLCID'])
    else:
        return row['species']

ssdf = sdf.copy()
ssdf['label'] = ssdf['species'].fillna('null')
ssdf['label'] = ssdf.apply(set_label, axis=1)

species_median = ssdf.groupby('label')['abundance'].median()
top10 = species_median.sort_values(ascending=False).head(10)
top10
top_abu = ssdf[ssdf['label'].isin(set(list(top10.index)))]
top_abu['label'] = pd.Categorical(top_abu['label'], categories=list(top10.index))

cdef = "#9dbcd4"

graph = Graph(top_abu)
graph.sns.boxplot(x='label', y='abundance', color=cdef)
graph.apply_xticklabels(rotation=25, ha='right', style='italic')
graph.ax.set_xlabel('')
graph.ax.set_ylabel('Relative abundance')

graph.set_size_inches(6, 5)
graph.tight_layout()<jupyter_output><ipython-input-14-2e40b85e1b20>:2: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  top_abu['label'] = pd.Categorical(top_abu['label'], categories=list(top10.index))
{'fraction_non_alphanumeric': 0.17775210430764152, 'fraction_numerical': 0.028222478956923584, 'mean_word_length': 4.872093023255814, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 14, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '196553', 'license': 'MIT', 'n_tokens_mistral': 2481, 'n_tokens_neox': 2382, 'n_words': 396, 'path': 'scripts/bracken_self_combine.ipynb', 'repo_name': 'Chiamh/SPMP', 'size': '6059'}