Load environment

import pandas as pd
import glob
import seaborn as sns
import matplotlib.pyplot as plt
import yaml
import numpy as np
from collections import defaultdict
from scipy.stats import variation, entropy

with open("../data/resources/rcParams.yaml") as f:
    rcParamsDict = yaml.full_load(f)
    for k in rcParamsDict["rcParams"]:
        print("{} {}".format(k,rcParamsDict["rcParams"][k]))
        plt.rcParams[k] = rcParamsDict["rcParams"][k]
    for k1 in set(list(rcParamsDict)).difference(set(["rcParams"])):
        print("{} {}".format(k1,rcParamsDict[k1]))

figure.dpi 80
savefig.dpi 500
figure.figsize [10, 10]
axes.facecolor None
figure.facecolor None
dotSize 20

first_tranche = glob.glob("../../CensusSeq_Feb/output/census/**.txt")
second_tranche = glob.glob("../../CensusSeq_Mar/output/census/**.txt")

all_files = first_tranche + second_tranche

"chimeric": color : "#EF6F6C" oldName : "chimeric" oldNameHeader : "chimeric" newName : "chimeric"

color_palette = {
 'CTL01A': '#DBB807',
    'CTL08A': '#0FB248',
    'CTL04E': '#FF0054',
    'CTL02A': '#7B00FF',
'H9': '#72190E',
 'H1': '#994F88',
 'CTL05A': '#1965B0',
 'CTL07C': '#437DBF',
 'CTL06F': '#CAE0AB',
 'CTL09A': '#FFFF00',
 'KTD8_2': '#E65518',
 'UCSFi001-A': '#7BAFDE'}

Load and format data

all_results = {}
for file in all_files:
    name = file.strip('.txt').split('/')[-1].split('.')[0]
    all_results[name] = pd.read_csv(file, skiprows = 2, sep = '\t')

all_results_df = pd.concat(all_results.values(), keys = all_results.keys()).reset_index()
donor_map_names = {i:j for i, j in zip(all_results_df['DONOR'], all_results_df['DONOR'])}
donor_map_names['CHD2WT'] = 'UCSFi001-A'
donor_map_names['CHD8WT'] = 'H9'
all_results_df['DONOR'] = all_results_df['DONOR'].map(donor_map_names)

all_results_df.head()

	level_0	level_1	DONOR	REPRESENTATION
0	M18	0	CTL01A	0.94993
1	M18	1	CTL02A	0.00285
2	M18	2	CTL05A	0.00750
3	M18	3	CTL07C	0.01086
4	M18	4	H1	0.00647

all_results_df.shape

(582, 4)

metadata = pd.read_excel('../../data/csv/CensusSeq_metadata_new.xlsx')
metadata

	Sample name	Mix	timepoint
0	M1	1	day -2
1	M2	2	day -2
2	M3	3	day -2
3	M4	4	day -2
4	M5	5	day -2
...	...	...	...
90	M91	6	day 50
91	M92	6	day 50
92	M93	8	day 50
93	M94	8	day 50
94	M95	8	day 50

95 rows × 3 columns

metadata.index = metadata['Sample name']

metadata['Mix'].unique()

array([1, 2, 3, 4, 5, 6, 7, 8])

len(metadata['Sample name'].unique())

95

len(all_results_df['level_0'].unique())

95

all_results_df = all_results_df.drop('level_1', axis  = 1)
all_results_df.columns = ['Sample name', 'DONOR', 'REPRESENTATION']
all_results_df

	Sample name	DONOR	REPRESENTATION
0	M18	CTL01A	0.94993
1	M18	CTL02A	0.00285
2	M18	CTL05A	0.00750
3	M18	CTL07C	0.01086
4	M18	H1	0.00647
...	...	...	...
577	M60	H9	0.30743
578	M60	CTL04E	0.01336
579	M60	CTL06F	0.24393
580	M60	CTL08A	0.02211
581	M60	CTL09A	0.40072

582 rows × 3 columns

all_results_df.index = all_results_df['Sample name']
all_results_df['MIX ID'] = all_results_df['Sample name'].map({i: j for i, j in zip(metadata['Sample name'], metadata['Mix'])})
all_results_df['Timepoint'] = all_results_df['Sample name'].map({i: j for i, j in zip(metadata['Sample name'], metadata['timepoint'])})
all_results_df = all_results_df[all_results_df['MIX ID'] != 7]
#all_results_df['MIX ID'] = all_results_df['MIX ID'].map({1:1, 2:2, 3:3,4:4,5:5,6:6,8:7})
all_results_df.to_csv('../../data/csv/CensusSeq_data.csv')

all_results_df.head()

	Sample name	DONOR	REPRESENTATION	MIX ID	Timepoint
Sample name
M18	M18	CTL01A	0.94993	4	day 5
M18	M18	CTL02A	0.00285	4	day 5
M18	M18	CTL05A	0.00750	4	day 5
M18	M18	CTL07C	0.01086	4	day 5
M18	M18	H1	0.00647	4	day 5

all_results_df.shape

(546, 5)

Plot contributions from each samples

#color_palette

order = ['day -2', 'day 5', 'day 12', 'day 25', 'day 50']
tp = all_results_df['MIX ID'].unique().tolist()
tp.sort()

fig, ax = plt.subplots(2, 4, figsize = (20, 20), gridspec_kw={'wspace': 0.4, 'hspace': 0.4})
ax = ax.flatten().T

for mix, ax in zip(tp, ax):
    
    sub = all_results_df[all_results_df['MIX ID'] == mix]
    #print(sub)
    sub_df_pivoted = pd.pivot(sub, index = 'Sample name', columns='DONOR', values='REPRESENTATION')
    
    sub_df_pivoted.index = sub_df_pivoted.index.map({i: j for i, j in zip(sub['Sample name'], sub['Timepoint'])})
    

    sub_df_pivoted.loc[[i for i in order if i in sub_df_pivoted.index]].plot(kind = 'bar', stacked = True, color = color_palette, ax = ax)
    ax.legend(bbox_to_anchor = (1,1))
    ax.set_title(f'Mix ID: {mix}')
    ax.set_xlabel('Time point')
    ax.set_ylabel('Proportion of identities')
#plt.savefig('censusSeq_results.png')
plt.tight_layout()
plt.show()

/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
/usr/local/lib/python3.8/dist-packages/pandas/plotting/_matplotlib/tools.py:400: MatplotlibDeprecationWarning: 
The is_first_col function was deprecated in Matplotlib 3.4 and will be removed two minor releases later. Use ax.get_subplotspec().is_first_col() instead.
  if ax.is_first_col():
<ipython-input-18-9284e24e05b2>:23: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Lineplots of contributions trends

all_results_df['Timepoint_int'] = all_results_df.Timepoint.apply(lambda x: int(x.strip('day ')))
all_results_df['Timepoint_int'].sort_values
mix_ordered = all_results_df['MIX ID'].unique().tolist()
mix_ordered.sort()

fig, ax = plt.subplots(4, 2, figsize = (20, 25))
ax = ax.flatten().T

ax[-1].set_axis_off()

for mix, ax in zip(mix_ordered, ax):
    
    sub = all_results_df[all_results_df['MIX ID'] == mix].reset_index(drop = True)

    sns.lineplot(data = sub, x = 'Timepoint_int', y = 'REPRESENTATION', hue = 'DONOR', marker = 'o', palette=color_palette, ax = ax)
            
    ax.set_title(f'Mix ID: {mix}', fontsize = 25)
    ax.set_ylabel('Fraction represented', fontsize = 20)
    ax.set_xlabel('Timepoint', fontsize = 20)
    #start, end = (-3, 52)
    ax.xaxis.set_ticks([-2, 5, 12, 25, 50])

    ax.tick_params(axis='both', labelsize=15)
    
    
plt.tight_layout()
plt.savefig('./figures/Mix_CensusSeq.svg', bbox_inches = 'tight')

<ipython-input-19-a668c17f281d>:1: 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
  all_results_df['Timepoint_int'] = all_results_df.Timepoint.apply(lambda x: int(x.strip('day ')))

all_results_df = all_results_df.drop('Sample name', axis = 1).reset_index()
max_rep = all_results_df[all_results_df.Timepoint == 'day 5'].groupby('Sample name').max('REPRESENTATION')['REPRESENTATION']
all_results_df[all_results_df.Timepoint == 'day 5'][all_results_df[all_results_df.Timepoint == 'day 5'].REPRESENTATION.isin(max_rep)]

	Sample name	DONOR	REPRESENTATION	MIX ID	Timepoint	Timepoint_int
0	M18	CTL01A	0.94993	4	day 5	5
28	M19	CTL01A	0.91013	4	day 5	5
56	M23	CTL01A	0.75915	5	day 5	5
95	M30	CTL01A	0.89443	8	day 5	5
112	M31	CTL01A	0.88202	8	day 5	5
143	M22	CTL01A	0.79156	5	day 5	5
159	M10	UCSFi001-A	0.54225	1	day 5	5
191	M21	CTL01A	0.78139	5	day 5	5
206	M20	CTL01A	0.91905	4	day 5	5
217	M32	CTL01A	0.87633	8	day 5	5
253	M11	UCSFi001-A	0.52613	1	day 5	5
262	M9	UCSFi001-A	0.53311	1	day 5	5
273	M12	CTL01A	0.62220	2	day 5	5
311	M25	CTL01A	0.42265	6	day 5	5
331	M26	CTL01A	0.44736	6	day 5	5
367	M17	H9	0.52285	3	day 5	5
372	M14	CTL01A	0.62952	2	day 5	5
377	M13	CTL01A	0.62534	2	day 5	5
430	M16	H9	0.59398	3	day 5	5
442	M24	CTL01A	0.42961	6	day 5	5
529	M15	H9	0.60331	3	day 5	5

all_results_df.Timepoint.unique()

array(['day 5', 'day 12', 'day 25', 'day 50', 'day -2'], dtype=object)

fig, ax = plt.subplots(figsize = (20,5))
sns.boxplot(data = all_results_df, x = 'DONOR', y = 'REPRESENTATION', hue = 'Timepoint', hue_order=['day -2', 'day 5', 'day 12', 'day 25', 'day 50'])
ax.legend(bbox_to_anchor = (1,1))

<matplotlib.legend.Legend at 0x7fe92c149460>

Compute entropy of representations

Shannon's entropy is computed for each timepoint and each sample. It quantifies the expected uncertainty inherent in the possible outcomes of a discrete random variable, therefore the higher its value and and the higher the balance in representation. In fact, we observe a higher entropy for -2 time point.

entropy_df = all_results_df.groupby(['Timepoint_int', 'Sample name'])['REPRESENTATION'].apply(entropy).reset_index()
entropy_df.head()

	Timepoint_int	Sample name	REPRESENTATION
0	-2	M1	1.598650
1	-2	M2	1.584976
2	-2	M3	1.777457
3	-2	M4	1.776391
4	-2	M5	1.775167

fig, ax = plt.subplots(figsize = (5,8))
sns.barplot(data = entropy_df, x = 'Timepoint_int', y = 'REPRESENTATION', ax = ax, color = '#2a9d8f')
#sns.swarmplot(data = std, x = 'MIX ID', y = 'REPRESENTATION', hue = 'Timepoint_int', ax = ax, dodge = True)
ax.set_xticklabels(ax.get_xticklabels(), rotation = 90)

ax.set_ylabel('Shannon entropy', fontsize = 20)
ax.set_xlabel('Day', fontsize = 20)
ax.tick_params(axis='both', which='major', labelsize=20)

#plt.legend(bbox_to_anchor = (1,1), title = 'Day')
plt.tight_layout()
plt.savefig('./figures/Entropy.svg', bbox_inches = 'tight')

Compute coefficient of variation of representations for each donor

Coefficient of variation is computed for each timepoint and each donor (so variation for each donor at a certain time point across different mixes). Time point -2 is excluded because no replicate i available for that.
Additionally, for each donor the absolute difference in CV between one time point and the previous one is computed.

std = all_results_df.groupby(['Timepoint_int', 'DONOR'])['REPRESENTATION'].apply(variation).reset_index()
std = std[std.Timepoint_int != -2]

std = std.sort_values(by=["DONOR", "Timepoint_int"])
std["Difference"] = std.groupby("DONOR")["REPRESENTATION"].diff()
#std

fig, ax = plt.subplots(figsize = (20,10))
sns.barplot(data = std, x = 'DONOR', y = 'REPRESENTATION', hue = 'Timepoint_int', ax = ax, palette = sns.cubehelix_palette())
ax.set_xticklabels(ax.get_xticklabels(), rotation = 90)
plt.legend(bbox_to_anchor = (1,1))
plt.savefig('./figures/CV_timepoint_donor.svg', bbox_inches = 'tight')

fig, ax = plt.subplots(figsize = (20,10))
std['Difference'] = np.abs(std['Difference'])
sns.barplot(data = std, x = 'DONOR', y = 'Difference', hue = 'Timepoint_int', ax = ax, palette = sns.cubehelix_palette())
ax.set_xticklabels(ax.get_xticklabels(), rotation = 90)
plt.legend(bbox_to_anchor = (1,1))

<matplotlib.legend.Legend at 0x7fe92cfa3d00>

Compute coefficient of variation the of representations for each donor in each mix

Coefficient of variation (CV) is computed for each timepoint, each mix and each donor (so variation across different replicates). Time point -2 is excluded because no replicate i available for that.
Additionally, for each donor the absolute difference in CV between one time point and the previous one is computed.

std = all_results_df.groupby(['Timepoint_int', 'MIX ID', 'DONOR'])['REPRESENTATION'].apply(variation).reset_index()
std = std[std.Timepoint_int != -2]

std = std.sort_values(by=["MIX ID", 'DONOR', "Timepoint_int"])
std["Difference"] = std.groupby("MIX ID")["REPRESENTATION"].diff()
std

	Timepoint_int	MIX ID	DONOR	REPRESENTATION	Difference
42	5	1	CTL04E	0.050608	NaN
84	12	1	CTL04E	0.144417	0.093809
122	25	1	CTL04E	0.178294	0.033876
164	50	1	CTL04E	0.060064	-0.118230
43	5	1	CTL05A	0.036617	-0.023446
...	...	...	...	...	...
162	25	8	CTL04E	0.134457	0.091223
204	50	8	CTL04E	0.240799	0.106342
83	5	8	CTL08A	0.100209	-0.140590
163	25	8	CTL08A	0.263623	0.163414
205	50	8	CTL08A	0.416809	0.153186

164 rows × 5 columns

We plot here the distribution of the CV in each mix for each timepoint (each point of the distribution would be a donor)

fig, ax = plt.subplots(figsize = (20,10))
sns.boxplot(data = std, x = 'MIX ID', y = 'REPRESENTATION', hue = 'Timepoint_int', ax = ax, palette = sns.cubehelix_palette())
#ax.set_xticklabels(ax.get_xticklabels(), rotation = 90)
ax.set_ylabel('Coefficient of variation', fontsize = 20)
ax.set_xlabel('MIX', fontsize = 20)

ax.tick_params(axis='both', which='major', labelsize=20)
plt.tight_layout()
plt.legend(bbox_to_anchor = (1,1), title = 'Day')
plt.savefig('./figures/CV_timepoint_mix_donor.svg', bbox_inches = 'tight')

std['Difference'] = np.abs(std['Difference'])

fig, ax = plt.subplots(figsize = (20,10))
sns.boxplot(data = std, x = 'MIX ID', y = 'Difference', hue = 'Timepoint_int', ax = ax, palette = sns.cubehelix_palette())
ax.set_xticklabels(ax.get_xticklabels(), rotation = 90)
ax.set_ylabel('Absolute CV difference')
plt.legend(bbox_to_anchor = (1,1), title = 'Time range')

<matplotlib.legend.Legend at 0x7fe92c171250>

Compute weighted rank and weighted normalized rank

Day 5

all_timepoints_norm_weighted = {}

ranked_lists = {}
for m in all_results_df[all_results_df['Timepoint'] == 'day 5']['MIX ID'].unique():
    df = all_results_df[((all_results_df['Timepoint'] == 'day 5') & (all_results_df['MIX ID'] == m)  )].sort_values(by = 'REPRESENTATION', ascending = False)
    #print(df.groupby('DONOR').sum().sort_values(by = 'REPRESENTATION', ascending = False))
    result = df.groupby('DONOR').sum().sort_values(by = 'REPRESENTATION', ascending = False).reset_index()
    ranked_lists[f'Mix {m}'] = result.DONOR.tolist()

ranked_lists

{'Mix 4': ['CTL01A', 'KTD8_2', 'CTL07C', 'CTL05A', 'H1', 'CTL02A'],
 'Mix 5': ['CTL01A', 'UCSFi001-A', 'CTL06F', 'CTL04E', 'CTL05A', 'H1'],
 'Mix 8': ['CTL01A', 'CTL08A', 'CTL04E', 'CTL02A'],
 'Mix 1': ['UCSFi001-A', 'CTL06F', 'CTL04E', 'CTL05A', 'H1'],
 'Mix 2': ['CTL01A', 'CTL09A', 'CTL08A', 'CTL07C', 'CTL02A'],
 'Mix 6': ['CTL01A',
  'CTL09A',
  'UCSFi001-A',
  'CTL08A',
  'CTL06F',
  'CTL04E',
  'CTL07C',
  'CTL05A',
  'H1',
  'CTL02A'],
 'Mix 3': ['H9', 'CTL09A', 'UCSFi001-A', 'CTL08A', 'CTL06F', 'CTL04E']}

data = list(ranked_lists.values())

Weighted mean - not including percentages

d = defaultdict(list)
d

for l in data:
    #print(len(l))
    for idx, value in enumerate(l):
        d[value].append( (idx + 1) * len(l)) # or (idx + 1) * len(l)? which one is better?

mean_d = {}
for l in d:
    mean_d[l] = np.mean(d[l])

mean_d_df = pd.DataFrame(mean_d.values(), mean_d.keys())
mean_d_df.columns = ['mean_weighted_rank']
#mean_d_df = mean_d_df.drop(1)
mean_d_df.sort_values(by = 'mean_weighted_rank', ascending = True)

	mean_weighted_rank
H9	6.00
CTL01A	6.20
KTD8_2	12.00
CTL09A	14.00
UCSFi001-A	16.25
CTL08A	21.75
CTL06F	27.00
CTL04E	29.40
CTL07C	36.00
CTL05A	38.50
CTL02A	44.25
H1	45.25

mean_d_df.sort_values(by = 'mean_weighted_rank').to_csv('../../data/csv/CensusSeq_weighted_rank_d5.csv')

Weighted mean - including percentages

all_results_df[(all_results_df['Timepoint'] == 'day 5')].groupby('DONOR').sum('REPRESENTATION').sort_values(by = 'REPRESENTATION')

	REPRESENTATION	MIX ID	Timepoint_int
DONOR
CTL02A	0.06150	60	60
H1	0.10079	48	60
KTD8_2	0.11159	12	15
CTL07C	0.21065	36	45
CTL05A	0.27002	48	60
CTL04E	0.92864	69	75
CTL08A	1.00095	57	60
CTL06F	1.12763	45	60
H9	1.72014	9	15
CTL09A	1.93780	33	45
UCSFi001-A	2.58963	45	60
CTL01A	10.94067	75	75

mean_d_df['mean_Representation'] = all_results_df[(all_results_df['Timepoint'] == 'day 5')].groupby('DONOR').mean('REPRESENTATION')['REPRESENTATION']
mean_d_df['combined_scores'] = mean_d_df['mean_weighted_rank'] * mean_d_df['mean_Representation']
mean_d_df.sort_values('combined_scores')

	mean_weighted_rank	mean_Representation	combined_scores
CTL02A	44.25	0.005125	0.226781
H1	45.25	0.008399	0.380062
KTD8_2	12.00	0.037197	0.446360
CTL07C	36.00	0.023406	0.842600
CTL05A	38.50	0.022502	0.866314
CTL08A	21.75	0.083413	1.814222
CTL04E	29.40	0.061909	1.820134
CTL06F	27.00	0.093969	2.537167
CTL09A	14.00	0.215311	3.014356
H9	6.00	0.573380	3.440280
UCSFi001-A	16.25	0.215802	3.506791
CTL01A	6.20	0.729378	4.522144

mean_d_df.sort_values('combined_scores').to_csv('../../data/csv/CensusSeq_combined_weighted_rank_d5.csv')

np.round(mean_d_df.sort_values('combined_scores'), 2)

	mean_weighted_rank	mean_Representation	combined_scores
CTL02A	44.25	0.01	0.23
H1	45.25	0.01	0.38
KTD8_2	12.00	0.04	0.45
CTL07C	36.00	0.02	0.84
CTL05A	38.50	0.02	0.87
CTL08A	21.75	0.08	1.81
CTL04E	29.40	0.06	1.82
CTL06F	27.00	0.09	2.54
CTL09A	14.00	0.22	3.01
H9	6.00	0.57	3.44
UCSFi001-A	16.25	0.22	3.51
CTL01A	6.20	0.73	4.52

all_timepoints_norm_weighted['day 5'] = mean_d_df['combined_scores']

Day 12

ranked_lists = {}
for m in all_results_df[all_results_df['Timepoint'] == 'day 12']['MIX ID'].unique():
    df = all_results_df[((all_results_df['Timepoint'] == 'day 12') & (all_results_df['MIX ID'] == m)  )].sort_values(by = 'REPRESENTATION', ascending = False)
    result = df.groupby('DONOR').sum().sort_values(by = 'REPRESENTATION', ascending = False).reset_index()
    ranked_lists[f'Mix {m}'] = result.DONOR.tolist()

data = list(ranked_lists.values())

Weighted mean - not including percentages

d = defaultdict(list)
d

for l in data:
    for idx, value in enumerate(l):
        d[value].append( (idx + 1) * len(l)) # or (idx + 1) * len(l)? which one is better?

mean_d = {}
for l in d:
    mean_d[l] = np.mean(d[l])

mean_d_df = pd.DataFrame(mean_d.values(), mean_d.keys())
mean_d_df.columns = ['mean_weighted_rank']
mean_d_df.sort_values(by = 'mean_weighted_rank', ascending = True)

	mean_weighted_rank
H9	6.000000
CTL09A	9.000000
CTL01A	10.500000
KTD8_2	12.000000
CTL06F	16.250000
CTL05A	26.500000
CTL08A	28.333333
UCSFi001-A	30.750000
CTL04E	36.500000
CTL07C	41.333333
H1	46.750000
CTL02A	51.666667

Weighted mean - including percentages

all_results_df[(all_results_df['Timepoint'] == 'day 12')].groupby('DONOR').sum('REPRESENTATION').sort_values(by = 'REPRESENTATION')

	REPRESENTATION	MIX ID	Timepoint_int
DONOR
CTL02A	0.03948	36	108
KTD8_2	0.06124	12	36
H1	0.06882	48	144
CTL07C	0.11447	36	108
CTL08A	0.26965	33	108
CTL05A	0.50789	48	144
CTL04E	0.57805	45	144
UCSFi001-A	0.68447	45	144
H9	1.39503	9	36
CTL06F	3.15496	45	144
CTL09A	4.78777	33	108
CTL01A	6.33816	51	144

mean_d_df['mean_Representation'] = all_results_df[(all_results_df['Timepoint'] == 'day 12')].groupby('DONOR').mean('REPRESENTATION')['REPRESENTATION']
mean_d_df['combined_scores'] = mean_d_df['mean_weighted_rank'] * mean_d_df['mean_Representation']
mean_d_df.sort_values('combined_scores')

	mean_weighted_rank	mean_Representation	combined_scores
CTL02A	51.666667	0.004387	0.226644
KTD8_2	12.000000	0.020413	0.244960
H1	46.750000	0.005735	0.268111
CTL07C	41.333333	0.012719	0.525714
CTL08A	28.333333	0.029961	0.848898
CTL05A	26.500000	0.042324	1.121590
UCSFi001-A	30.750000	0.057039	1.753954
CTL04E	36.500000	0.048171	1.758235
H9	6.000000	0.465010	2.790060
CTL06F	16.250000	0.262913	4.272342
CTL09A	9.000000	0.531974	4.787770
CTL01A	10.500000	0.528180	5.545890

mean_d_df.sort_values('combined_scores').to_csv('../../data/csv/CensusSeq_combined_weighted_rank_d12.csv')

np.round(mean_d_df.sort_values('combined_scores'), 2)

	mean_weighted_rank	mean_Representation	combined_scores
CTL02A	51.67	0.00	0.23
KTD8_2	12.00	0.02	0.24
H1	46.75	0.01	0.27
CTL07C	41.33	0.01	0.53
CTL08A	28.33	0.03	0.85
CTL05A	26.50	0.04	1.12
UCSFi001-A	30.75	0.06	1.75
CTL04E	36.50	0.05	1.76
H9	6.00	0.47	2.79
CTL06F	16.25	0.26	4.27
CTL09A	9.00	0.53	4.79
CTL01A	10.50	0.53	5.55

all_timepoints_norm_weighted['day 12'] = mean_d_df['combined_scores']

Day 25

ranked_lists = {}
for m in all_results_df[all_results_df['Timepoint'] == 'day 25']['MIX ID'].unique():
    df = all_results_df[((all_results_df['Timepoint'] == 'day 25') & (all_results_df['MIX ID'] == m)  )].sort_values(by = 'REPRESENTATION', ascending = False)
    result = df.groupby('DONOR').sum().sort_values(by = 'REPRESENTATION', ascending = False).reset_index()
    ranked_lists[f'Mix {m}'] = result.DONOR.tolist()

data = list(ranked_lists.values())

Weighted mean - not including percentages

d = defaultdict(list)
d

for l in data:
    for idx, value in enumerate(l):
        d[value].append( (idx + 1) * len(l)) # or (idx + 1) * len(l)? which one is better?

mean_d = {}
for l in d:
    mean_d[l] = np.mean(d[l])

mean_d_df = pd.DataFrame(mean_d.values(), mean_d.keys())
mean_d_df.columns = ['mean_weighted_rank']
mean_d_df.sort_values(by = 'mean_weighted_rank', ascending = True)

	mean_weighted_rank
CTL09A	7.000000
CTL01A	9.200000
H9	12.000000
CTL06F	16.250000
CTL05A	22.500000
KTD8_2	24.000000
CTL08A	24.250000
CTL04E	27.200000
CTL02A	37.750000
UCSFi001-A	37.750000
CTL07C	42.666667
H1	49.250000

Weighted mean - including percentages

all_results_df[(all_results_df['Timepoint'] == 'day 25')].groupby('DONOR').sum('REPRESENTATION').sort_values(by = 'REPRESENTATION')

	REPRESENTATION	MIX ID	Timepoint_int
DONOR
KTD8_2	0.02428	12	75
H1	0.05101	48	300
CTL07C	0.08319	36	225
CTL02A	0.09660	84	375
UCSFi001-A	0.31279	45	300
CTL08A	0.35909	81	375
CTL05A	0.51931	48	300
CTL04E	0.82397	93	450
H9	1.07928	9	75
CTL06F	3.85586	45	300
CTL09A	5.24697	33	225
CTL01A	11.54762	99	450

mean_d_df['mean_Representation'] = all_results_df[(all_results_df['Timepoint'] == 'day 25')].groupby('DONOR').mean('REPRESENTATION')['REPRESENTATION']
mean_d_df['combined_scores'] = mean_d_df['mean_weighted_rank'] * mean_d_df['mean_Representation']
mean_d_df.sort_values('combined_scores')

	mean_weighted_rank	mean_Representation	combined_scores
KTD8_2	24.000000	0.008093	0.194240
H1	49.250000	0.004251	0.209354
CTL02A	37.750000	0.006440	0.243110
CTL07C	42.666667	0.009243	0.394382
CTL08A	24.250000	0.023939	0.580529
CTL05A	22.500000	0.043276	0.973706
UCSFi001-A	37.750000	0.026066	0.983985
CTL04E	27.200000	0.045776	1.245110
CTL09A	7.000000	0.582997	4.080977
H9	12.000000	0.359760	4.317120
CTL06F	16.250000	0.321322	5.221477
CTL01A	9.200000	0.641534	5.902117

np.round(mean_d_df.sort_values('combined_scores'), 2)

	mean_weighted_rank	mean_Representation	combined_scores
KTD8_2	24.00	0.01	0.19
H1	49.25	0.00	0.21
CTL02A	37.75	0.01	0.24
CTL07C	42.67	0.01	0.39
CTL08A	24.25	0.02	0.58
CTL05A	22.50	0.04	0.97
UCSFi001-A	37.75	0.03	0.98
CTL04E	27.20	0.05	1.25
CTL09A	7.00	0.58	4.08
H9	12.00	0.36	4.32
CTL06F	16.25	0.32	5.22
CTL01A	9.20	0.64	5.90

mean_d_df.sort_values('combined_scores').to_csv('../../data/csv/CensusSeq_combined_weighted_rank_d25.csv')

all_timepoints_norm_weighted['day 25'] = mean_d_df['combined_scores']

Day 50

ranked_lists = {}
for m in all_results_df[all_results_df['Timepoint'] == 'day 50']['MIX ID'].unique():
    df = all_results_df[((all_results_df['Timepoint'] == 'day 50') & (all_results_df['MIX ID'] == m)  )].sort_values(by = 'REPRESENTATION', ascending = False)
    result = df.groupby('DONOR').sum().sort_values(by = 'REPRESENTATION', ascending = False).reset_index()
    ranked_lists[f'Mix {m}'] = result.DONOR.tolist()

data = list(ranked_lists.values())

Weighted mean - not including percentages

d = defaultdict(list)
d

for l in data:
    for idx, value in enumerate(l):
        d[value].append( (idx + 1) * len(l)) # or (idx + 1) * len(l)? which one is better?

mean_d = {}
for l in d:
    mean_d[l] = np.mean(d[l])

mean_d_df = pd.DataFrame(mean_d.values(), mean_d.keys())
mean_d_df.columns = ['mean_weighted_rank']
mean_d_df.sort_values(by = 'mean_weighted_rank', ascending = True)

	mean_weighted_rank
H9	6.000000
CTL09A	9.000000
CTL01A	9.200000
KTD8_2	18.000000
CTL06F	19.000000
CTL04E	20.600000
CTL08A	23.250000
CTL05A	35.000000
UCSFi001-A	36.250000
CTL02A	38.250000
CTL07C	41.333333
H1	46.750000

Weighted mean - including percentages

all_results_df[(all_results_df['Timepoint'] == 'day 50')].groupby('DONOR').sum('REPRESENTATION').sort_values(by = 'REPRESENTATION')

	REPRESENTATION	MIX ID	Timepoint_int
DONOR
KTD8_2	0.01621	12	150
CTL02A	0.04337	60	600
H1	0.05719	48	600
CTL07C	0.08289	36	450
CTL05A	0.18000	48	600
CTL08A	0.29861	57	600
UCSFi001-A	0.34944	45	600
CTL06F	1.42171	45	600
H9	2.05306	9	150
CTL04E	2.20840	69	750
CTL09A	3.99528	33	450
CTL01A	10.29388	75	750

mean_d_df['mean_Representation'] = all_results_df[(all_results_df['Timepoint'] == 'day 50')].groupby('DONOR').mean('REPRESENTATION')['REPRESENTATION']
mean_d_df['combined_scores'] = mean_d_df['mean_weighted_rank'] * mean_d_df['mean_Representation']
mean_d_df.sort_values('combined_scores')

	mean_weighted_rank	mean_Representation	combined_scores
KTD8_2	18.000000	0.005403	0.097260
CTL02A	38.250000	0.003614	0.138242
H1	46.750000	0.004766	0.222803
CTL07C	41.333333	0.009210	0.380680
CTL05A	35.000000	0.015000	0.525000
CTL08A	23.250000	0.024884	0.578557
UCSFi001-A	36.250000	0.029120	1.055600
CTL06F	19.000000	0.118476	2.251041
CTL04E	20.600000	0.147227	3.032869
CTL09A	9.000000	0.443920	3.995280
H9	6.000000	0.684353	4.106120
CTL01A	9.200000	0.686259	6.313580

mean_d_df.sort_values('combined_scores').to_csv('../../data/csv/CensusSeq_combined_weighted_rank_d50.csv')

np.round(mean_d_df.sort_values('combined_scores'), 2)

	mean_weighted_rank	mean_Representation	combined_scores
KTD8_2	18.00	0.01	0.10
CTL02A	38.25	0.00	0.14
H1	46.75	0.00	0.22
CTL07C	41.33	0.01	0.38
CTL05A	35.00	0.02	0.52
CTL08A	23.25	0.02	0.58
UCSFi001-A	36.25	0.03	1.06
CTL06F	19.00	0.12	2.25
CTL04E	20.60	0.15	3.03
CTL09A	9.00	0.44	4.00
H9	6.00	0.68	4.11
CTL01A	9.20	0.69	6.31

all_timepoints_norm_weighted['day 50'] = mean_d_df['combined_scores']

all_results_df[((all_results_df['MIX ID'] == 6) & (all_results_df['Timepoint'] == 'day 50'))].groupby(['DONOR']).mean()

	REPRESENTATION	MIX ID	Timepoint_int
DONOR
CTL01A	0.359893	6	50
CTL02A	0.001277	6	50
CTL04E	0.029830	6	50
CTL05A	0.018230	6	50
CTL06F	0.096497	6	50
CTL07C	0.010813	6	50
CTL08A	0.038487	6	50
CTL09A	0.423120	6	50
H1	0.007097	6	50
UCSFi001-A	0.014763	6	50