In [1]:
import glob
import os
import json
from io import StringIO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors
from utils import TREATMENTS, COLORS, DATASETS

genealogies = pd.concat([pd.read_csv(x) for x in DATASETS])
PATH_IMPORT = "raw/24-03-22/**/summary.json"
ignore = [
"/Breseq/WT_s/WT_a7/WT_e7/data", # Duplicated data
"/Breseq/17_s/Line_17/data",     # Duplicated data
"/Breseq/PX-10/",                # Unrelated data
"/Breseq/PE-5-WT-L-E3",          # Low coverage
"/Breseq/PE-4-WT-L-E3"           # Re-sequenced
]

EXPORT_PATH = '04_mutations'
EXPORT_TREE_PATH = '04_mutations/trees'

if not os.path.exists(EXPORT_PATH):
    os.mkdir(EXPORT_PATH)
if not os.path.exists(EXPORT_TREE_PATH):
    os.mkdir(EXPORT_TREE_PATH)
In [2]:
def read_breseq_summary(file):
    js = json.load(open(file,'r'))
    assert len(js['references']['reference'].keys()) == 1
    ref = list(js['references']['reference'].keys())[0]
    return dict(folder=file.replace('/data/summary.json','').replace('raw/24-03-22/',''),
                 name=file.split("/")[-3],
                 ref=ref,
                 coverage_ref=js['references']['reference'][ref]['coverage_average'],
                 fraction_aligned_reads=js['reads']['total_fraction_aligned_reads'],
                 n_reads=js['reads']['total_reads'],
                 read_files = ", ".join(js['reads']['read_file'].keys()),
                 n_read_files=len(js['reads']['read_file'].keys()))
        
def read_breseq_mutation_table(file):
    with open(file, 'r') as fh:
        html = fh.read()
    read = pd.read_html(StringIO(html.split('<!--Output Html_Mutation_Table_String-->')[1].split('</body>')[0]))[0]
    read.columns = [x[1] for x in read.columns]
    read['file'] = file
    return read

def parse_PE(x):
    try:
        x = x.replace('-','_')
        splitted = x.split('_')
        if x.startswith("Line_17"):
            geno = "LCS-" if "mutS" in x else "LCS+"
            return(dict(pe=x,
                        cycle=0, size="",
                        ancestor=geno, treatment="",
                        tube="", genotype=f"{geno} ancestor"))

        if "PE" in x and "Line" not in x:
            cycle = int(splitted[1])+1
            large = splitted[-1][0].upper() == splitted[-1][0] 
            tube = splitted[-1]
            geno = 'LCS-' if 'WT_' in x else 'LCS+'
            size = 'L' if large else 'S'

            try:
                name = f"{size}-{geno}_{cycle}-I-{tube}"
                parent = genealogies.query(f"name=='{name}'").parent.values[0].replace("LLCS_ROOT","LCS- ancestor")
            except Exception:
                parent = 'UNK'
            id_geno = f"{parent.replace('_',':')}→{cycle}-I-{tube}"

        elif splitted[0] in ("WT","17"):
            cycle = 11
            large = splitted[-1][0].upper() == splitted[-1][0] 
            size = 'L' if large else 'S'
            tube  = splitted[-1]
            geno  = {'WT':'LCS-',"17":"LCS+"}[splitted[0]]
            id_geno = f"{size}-{geno}: 10-II-{tube}→"
            parent = f"{size}-{geno}_10-II-{tube}"
        else:
            return {"pe":x, "microcosm":x}
        exp = size + '-' + geno
        return(dict(pe=x,
                    cycle=cycle, size=size,
                    parent=parent,
                        ancestor=geno, treatment=exp,
                        tube=tube, genotype=id_geno))

    except Exception as ex:
          return {'error':ex, "pe":x} 
In [3]:
metadata = []
data = []
for s in glob.glob(PATH_IMPORT, recursive=True):
    if not any([x in s for x in ignore]):
        metadata.append(read_breseq_summary(s))
        metadata[-1].update(parse_PE(metadata[-1]['name'].replace('org','')))
        h = s.replace("data/summary.json", 'output/index.html')
        data.append(read_breseq_mutation_table(h))
        try:
            data[-1]['genotype'] = metadata[-1]['genotype']
        except: 
            print(metadata[-1])
            raise 
    else:
        print("Ignored:", s)
metadata = pd.DataFrame(metadata)
data = pd.concat(data)
Ignored: raw/24-03-22/Breseq/PE-4-WT-L-E3org/data/summary.json
Ignored: raw/24-03-22/Breseq/PX-10/PX_17_a6/data/summary.json
Ignored: raw/24-03-22/Breseq/PX-10/PX_17_b3/data/summary.json
Ignored: raw/24-03-22/Breseq/PX-10/PX_17_c3/data/summary.json
Ignored: raw/24-03-22/Breseq/PX-10/PX_17_e5/data/summary.json
Ignored: raw/24-03-22/Breseq/17_s/Line_17/data/summary.json
Ignored: raw/24-03-22/Breseq/WT_s/WT_a7/WT_e7/data/summary.json
Ignored: raw/24-03-22/Breseq/PE-4-WT-L-E3/data/summary.json
Ignored: raw/24-03-22/Breseq/PE-5-WT-L-E3/data/summary.json
In [4]:
# We consider that two mutations are the same if 
# ["position","mutation","annotation","gene","description"] are equal. 
# Two entries with the same values in those columns but different evidence are 
# considered the same mutation.
# Example: At position 1354000 there is a Δ9 bp deletion in wspB that is sometimes detected
# with JC sometimes with MC and JC
data.fillna({"annotation":''}, inplace=True)
cols = ["position","mutation","annotation","gene","description"]
count = data.groupby(cols).count().sort_values('genotype').genotype
count.name = "time_sequenced"

mutations = pd.pivot(data, index=cols, columns='genotype', values='evidence').fillna(0)
mutations['time_sequenced'] = count
mutations = mutations[sorted(mutations.columns, reverse=True)].reset_index()
mutations.columns.name = ''
mutations['mut_id'] = [f"MUT_{x:04}" for x in range(mutations.shape[0])]
mutations['name'] = mutations['mut_id']

# Checks that the mutations are unique
assert mutations.shape[0] == mutations.loc[:,cols].drop_duplicates().shape[0]
In [5]:
lcsp_ancestor = set(x.mut_id for k,x in mutations.iterrows() if x['LCS+ ancestor'])
lcsm_ancestor = set(x.mut_id for k,x in mutations.iterrows() if x['LCS- ancestor'])
all_sequences = set(mutations[mutations.time_sequenced==data.file.nunique()].mut_id)
any_ancestors = set.union(lcsm_ancestor, lcsp_ancestor)
both_ancestors = set.intersection(lcsm_ancestor, lcsp_ancestor)

print(f"There are:\n\t {mutations.shape[0]} unique mutations")
print(f"\t{len(lcsp_ancestor)} in LCS+ ancestor")
print(f"\t{len(lcsm_ancestor)} in LCS- ancestor")
print(f"\t{len(all_sequences)} in all sequencing")
print(f"\t{len(any_ancestors)} in any of the two ancestors")
print(f"\t{len(both_ancestors)} in any of the two ancestors")
There are:
	 9271 unique mutations
	58 in LCS+ ancestor
	74 in LCS- ancestor
	28 in all sequencing
	78 in any of the two ancestors
	54 in any of the two ancestors
In [6]:
# Count the number of mutations in metadata
cmut = data.groupby('genotype').count().mutation
metadata.set_index('genotype',inplace=True)
metadata["n_mutations"] = cmut 
In [7]:
mut_mutations = mutations[['mut' in x.gene for k,x in mutations.iterrows()]].copy()
mut_mutations['mutation_name'] = [x.gene.split()[0]+' '+x.annotation.split()[0] for k,x in mut_mutations.iterrows()]
col_mut = cols+["time_sequenced"]
mut_mutations = mut_mutations.drop(columns=col_mut).set_index('mutation_name').drop(columns=['mut_id','name']).transpose()
In [8]:
# EXPORT 
metadata.to_csv(os.path.join(EXPORT_PATH, "sequenced_cultures.csv"))
mutation_description = mutations.loc[:, ['mut_id', 'position', 'mutation',
                                         'annotation', 'gene', 'description',
                                         'time_sequenced', 'name']].set_index('mut_id')
mutation_description.to_csv(os.path.join(EXPORT_PATH,"mutation_description.csv"))

long_table = pd.merge(data, mutations.loc[:,cols+['mut_id']], left_on=cols, right_on=cols, )
long_table.to_csv(os.path.join(EXPORT_PATH,"sequencing_long_table.csv"))

wide_table = mutations.set_index('mut_id')
wide_table.to_csv(os.path.join(EXPORT_PATH,"sequencing_wide_table.csv"))

mut_mutations.to_csv(os.path.join(EXPORT_PATH,"mut_mutations_wide.csv"))

all_data = pd.merge(long_table, metadata, left_on="genotype", right_index=True)
for treat, df in all_data.groupby('treatment'):
    if treat:
        d = df.loc[:, ['parent', 'mut_id', 'cycle']].rename(columns={"parent":"node", "mut_id":"mutation"})
        anc = f"{treat.split('-',1)[1]} ancestor"
        r = all_data[all_data.genotype==anc].loc[:, ['genotype', 'mut_id']].rename(columns={"genotype":"node", "mut_id":"mutation"})
        r.replace(anc, treat+"_"+"ROOT", inplace=True)
        out = pd.concat([d,r]).reset_index(drop=True).fillna({'cycle':0})
        out['endpoint'] = [c == 11 or c == 0 for c in out.cycle]
        out.loc[:, ['node', 'mutation']].drop_duplicates().to_csv(os.path.join(EXPORT_TREE_PATH,treat+"_mutations.csv"))
        out.loc[out.endpoint==True, ['node', 'mutation']].drop_duplicates().to_csv(os.path.join(EXPORT_TREE_PATH,treat+"_mutations_endpoints.csv"))
In [9]:
import sqlite3
with sqlite3.connect("lce_data.sqlite") as database:
    metadata.drop(columns=['cycle','size','ancestor','treatment','tube','name','pe']).to_sql("sequenced_cultures", database, if_exists="replace", dtype={'genotype':"TEXT PRIMARY KEY"})
    mutation_description.drop(columns=['name']).to_sql("mutations", database, if_exists="replace", dtype={'mut_id':"TEXT PRIMARY KEY"})
    long_table.loc[:,["evidence","mut_id","genotype"]].to_sql("sequencing", database, if_exists="replace", index=False)
In [10]:
fig, ax = plt.subplots(1,1, figsize=(7,4))
for k, v in metadata.groupby(["ancestor", "size"]):
    ax.scatter(v.cycle-(1 if k[1] else 0)+(np.random.random(v.shape[0])-0.5)*0.8, v.n_mutations, 
               label=f"{k[1]}-{k[0]}" if k[1] else k[0]+' ancestor',
               color=COLORS[f"{k[1]}-{k[0]}"] if k[1] else 'k',
               marker='.' if k[1] else {"LCS+":"x","LCS-":"+"}[k[0]])
ax.legend()
ax.set(xlabel='Cycle', 
       xlim=(-0.5,10.5), 
       ylim=(0,420),
       ylabel="Number of mutations identified",
       xticks=np.arange(11))
for x in np.arange(11)+0.5:
    plt.axvline(x,color='k',ls='--', lw=1)
fig.savefig(os.path.join(EXPORT_PATH, 'number_of_mutations.svg'),bbox_inches='tight')
fig.savefig(os.path.join(EXPORT_PATH, 'number_of_mutations.pdf'),bbox_inches='tight')
No description has been provided for this image

Filtering out the mutations from the ancestors¶

In [11]:
filtered_out = long_table.copy()
print(f"{filtered_out.shape[0]} entries")
for anc,anc_name in (["LCS+ ancestor", "LCS+"], ["LCS- ancestor", "LCS-"]):
    anc_dt = frozenset(filtered_out[filtered_out.genotype==anc].mut_id)
    genotypes = frozenset(metadata[metadata.ancestor==anc_name].index)
    print(f"Ancestor {anc}: {len(anc_dt)} mutations to be filtered out of {len(genotypes)} datasets")
    
    idx = np.logical_not(np.logical_and([x in genotypes for x in filtered_out.genotype],
                                        [x in anc_dt for x in filtered_out.mut_id]))
    filtered_out = filtered_out.loc[idx]
    print(filtered_out.shape[0])
    
nmut2 = filtered_out.groupby("genotype").count().mut_id
metadata2 = metadata.copy()
metadata2['n_mutations'] = nmut2

fig, ax = plt.subplots(1,1, figsize=(7,4))
for k, v in metadata2.groupby(["ancestor", "size"]):
    if k[1]:
        ax.scatter(v.cycle-(1 if k[1] else 0)+(np.random.random(v.shape[0])-0.5)*0.8, v.n_mutations, 
                   label=f"{k[1]}-{k[0]}" if k[1] else k[0]+' ancestor',
                   color=COLORS[f"{k[1]}-{k[0]}"] if k[1] else 'k',
                   marker='.' if k[1] else {"LCS+":"x","LCS-":"+"}[k[0]])
ax.legend()
ax.set(xlabel='Cycle', 
       xlim=(-0.5,10.5), 
       #ylim=(-1,420),
       ylabel="Number of mutations identified",
       xticks=np.arange(11))
for x in np.arange(11)+0.5:
    plt.axvline(x,color='k',ls='--', lw=1)
plt.axhline(0,color='k',ls='--', lw=1)
fig.savefig(os.path.join(EXPORT_PATH, 'filtered_number_of_mutations.svg'),bbox_inches='tight')
fig.savefig(os.path.join(EXPORT_PATH, 'filtered_number_of_mutations.pdf'),bbox_inches='tight')
37819 entries
Ancestor LCS+ ancestor: 58 mutations to be filtered out of 108 datasets
32252
Ancestor LCS- ancestor: 74 mutations to be filtered out of 90 datasets
25965
No description has been provided for this image

A table for all sequencing data¶

In [12]:
cat = {'folder': 'Breseq', 'name': 'ID',
       'ref': 'Reference', 'coverage_ref': 'Reference', 
       'fraction_aligned_reads': 'Reads', 'n_reads': 'Reads',
       'read_files': 'Reads', 'n_read_files': 'Reads', 
       'cycle': 'ID', 'size': 'ID', 'ancestor': 'ID', 'treatment': 'ID',
       'tube': 'ID', 'genotype': 'ID', 'n_mutations':'Mutations'}
metadisp = metadata.copy().reset_index()
metadisp = pd.merge(metadisp, mut_mutations, left_on='genotype',right_index=True)
cols_mut = list(mut_mutations.columns) 
metadisp = metadisp[["name","genotype","treatment", 
        "n_mutations",
        "n_read_files", "n_reads", "fraction_aligned_reads",
        "ref","coverage_ref"]+cols_mut ]
metadisp.set_index('genotype',inplace=True)
metadisp.sort_values('n_mutations', inplace=True)
metadisp.columns = pd.MultiIndex.from_tuples([(cat[x], x)
                                              if x in cat 
                                              else (x.split()[0], x)
                                              for x in metadisp.columns
                                              ])
C0 = matplotlib.colors.to_hex('C0')
border =  [{'selector': 'th', 'props': 'border-left: 1px solid black'},
           {'selector': 'td', 'props': 'border-left: 1px solid black'}]
html = metadisp.style.set_table_styles({
    ('ID','name'): border
}, overwrite=False, axis=0)\
.set_sticky(axis="index")\
.set_sticky(axis="columns")\
.bar(subset=[('Reference','coverage_ref'), ('Reads','fraction_aligned_reads'),('Reads',"n_reads"),('Mutations',"n_mutations")], color=matplotlib.colors.to_hex('C0'))\
.map(lambda x: (f"background-color:{matplotlib.colors.to_hex(COLORS[x])}" if x in COLORS else ''), subset=[('ID','treatment')])\
.map(lambda x: (f"background-color:{C0};color:{C0};border-color:black" if x else 'color:white;border-color:black'), 
     subset=[(x.split()[0],x) for x in cols_mut])\
.map_index(lambda v: "max-width:3em;font-size:0.8em;" if 'mut' in v and len(v.split())>1 else '', axis=1)\
.to_html()

with open(os.path.join(EXPORT_PATH,'sequencing.html'),'w') as f:
    f.write("<html><head><meta charset=\"utf-8\"><title>Sequencing Data</title><style> table, th, td {border: 1px solid} table{border-collapse: collapse;} thead{background-color:white;} </style></head><body>")
    f.write(html)
    f.write("</body>")
In [13]:
key_col = lambda x: ((0,0) if ("ancestor" in x or 'ROOT' in x)
                     else (0,-1) if ("name" in x or "All" in x)
                     else ({'+':0,'-':1}[x[5]],{'S':0,'L':1}[x[0]], int(x.split(":")[1].split('-')[0])))

common_mut = mutations.set_index('mut_id').loc[list(any_ancestors),:].reset_index()
common_mut['All Sequencing'] = [m in all_sequences for m in common_mut.mut_id]
common_mut = common_mut.set_index(cols+['time_sequenced','mut_id']).drop(columns='name')
common_columns = sorted(list(common_mut.columns), key=key_col)
common_mut = common_mut.loc[:,common_columns].reset_index()

common_mut = common_mut.sort_values(['All Sequencing','LCS- ancestor', 'LCS+ ancestor','time_sequenced'], ascending=False).reset_index(drop=True)
html = common_mut.style.set_sticky(axis="index")\
.map(lambda x: (f"background-color:{C0};color:{C0};border-color:black" if x else 'color:white;border-color:black'), 
     subset=common_columns)\
.set_sticky(axis="columns")\
.to_html()

with open(os.path.join(EXPORT_PATH,'common_mutations.html'),'w') as f:
    f.write("<html><head><meta charset=\"utf-8\"><title>Sequencing Data</title><style> table, th, td {border: 1px solid} table{border-collapse: collapse;} thead{background-color:white;} </style></head><body>")
    f.write(html)
    f.write("</body>")
In [14]:
m_mut = mutations[['mut' in x for x in mutations.gene]].copy()
m_mut['mutation_name'] = [x.gene.split()[0]+' '+x.annotation.split()[0] for k,x in m_mut.iterrows()]

m_mut = m_mut.set_index(cols+['time_sequenced','mut_id',]).drop(columns=['name'])
sorted_cols = sorted(list(m_mut.columns), key=key_col)
m_mut = m_mut.loc[:,sorted_cols].reset_index().drop(columns='mut_id').set_index(cols+['mutation_name','time_sequenced'])

print("Latex Table:")
print(80*"-")
print(m_mut.reset_index().loc[:,cols].rename(columns=lambda x: x.capitalize().replace("_"," ").replace("Time", r"\# times")).set_index('Gene').to_latex().replace('_',r"\_"))
print(80*"-")
Latex Table:
--------------------------------------------------------------------------------
\begin{tabular}{lrlll}
\toprule
 & Position & Mutation & Annotation & Description \\
Gene &  &  &  &  \\
\midrule
mutY ← & 355688 & C→T & G153D (GGC→GAC) & A/G‑specific adenine glycosylase \\
mutL → & 586966 & C→T & S5L (TCG→TTG) & DNA mismatch repair protein \\
mutL → & 587050 & A→G & E33G (GAG→GGG) & DNA mismatch repair protein \\
mutL → & 587260 & T→G & L103R (CTG→CGG) & DNA mismatch repair protein \\
mutL → & 587555 & C→T & L201L (CTC→CTT) & DNA mismatch repair protein \\
mutL → & 587737 & A→G & N262S (AAC→AGC) & DNA mismatch repair protein \\
mutL → & 588126 & (G)5→6 & coding (1174/1902 nt) & DNA mismatch repair protein \\
mutL → & 588551 & C→T & A533A (GCC→GCT) & DNA mismatch repair protein \\
mutS ← & 1300356 & C→T & G707S (GGC→AGC) & DNA mismatch repair protein \\
mutS ← & 1300673 & +A & coding (1802/2592 nt) & DNA mismatch repair protein \\
mutS ← & 1300783 & A→G & R564R (CGT→CGC) & DNA mismatch repair protein \\
mutS ← & 1300986 & T→G & T497P (ACC→CCC) & DNA mismatch repair protein \\
mutS ← & 1301392 & A→G & R361R (CGT→CGC) & DNA mismatch repair protein \\
mutS ← & 1302049 & A→G & R142R (CGT→CGC) & DNA mismatch repair protein \\
mutS ← & 1302080 & T→C & N132S (AAC→AGC) & DNA mismatch repair protein \\
mutS ← / ← PFLU\_1165 & 1302576 & G→A & intergenic (‑102/+68) & DNA mismatch repair protein/putative membrane protein \\
mutM ← & 6343361 & A→G & S46P (TCC→CCC) & formamidopyrimidine‑DNA glycosylase \\
\bottomrule
\end{tabular}

--------------------------------------------------------------------------------
In [15]:
m_mut.sort_values('time_sequenced',ascending=False, inplace=True)
yticks = [f"{x[5]:>15} ({x[6]})" for x in m_mut.index]
xticks = [x for x in m_mut.columns]

fig, ax = plt.subplots(1,1, figsize=(30,5))
ax.matshow(m_mut.values!=0, cmap="bone_r")
ax.set(yticks=range(len(yticks)),
       yticklabels=yticks,
       xticks=range(len(xticks)),
       xticklabels=xticks);
ax.set_xticklabels(xticks, rotation=90);
fig.savefig(os.path.join(EXPORT_PATH,"sequenced_mut_mutations.pdf"), bbox_inches='tight')
No description has been provided for this image
In [ ]: