<jupyter_start><jupyter_text>*Running the nb with 08 threshold for replicate correlation - 2021-02-16* Table of Contents1&nbsp;&nbsp;Setup1.1&nbsp;&nbsp;QC Thresholds1.2&nbsp;&nbsp;Inputs1.3&nbsp;&nbsp;Load expression data1.4&nbsp;&nbsp;Load QC data1.4.1&nbsp;&nbsp;Note that there are 6 genes missing in the log_tpm expression data. Those are:  SRX1652180, SRX1655390, SRX1655391,  SRX1655392, SRX1655393, SRX16553941.5&nbsp;&nbsp;Load metadata1.6&nbsp;&nbsp;Remove extra sample rows2&nbsp;&nbsp;Check statistics2.1&nbsp;&nbsp;FastQC stats2.2&nbsp;&nbsp;Number of aligned reads3&nbsp;&nbsp;Examine Global Correlations3.1&nbsp;&nbsp;Hierarchical Clustering4&nbsp;&nbsp;Remove failed samples5&nbsp;&nbsp;Save Metadata for Curation Process6&nbsp;&nbsp;Metadata Curation6.1&nbsp;&nbsp;BioProject counts of passing metadata6.2&nbsp;&nbsp;BioProject counts of all metadata7&nbsp;&nbsp;Correlations between replicates7.1&nbsp;&nbsp;Compute Pearson R Score7.2&nbsp;&nbsp;Investigating Individual Correlation7.3&nbsp;&nbsp;Investigating rpaA7.4&nbsp;&nbsp;Investigating ppGpp with theo7.5&nbsp;&nbsp;Clock7.6&nbsp;&nbsp;S77.7&nbsp;&nbsp;Drop Samples with no replicates7.8&nbsp;&nbsp;Calculating failed samples for Replicate correlation7.9&nbsp;&nbsp;Calculating total failed samples (FastQC to Replicate correlation)7.10&nbsp;&nbsp;Calculating stats with adding failed time-series data7.10.0.1&nbsp;&nbsp;1. Adding the failed time-series data back to metadata_curated DataFrame7.10.0.2&nbsp;&nbsp;1. Changing "false" to "true" for these projects t ocalculate finale pie chart8&nbsp;&nbsp;Manually add failed time-series samples to curated metadata and save9&nbsp;&nbsp;Save final datasets10&nbsp;&nbsp;Plots This is a template notebook for performing preliminary quality control on your organism's expression data.<jupyter_code>import pandas as pd
import os
from os import path
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
#sns.set_style('ticks')
sns.set_style('whitegrid')<jupyter_output><empty_output><jupyter_text>Setup  QC Thresholds<jupyter_code>min_mrna_reads = 5e5 # Minimum number of reads mapped to mRNA<jupyter_output><empty_output><jupyter_text>Inputs Show files in the data directory<jupyter_code>data_dir = path.join('/home/tahani/Documents/elongatus/data/')
os.listdir(data_dir)<jupyter_output><empty_output><jupyter_text>Enter name of curated metadata file here<jupyter_code>metadata_file = "0_elongatus_2020-08-20.tsv" # Enter metadata filename here<jupyter_output><empty_output><jupyter_text>Load expression data<jupyter_code>DF_log_tpm = pd.read_csv(path.join(data_dir,'0_log_tpm.csv'),index_col=0).fillna(0)
print('Number of genes:',DF_log_tpm.shape[0])
print('Number of samples:',DF_log_tpm.shape[1])

#note that the old log_tpm file contained 2714 genes (including the pseudogenes with rockhopper method)
DF_log_tpm.head()<jupyter_output><empty_output><jupyter_text>Load QC dataThere may be some datasets that failed along the processing pipeline, so the number of samples with QC data may be higher than the number of samples with expression data.<jupyter_code>DF_qc_stats = pd.read_csv(path.join(data_dir,'0_multiqc_stats.tsv'),index_col=0, sep='\t')
print('Number of samples with QC data:',DF_qc_stats.shape[0])<jupyter_output>Number of samples with QC data: 323<jupyter_text>Note that there are 6 genes missing in the log_tpm expression data. Those are:  SRX1652180, SRX1655390, SRX1655391,  SRX1655392, SRX1655393, SRX1655394<jupyter_code>DF_qc_stats.fillna(0,inplace=True)
DF_qc_stats.head()<jupyter_output><empty_output><jupyter_text>Load metadata<jupyter_code>DF_metadata = pd.read_csv(path.join(data_dir,metadata_file),index_col=0,sep='\t')
print('Number of samples with metadata:',DF_metadata.shape[0])
DF_metadata.head()<jupyter_output>Number of samples with metadata: 323<jupyter_text>Remove extra sample rows<jupyter_code># Ensure that metadata and qc_stats data contain all log_tpm sample information
assert(set(DF_log_tpm.columns) - set(DF_metadata.index) == set())
assert(set(DF_log_tpm.columns) - set(DF_qc_stats.index) == set())
DF_metadata = DF_metadata.loc[DF_log_tpm.columns]
DF_qc_stats = DF_qc_stats.loc[DF_log_tpm.columns]
print('Number of samples with QC data:',DF_qc_stats.shape[0])
print('Number of samples in log_tpm:',DF_log_tpm.shape[1])<jupyter_output>Number of samples with QC data: 317
Number of samples in log_tpm: 317<jupyter_text>Check statistics From here, create a new spreadsheet where you can flag samples based on various QC statistics  FastQC stats<jupyter_code>fastqc_cols = ['per_base_sequence_quality',
       'per_tile_sequence_quality', 'per_sequence_quality_scores',
       'per_base_sequence_content', 'per_sequence_gc_content',
       'per_base_n_content', 'sequence_length_distribution',
       'sequence_duplication_levels', 'overrepresented_sequences',
       'adapter_content']
DF_fastqc = DF_qc_stats[fastqc_cols]
ax = sns.heatmap(DF_fastqc.replace('pass',1).replace('warn',0).replace('fail',-1),
            cmap='RdYlBu',vmax=1.3,vmin=-1.3)
cbar = ax.collections[0].colorbar
cbar.set_ticks([-1,0,1])
cbar.set_ticklabels(['fail','warn','pass'])<jupyter_output><empty_output><jupyter_text>The following four categories are the most important:    - per_base_sequence_quality    - per_sequence_quality_scores    - per_base_n_content    - adapter_content    If a sample does not pass any of these four categories, discard the sample<jupyter_code>fastqc_fail_cols = ['per_base_sequence_quality','per_sequence_quality_scores','per_base_n_content','adapter_content']
DF_failed_fastqc = DF_fastqc[fastqc_fail_cols][(DF_fastqc[fastqc_fail_cols] != 'pass').any(axis=1)]
DF_failed_fastqc[fastqc_fail_cols]<jupyter_output><empty_output><jupyter_text>Mark samples as failed.<jupyter_code>DF_metadata['passed_fastqc'] = ~DF_metadata.index.isin(DF_failed_fastqc.index)
print('#samples that failed FastQC', len([i for i in DF_metadata['passed_fastqc'] if i == False]))<jupyter_output>#samples that failed FastQC 3<jupyter_text>Number of aligned reads The following histogram shows how many reads map to coding sequences (i.e. mRNA). Too few aligned reads reduces the sensitivity of the resulting data.<jupyter_code>min_mrna_reads
fig,ax = plt.subplots()
ax.hist(DF_qc_stats['Assigned']/1e6,bins=50,alpha=0.8)
ymin,ymax = ax.get_ylim()
ax.vlines(min_mrna_reads/1e6,ymin,ymax,color='r')
ax.set_ylim((ymin,ymax))
ax.set_xlabel('# Reads (M)',fontsize=14)
ax.set_ylabel('# Samples',fontsize=14)
ax.set_title('Number of reads mapped to CDS',fontsize=16)<jupyter_output><empty_output><jupyter_text>Identify samples with poor read depth:<jupyter_code>DF_failed_mrna = DF_qc_stats[DF_qc_stats['Assigned'] < min_mrna_reads].sort_values('Assigned')
DF_failed_mrna

#Note: 37 were failed in the previous data and now its 44<jupyter_output><empty_output><jupyter_text>Mark samples as failed.<jupyter_code>DF_metadata['passed_reads_mapped_to_CDS'] = ~DF_metadata.index.isin(DF_failed_mrna.index)
print('#samples that failed number of aligned reads', len([i for i in DF_metadata['passed_reads_mapped_to_CDS'] if i == False]))<jupyter_output>#samples that failed number of aligned reads 44<jupyter_text>Examine Global Correlations  Hierarchical Clustering  A clustermap is a great way to visualize the global correlations between one sample and all others. The following code uses hierarchical clustering to identify specific clusters in the clustermap      To increase the number of clusters, decrease the value of `thresh`. To decrease the number of clusters, increase the value of `thresh`<jupyter_code>import scipy.cluster.hierarchy as sch
import matplotlib.patches as patches
# change this to get different number of clusters
thresh = .3

# retrieve clusters using fcluster 
corr = DF_log_tpm.corr()
corr.fillna(0,inplace=True)
dist = sch.distance.pdist(corr)
link = sch.linkage(dist, method='complete')
clst = pd.DataFrame(index=DF_log_tpm.columns)
clst['cluster'] = sch.fcluster(link, thresh * dist.max(), 'distance')

#get colors for each cluster
cm = plt.cm.get_cmap('tab20')
clr = dict(zip(clst.cluster.unique(), cm.colors))
clst['color'] = clst.cluster.map(clr)

print('Number of cluster: ', len(clr))<jupyter_output>Number of cluster:  5<jupyter_text>To view sample IDs in the clustermap, set `xticklabels` and `yticklabels` to `True`. You can increase the `size` variable to improve readability of sample IDs<jupyter_code>size = 12

legend_TN = [patches.Patch(color=c, label=l) for l,c in clr.items()]

sns.set(rc={'figure.facecolor':'white'})
g = sns.clustermap(DF_log_tpm.corr(), figsize=(size,size), 
                   row_linkage=link, col_linkage=link, col_colors=clst.color,
                  yticklabels=False, xticklabels=False, vmin=0, vmax=1)

l2=g.ax_heatmap.legend(loc='upper left', bbox_to_anchor=(1.01,0.85), handles=legend_TN,frameon=True)
l2.set_title(title='Clusters',prop={'size':10})<jupyter_output><empty_output><jupyter_text>Select clusters to remove.<jupyter_code>#indicate which clusters you want to remove
remove_clst = [5,4]
failed_global_corr = clst[clst.cluster.isin(remove_clst)].index
failed_global_corr 
DF_metadata['passed_global_correlation'] = ~DF_metadata.index.isin(failed_global_corr)
print('#samples that failed global_correlation ', len([i for i in DF_metadata['passed_global_correlation'] if i == False]))
DF_metadata.head()<jupyter_output><empty_output><jupyter_text>Remove failed samples<jupyter_code>qc_columns = ['passed_fastqc',
              'passed_reads_mapped_to_CDS',
              'passed_global_correlation']
pass_qc = DF_metadata[qc_columns].all(axis=1)
DF_metadata_passed = DF_metadata[pass_qc]
DF_metadata_passed.shape
_,_,pcts = plt.pie(pass_qc.value_counts().sort_values(),
        labels = ['Failed','Passed'],
        colors=['tab:red','tab:blue'],
        autopct='%.0f%%',textprops={'size':16});

# Colors percents white
for pct in pcts:
    pct.set_color('white')<jupyter_output><empty_output><jupyter_text>a total of 63 samples failed out of 317 samples in total gives a percentage of failed sames to be 17%<jupyter_code>DF_metadata_passed.head()<jupyter_output><empty_output><jupyter_text>Save Metadata for Curation Process<jupyter_code># DF_metadata_passed.to_csv(os.path.join(data_dir,'BRUH.csv'))<jupyter_output><empty_output><jupyter_text>Metadata Curation The following sections can only be run after metadata curation is complete. To enable metadata curation, the code in this section sorts data by BioProject.   BioProject counts of passing metadata<jupyter_code>DF_metadata_passed.BioProject.value_counts().sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>BioProject counts of all metadata<jupyter_code>DF_metadata.BioProject.value_counts().sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Correlations between replicates  First, get a full sample name<jupyter_code># Turn off pesky warning
pd.set_option('mode.chained_assignment', None)
curated_metadata_file= "2_metadata_passed_curated_262_INPUT.csv" #gives 81 
# curated_metadata_file= "2_metadata_passed_262_UPDATED.csv" #gives 81 

DF_metadata_curated = pd.read_csv(path.join(data_dir,curated_metadata_file),index_col=0,sep=',')
print('Number of samples with metadata:',DF_metadata_passed.shape[0])
DF_metadata_curated.head()
DF_metadata_curated['full_name'] = DF_metadata_curated.loc[:,'project'].str.cat(DF_metadata_curated.loc[:,'condition'],sep=':')
DF_metadata_curated['full_name']<jupyter_output><empty_output><jupyter_text>Compute Pearson R Score Once you have updated your metadata files with the sample information, we can investigate correlations between biological replicates. We require biological replicates to have a Pearson R correlation above 0.95. For samples with more than 2 replicates, the replicates must have R >= 0.95 with at least one other replicate or it will be dropped.<jupyter_code>from tqdm.notebook import tqdm
import itertools
from scipy import stats
import numpy as np
rep_corrs = {}
rand_corrs = {}

num_comparisons = len(DF_metadata_curated)*(len(DF_metadata_curated)-1)/2
for exp1,exp2 in tqdm(itertools.combinations(DF_metadata_curated.index,2),total=num_comparisons):
    if DF_metadata_curated.loc[exp1,'full_name'] == DF_metadata_curated.loc[exp2,'full_name']:
        rep_corrs[(exp1,exp2)] = stats.pearsonr(DF_log_tpm[exp1],DF_log_tpm[exp2])[0]
    else:
        rand_corrs[(exp1,exp2)] = stats.pearsonr(DF_log_tpm[exp1],DF_log_tpm[exp2])[0]
sns.set_style('ticks')
fig,ax = plt.subplots(figsize=(5,5))
ax2 = ax.twinx()
ax2.hist(rep_corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='green',linewidth=0)
ax.hist(rand_corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='blue',linewidth=0)
ax.set_title('Pearson R correlation between experiments',fontsize=14)
ax.set_xlabel('Pearson R correlation',fontsize=14)
ax.set_ylabel('Different Conditions',fontsize=14)
ax2.set_ylabel('Known Replicates',fontsize=14)

med_corr = np.median([v for k,v in rep_corrs.items()])
print('Median Pearson R between replicates: {:.2f}'.format(med_corr))
fig.savefig('myimage.png', format='png', bbox_inches='tight', dpi=1200)
#pearson r cutoff for replicates
rcutoff = 0.9

#for each sample get max correlation between replicates
dissimilar = []
for idx, grp in DF_metadata_curated.groupby('full_name'):
    ident = np.identity(len(grp))
    corrs = (DF_log_tpm[grp.index].corr() - ident).max()
    dissimilar.extend(corrs[corrs<rcutoff].index)
DF_metadata['passed_replicate_corr'] = ~DF_metadata.index.isin(dissimilar)
DF_metadata_curated['passed_similar_replicates'] = ~DF_metadata_curated.index.isin(dissimilar)

print('#samples that failed replicates', len([i for i in DF_metadata_curated['passed_similar_replicates']  if i == False]))
# marking all 85 samples with a false 
x = DF_metadata_curated[DF_metadata_curated['passed_similar_replicates'] == False]
x.to_csv(os.path.join(data_dir,'3_samples_failed_rep_corr.csv'))<jupyter_output><empty_output><jupyter_text>Investigating Individual Correlation<jupyter_code># rpaA = DF_metadata_curated[DF_metadata_curated.BioProject=='PRJNA354335']
# rpaA_corrs = {}
# #rpaA = DF_metadata_passed[DF_metadata_passed.BioProject=='PRJNA280002'].index
# num_comparisons = len(rpaA)*(len(rpaA)-1)/2
# for exp1,exp2 in tqdm(itertools.combinations(rpaA.index,2),total=num_comparisons):
#     rpaA_corrs[(exp1,exp2)] = stats.pearsonr(DF_log_tpm[exp1],DF_log_tpm[exp2])[0]
# fig,ax = plt.subplots(figsize=(5,5))
# #ax2 = ax.twinx()
# #ax2.hist(rep_corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='green',linewidth=0)
# ax.hist(rpaA_corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='blue',linewidth=0)
# ax.set_title('Pearson R correlation between experiments',fontsize=14)
# ax.set_xlabel('Pearson R correlation',fontsize=14)
# ax.set_ylabel('Different Conditions',fontsize=14)
# #ax2.set_ylabel('Known Replicates',fontsize=14)

# med_corr = np.median([v for k,v in rpaA_corrs.items()])
# print('Median Pearson R between replicates: {:85.2f}'.format(med_corr))<jupyter_output><empty_output><jupyter_text>Cannot tell clearly establish a threshold since there is no clear seperation. So we create lesser combinations as shown in the section below.   Investigating rpaA<jupyter_code>wt_dusk = ['SRX2356912','SRX2356914','SRX2356915']
wt_dawn = ['SRX2356921', 'SRX2356922','SRX2356923']
ko_dusk = ['SRX2356924', 'SRX2356925','SRX2356926', 'SRX2356927', 'SRX2356928', 'SRX2356929']
ko_dawn = ['SRX2356933', 'SRX2356934' ,'SRX2356935']
corrs = {}
for i in range(0,2):
    samp1 = ko_dawn[i] 
    samp2 = ko_dawn[i+1]
    corrs[(samp1,samp2)]= stats.pearsonr(DF_log_tpm[samp1],DF_log_tpm[samp2])[0]
    #corrs.append(stats.pearsonr(DF_log_tpm.sample1, DF_log_tpm.sample2))
    #print(sample1,sample2)
corrs<jupyter_output><empty_output><jupyter_text>Investigating ppGpp with theo<jupyter_code>wt_dusk = ['SRX3172761','SRX3172762','SRX3172763']
ko_dusk = ['SRX3172764','SRX3172765','SRX3172766']
corrs = {}
for i in range(0,2):
    samp1 = ko_dusk[i] 
    samp2 = ko_dusk[i+1]
    corrs[(samp1,samp2)]= stats.pearsonr(DF_log_tpm[samp1],DF_log_tpm[samp2])[0]
    #corrs.append(stats.pearsonr(DF_log_tpm.sample1, DF_log_tpm.sample2))
    #print(sample1,sample2)
corrs<jupyter_output><empty_output><jupyter_text>Clock<jupyter_code>wt = ['SRX359366','SRX359367','SRX359368','SRX359369','SRX359370','SRX359371']
KaiBC_IPTG = ['SRX359372','SRX359373','SRX359374','SRX359375','SRX359376','SRX359377']
KaiBC = ['SRX359378','SRX359379','SRX359380']
empty = ['SRX359381','SRX359382','SRX359383']

corrs = {}
for i in range(0,2):
    samp1 = empty[i] 
    samp2 = empty[i+1]
    corrs[(samp1,samp2)]= stats.pearsonr(DF_log_tpm[samp1],DF_log_tpm[samp2])[0]
    #corrs.append(stats.pearsonr(DF_log_tpm.sample1, DF_log_tpm.sample2))
    #print(sample1,sample2)
corrs
# fig,ax = plt.subplots(figsize=(5,5))
# #ax2 = ax.twinx()
# #ax2.hist(rep_corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='green',linewidth=0)
# ax.hist(corrs.values(),bins=50,range=(0.2,1),alpha=0.8,color='blue',linewidth=0)
# ax.set_title('Pearson R correlation between experiments',fontsize=14)
# ax.set_xlabel('Pearson R correlation',fontsize=14)
# ax.set_ylabel('Different Conditions',fontsize=14)
# #ax2.set_ylabel('Known Replicates',fontsize=14)

# med_corr = np.median([v for k,v in rpaA_corrs.items()])
# print('Median Pearson R between replicates: {:85.2f}'.format(med_corr))<jupyter_output><empty_output><jupyter_text>S7<jupyter_code>ph = ['SRX359366','SRX359367','SRX359368','SRX359369','SRX359370','SRX359371']
saline = ['SRX359372','SRX359373','SRX359374','SRX359375','SRX359376','SRX359377']
co2 = ['SRX359378','SRX359379','SRX359380']
temp = ['SRX359381','SRX359382','SRX359383']
light = ['SRX359381','SRX359382','SRX359383']
U = ['SRX359381','SRX359382','SRX359383']

corrs = {}
for i in range(0,2):
    samp1 = empty[i] 
    samp2 = empty[i+1]
    corrs[(samp1,samp2)]= stats.pearsonr(DF_log_tpm[samp1],DF_log_tpm[samp2])[0]
    #corrs.append(stats.pearsonr(DF_log_tpm.sample1, DF_log_tpm.sample2))
    #print(sample1,sample2)
corrs<jupyter_output><empty_output><jupyter_text>between samples 15-21, the correlaitons is 0.55 and between 14-15 or 21-22 the correlation is way higher so in this case we passed all the samples from this project (Tahani/anand 9/11/2020)  Drop Samples with no replicates cond_counts = DF_metadata_curated.condition.value_counts()drop_conds = cond_counts[cond_counts < 2].indexdrop_conds DF_metadata_curated['passed_number_replicates'] = ~DF_metadata_curated.condition.isin(drop_conds)print('samples that failed replicates', len([i for i in DF_metadata_curated['passed_number_replicates']  if i == False])) DF_metadata_curated = DF_metadata_curated[DF_metadata_curated['passed_number_replicates']]  Calculating failed samples for Replicate correlation<jupyter_code>DF_metadata['passed_similar_replicates'] = ~DF_metadata.index.isin(dissimilar)
print('#samples that failed similar replicates', len([i for i in DF_metadata['passed_similar_replicates']  if i == False]))
qc_columns2 = ['passed_similar_replicates']
pass_qc2 = DF_metadata[qc_columns2].all(axis=1)
# pass_qc2 = DF_metadata_curated[qc_columns2].all(axis=1)
# DF_metadata_curated = DF_metadata_curated[pass_qc2]
# DF_metadata_curated.shape
# # 262 - 8 failed samples 
_,_,pcts = plt.pie(pass_qc2.value_counts().sort_values(),
        labels = ['Failed','Passed'],
        colors=['tab:red','tab:blue'],
        autopct='%.0f%%',textprops={'size':16});

# Colors percents white
for pct in pcts:
    pct.set_color('white')<jupyter_output><empty_output><jupyter_text>Total samples failing replicate correlation is 26%  (81 failed samples from the original 317)  Calculating total failed samples (FastQC to Replicate correlation)<jupyter_code>qc_columns3 = ['passed_fastqc',
              'passed_reads_mapped_to_CDS',
              'passed_global_correlation','passed_similar_replicates']
pass_qc3 = DF_metadata[qc_columns3].all(axis=1)
_,_,pcts = plt.pie(pass_qc3.value_counts().sort_values(),
        labels = ['Failed','Passed'],
        colors=['tab:red','tab:blue'],
        autopct='%.0f%%',textprops={'size':16});

# Colors percents white
for pct in pcts:
    pct.set_color('white')
pass_qc3.value_counts().sort_values()<jupyter_output><empty_output><jupyter_text>Total samples failing ALL METRICS is 136 samples (34%) of the total 317 samples  Calculating stats with adding failed time-series data  1. Adding the failed time-series data back to metadata_curated DataFrame<jupyter_code>#before 
DF_metadata_curated.shape
rpaA = DF_metadata_curated[DF_metadata_curated.BioProject=='PRJNA354335']
ppGpp_theo = DF_metadata_curated[DF_metadata_curated.BioProject=='PRJNA403840']
clock = DF_metadata_curated[DF_metadata_curated.BioProject=='PRJNA221220']
# After
DF_metadata_curated.shape<jupyter_output><empty_output><jupyter_text>1. Changing "false" to "true" for these projects t ocalculate finale pie chart<jupyter_code>DF_metadata.to_csv(os.path.join(data_dir,'temp.csv'))
DF_metadata = pd.read_csv(path.join(data_dir,'3_after_adding_True_to_failed_samples.csv'),index_col=0)
qc_columns4 = ['passed_fastqc',
              'passed_reads_mapped_to_CDS',
              'passed_global_correlation','passed_similar_replicates']
pass_qc4 = DF_metadata[qc_columns4].all(axis=1)
pass_qc4.value_counts().sort_values()
_,_,pcts = plt.pie(pass_qc4.value_counts().sort_values(),
        labels = ['Failed','Passed'],
        colors=['tab:red','tab:blue'],
        autopct='%.0f%%',textprops={'size':16});

# Colors percents white
for pct in pcts:
    pct.set_color('white')<jupyter_output><empty_output><jupyter_text>stats after reviving the failed samples from replicate correlation (ppGpp_theo, rpaA and clock projects)  Manually add failed time-series samples to curated metadata and save  1. get updated metadata with all pased replicate correlation samples<jupyter_code>qc_columns = ['passed_similar_replicates']
DF_metadata_final = DF_metadata_curated[DF_metadata_curated[qc_columns].all(axis=1)]
DF_metadata_final.shape<jupyter_output><empty_output><jupyter_text>2. Add manually the failed time-series samples<jupyter_code>DF_metadata_final = rpaA.append(DF_metadata_final) 
DF_metadata_final = ppGpp_theo.append(DF_metadata_final)
DF_metadata_final = clock.append(DF_metadata_final)
DF_metadata_final.shape
DF_log_tpm_final = DF_log_tpm[DF_metadata_final.index]
DF_log_tpm_final.shape
DF_metadata_final.BioProject.value_counts().sort_values(ascending=False)<jupyter_output><empty_output><jupyter_text>Save final datasets<jupyter_code>DF_log_tpm_final.to_csv(os.path.join(data_dir,'3_log_tpm_qc_ICA.csv'))
DF_metadata_final.to_csv(os.path.join(data_dir,'3_metadata_qc_ICA.csv'))<jupyter_output><empty_output><jupyter_text>Plots<jupyter_code>DF_metadata_final.BioProject.value_counts().sort_values(ascending=False)
data = DF_metadata_final.project.value_counts().sort_values(ascending=False)
data
# Libraries
import matplotlib.pyplot as plt
 
# Make data: I have 3 groups and 7 subgroups
group_names = ['light', 'ppGpp', 'sigFctr', 'FFA', 'rpaA']
group_size=[60,40,35,15,15]
subgroup_names=['PRJNA412032','PRJNA401777','PRJNA404081','PRJNA401742','PRJNA415380','PRJNA472248', 'PRJNA196229', 'PRJNA354335'  ]
subgroup_size=[60,8,12,8,12,35,15,15]
 
# Create colors
a, b, c, d, e =[plt.cm.Blues, plt.cm.Reds, plt.cm.Greens, plt.cm.Purples, plt.cm.Wistia]
 
# First Ring (outside)
fig, ax = plt.subplots()
ax.axis('equal')
mypie, _ = ax.pie(group_size, radius=1.7, labels=group_names, colors=[a(0.6), b(0.6), c(0.6), d(0.6), e(0.6)], labeldistance = 1.1)
plt.setp( mypie, width=0.8, edgecolor='white')
 
# Second Ring (Inside)
mypie2, _ = ax.pie(subgroup_size, radius=1.3-0.3, labels=subgroup_names, labeldistance=0.7, colors=[a(0.4), b(0.5), b(0.4), b(0.3), b(0.2), c(0.4), d(0.4), e(0.4)],rotatelabels= True)
plt.setp( mypie2, width=0.4, edgecolor='white')
plt.margins(0,0)
 
# show it
plt.show()
my_circle = plt.Circle((0, 0), 0.7, color='white')

d = plt.pie(group_size, labels=group_names, autopct='%1.1f%%',
            startangle=9, labeldistance=1.05)
plt.axis('equal')
plt.gca().add_artist(my_circle)
            
plt.show()
path = '/home/tahani/Documents/elongatus/data/1_iModulon/pie_cateogries.csv'
import plotly.express as px 
  
df = pd.read_csv(path, index_col = 0)
  
fig = px.sunburst(df, path=['core', 'inner', 'outer'],
                  color_discrete_sequence = ["LightSkyBlue", "DarkSalmon", "YellowGreen", "Gold", "Violet"])
#                   px.colors.qualitative.Set2)
fig.show()

fig.write_image("fig1.pdf")
fig = px.sunburst(df, path=['core', 'inner', 'outer'],
                  color_discrete_sequence = px.colors.qualitative.Set2)
fig.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13162587129194359, 'fraction_numerical': 0.04550980710001621, 'mean_word_length': 7.5565187239944525, 'pattern_counts': {'":': 0, '<': 150, '<?xml version=': 0, '>': 148, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 2, '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': '38', 'dir': 'jupyter-structured-clean-dedup', 'id': '309977', 'license': 'MIT', 'n_tokens_mistral': 9882, 'n_tokens_neox': 8973, 'n_words': 1607, 'path': 'notebooks/QC/2_expression_QC_SOP.ipynb', 'repo_name': 'talbulus/Synechococcus_elongatus', 'size': '24676'}