Chapter IV - Reproducing Kiefl et al, 2022

#! /usr/bin/env python

import pandas as pd
from pathlib import Path

tables_dir = Path('WW_TABLES')
tables_dir.mkdir(exist_ok=True)

sample_metadata = pd.read_csv("00_SAMPLE_INFO_FULL.txt", sep='\t')
sample_metadata = sample_metadata[[col for col in sample_metadata if "Used_" not in col]]
sample_metadata = sample_metadata[~sample_metadata["sample_id"].isnull()]
ftp_links = pd.read_csv("00_FTP_LINKS", sep='\t', names=["link"])
cov = pd.read_csv("07_SUMMARY_ALL/bins_across_samples/mean_coverage.txt", sep='\t')
q2q3 = pd.read_csv("07_SUMMARY_ALL/bins_across_samples/mean_coverage_Q2Q3.txt", sep='\t')
det = pd.read_csv("07_SUMMARY_ALL/bins_across_samples/detection.txt", sep='\t')
per = pd.read_csv("07_SUMMARY_ALL/bins_across_samples/bins_percent_recruitment.txt", sep='\t')
himb083_genes = pd.read_csv("07_SUMMARY_ALL/bin_by_bin/HIMB083/HIMB083-gene_coverages.txt", sep='\t')

with pd.ExcelWriter(tables_dir/'RR.xlsx') as writer:
    sample_metadata.to_excel(writer, sheet_name='Sample identifiers')
    ftp_links.to_excel(writer, sheet_name='Sample FTP links')
    cov.to_excel(writer, sheet_name='Coverage')
    q2q3.to_excel(writer, sheet_name='Coverage Q2Q3')
    det.to_excel(writer, sheet_name='Detection')
    per.to_excel(writer, sheet_name='% recruitment')
    himb083_genes.to_excel(writer, sheet_name='HIMB083 gene coverages')

#! /usr/bin/env python

import numpy as np
import pandas as pd
import argparse
import anvio.dbops as dbops
import anvio.bamops as bamops

ap = argparse.ArgumentParser()

ap.add_argument('--contigs-db', help="Filepath to the contigs db", required=True)
ap.add_argument('--goi', help="Genes of interest (file of gene caller ids)", required=True)
ap.add_argument('--bams', help="File of bampaths. First column is sample name, second column is bam path. Header should be sample_id\tpath", required=True)
ap.add_argument('--output', help="Output text file", required=True)

args = ap.parse_args()

# ---------------------------------------------------

bams = pd.read_csv(args.bams, sep='\t')
bams = dict(zip(bams['sample_id'], bams['path']))

cdb = dbops.ContigsDatabase(args.contigs_db)
gene_calls = cdb.db.get_table_as_dataframe('genes_in_contigs').set_index('gene_callers_id')
gene_calls = gene_calls[gene_calls.index.isin([int(x.strip()) for x in open(args.goi).readlines()])]

d = {sample_id: [] for sample_id in bams.keys()}
i = 0
for sample_id, bam_path in bams.items():
    bam = bamops.BAMFileObject(bam_path)

    j = 0
    for gene_id, row in gene_calls.iterrows():
        print(f"Sample {sample_id} ({i}/{len(bams)}); Gene {gene_id} ({j}/{gene_calls.shape[0]})")

        percent_ids = []
        for read in bam.fetch_and_trim(row['contig'], row['start'], row['stop']):
            read.vectorize()
            percent_ids.append(read.get_percent_id())

        percent_ids = np.array(percent_ids)
        d[sample_id].append(np.mean(percent_ids))

        j += 1
    i += 1

    bam.close()

d['gene_callers_id'] = gene_calls.index.tolist()
pd.DataFrame(d).to_csv(args.output, sep='\t', index=False)

#! /usr/bin/env python

import pandas as pd
from pathlib import Path

bam_dir = Path('04_MAPPING/SAR11_clade')
soi = [sample.strip() for sample in open('soi').readlines()]

bam_paths = dict(
    sample_id = [],
    path = [],
)

for bam_path in bam_dir.glob("*.bam"):
    sample_id = bam_path.stem
    if sample_id not in soi:
        continue
    bam_paths['sample_id'].append(sample_id)
    bam_paths['path'].append(str(bam_path))

pd.DataFrame(bam_paths).to_csv("mgx_bam_paths", sep='\t', index=False)

#! /usr/bin/env python

import pandas as pd

goi = [int(x.strip()) for x in open('goi').readlines()]
df = pd.read_csv("07_SUMMARY_PAN/SAR11_gene_clusters_summary.txt", sep='\t')

gc = df.loc[(df['genome_name'] == 'HIMB083') & (df['gene_callers_id'].isin(goi)), 'gene_cluster_id']
gc = gc.value_counts()[gc.value_counts() == 1].index.tolist()

with open('gcoi', 'w') as f:
    f.write("\n".join([str(x) for x in gc]))

#! /usr/bin/env python

import pandas as pd

from prody.sequence.msafile import parseMSA
from pathlib import Path

genome_name = 'HIMB083'
msa_paths = Path('18_HIMB83_CORE_COMPARED_TO_PANGENOME/').glob('*.fa')

dd = {
    'gene_callers_id': [],
    'percent_id': [],
}

for path in msa_paths:
    msa = parseMSA(str(path), format='FASTA')

    for i, label in enumerate(msa._labels):
        if genome_name in label:
            ref_id = i
    is_genome_of_interest = [genome_name in defline for defline in msa._labels]
    reference = msa._msa[is_genome_of_interest,:][0]
    ref_label = msa._labels[ref_id]
    ref_gene_id = int(ref_label.split('|')[-1].split(':')[-1])

    labels = [defline for defline in msa._labels if genome_name not in defline]
    sequences = msa._msa[[not x for x in is_genome_of_interest],:]
    sequences = dict(zip(labels, sequences))

    d = {
        'gene_callers_id': [],
        'match': [],
        'mismatch': [],
        'total': [],
    }

    total_nucleotides = 0
    for label, sequence in sequences.items():
        match, mismatch = 0, 0

        for ref, seq in zip(reference, sequence):
            ref, seq = ref.decode('utf-8'), seq.decode('utf-8')

            if ref == '-' or seq == '-':
                # only compare fully aligned nucleotides
                continue

            if ref == seq:
                match += 1
            else:
                mismatch += 1

        d['gene_callers_id'].append(ref_gene_id)
        d['match'].append(match)
        d['mismatch'].append(mismatch)
        d['total'].append(match + mismatch)

        total_nucleotides += match + mismatch

    df = pd.DataFrame(d)
    df['percent_id'] = df['match']/df['total']
    df['weight'] = df['total']/df['total'].sum()

    avg_percent_id = (df['percent_id']*df['weight']).sum() * 100

    dd['gene_callers_id'].append(ref_gene_id)
    dd['percent_id'].append(avg_percent_id)

df = pd.DataFrame(dd)
df.to_csv("18_PERCENT_ID_PANGENOME.txt", sep="\t", index=False)

#! /usr/bin/env python

from pathlib import Path
from prody.database.pfam import fetchPfamMSA
from prody.sequence.msafile import parseMSA

import numpy as np
import pandas as pd
import argparse
import subprocess
import anvio.utils as utils
import anvio.dbops as dbops

ap = argparse.ArgumentParser()

ap.add_argument('--contigs-db', help="Filepath to the contigs db", required=True)
ap.add_argument('--genome-name', help="e.g. HIMB083", required=True)

args = ap.parse_args()

# --------------------------------------------------------------------------------------

cdb = dbops.ContigsDatabase(args.contigs_db)
functions = cdb.db.get_table_as_dataframe('gene_functions')
cdb.disconnect()

def get_best_pfam(gene_id):
    try:
        return functions.\
            query(f'gene_callers_id == {gene_id} & source == "Pfam"').\
            sort_values(by='e_value').\
            iloc[0, :]\
            ['accession']
    except IndexError:
        return None

# --------------------------------------------------------------------------------------

dd = {
    'gene_callers_id': [],
    'percent_id': [],
}

goi = [int(x.strip()) for x in open("goi").readlines()]
contigs_db = dbops.ContigsSuperclass(args)
for i, gene_id in enumerate(goi):
    print(f"Gene {gene_id} ({i}/{len(goi)})")
    pfam = get_best_pfam(gene_id)
    if not pfam:
        continue

    Path("18_HIMB83_CORE_COMPARED_TO_PFAM").mkdir(exist_ok=True)
    gene_fasta = Path("18_HIMB83_CORE_COMPARED_TO_PFAM") / f"{gene_id}.fa"
    pfam_msa_path = Path('18_HIMB83_CORE_COMPARED_TO_PFAM') / f"{pfam}_msa_for_{gene_id}.fa"
    default_dump_path = Path(f'{pfam}_seed.fasta')

    # Get the Pfam MSA
    try:
        pfam_versionless = ''.join(pfam.split('.')[:-1])
        url = f"https://pfam.xfam.org/family/{pfam_versionless}/alignment/seed/format?format=fasta&alnType=seed&order=t&case=u&gaps=dashes&download=1"
        utils.download_file(url=url, output_file_path=str(default_dump_path), check_certificate=False)
    except:
        print("FAILED")
        continue

    default_dump_path.replace(pfam_msa_path)

    # Get the gene AA sequence
    contigs_db.get_sequences_for_gene_callers_ids( gene_caller_ids_list=[gene_id], output_file_path=str(gene_fasta), report_aa_sequences=True)

    # massage gene into the MSA
    cmdline = f'muscle -profile -in1 {gene_fasta} -in2 {pfam_msa_path} -out {pfam_msa_path}'
    muscle_stdout = subprocess.call(cmdline, shell=True)

    # load the final MSA
    msa = parseMSA(str(pfam_msa_path), format='FASTA')

    for i, label in enumerate(msa._labels):
        if args.genome_name in label:
            ref_id = i
    is_genome_of_interest = [args.genome_name in defline for defline in msa._labels]
    reference = msa._msa[is_genome_of_interest,:][0]

    labels = [defline for defline in msa._labels if args.genome_name not in defline]
    sequences = msa._msa[[not x for x in is_genome_of_interest],:]
    sequences = dict(zip(labels, sequences))

    d = {
        'gene_callers_id': [],
        'match': [],
        'mismatch': [],
        'total': [],
    }

    total_nucleotides = 0
    for label, sequence in sequences.items():
        match, mismatch = 0, 0

        for ref, seq in zip(reference, sequence):
            ref, seq = ref.decode('utf-8'), seq.decode('utf-8')

            if ref == '-' or seq == '-':
                # only compare fully aligned nucleotides
                continue

            if ref == seq:
                match += 1
            else:
                mismatch += 1

        d['gene_callers_id'].append(gene_id)
        d['match'].append(match)
        d['mismatch'].append(mismatch)
        d['total'].append(match + mismatch)

        total_nucleotides += match + mismatch

    df = pd.DataFrame(d)
    df['percent_id'] = df['match']/df['total']
    df['weight'] = df['total']/df['total'].sum()

    avg_percent_id = (df['percent_id']*df['weight']).sum() * 100

    dd['gene_callers_id'].append(gene_id)
    dd['percent_id'].append(avg_percent_id)

df = pd.DataFrame(dd)
df.to_csv("18_PERCENT_ID_PFAM.txt", sep="\t", index=False)

#! /usr/bin/env python

import pandas as pd
from pathlib import Path

tables_dir = Path('WW_TABLES')
tables_dir.mkdir(exist_ok=True)

# ---------------------------------------------------------

df = pd.read_csv("11_SCVs.txt", sep='\t')

per_gene_per_sample = df.\
    groupby(['corresponding_gene_call', 'sample_id'])\
    [['pN_popular_consensus', 'pS_popular_consensus']].\
    mean().\
    reset_index().\
    rename(columns={'pN_popular_consensus':'mean(pN_site)', 'pS_popular_consensus':'mean(pS_site)'})

per_gene = df.\
    groupby(['corresponding_gene_call'])\
    [['pN_popular_consensus', 'pS_popular_consensus']].\
    mean().\
    reset_index().\
    rename(columns={'pN_popular_consensus':'mean(pN_site)', 'pS_popular_consensus':'mean(pS_site)'})

per_sample = df.\
    groupby(['sample_id'])\
    [['pN_popular_consensus', 'pS_popular_consensus']].\
    mean().\
    reset_index().\
    rename(columns={'pN_popular_consensus':'mean(pN_site)', 'pS_popular_consensus':'mean(pS_site)'})

overall = df\
    [['pN_popular_consensus', 'pS_popular_consensus']].\
    mean().\
    reset_index().\
    rename(columns={'index':'rate', 0:'value'})
overall['rate'] = overall['rate'].map({'pN_popular_consensus': 'mean(pN_site)', 'pS_popular_consensus': 'mean(pS_site)'})

with pd.ExcelWriter(tables_dir/'PNPS_SUMS.xlsx') as writer:
    overall.to_excel(writer, sheet_name='Overall average per site rates')
    per_gene_per_sample.to_excel(writer, sheet_name='Averaged each gene-sample pair')
    per_gene.to_excel(writer, sheet_name='Averaged for each gene')
    per_sample.to_excel(writer, sheet_name='Averaged for each sample')

#! /usr/bin/env python

import argparse
import anvio.structureops as sops
import anvio.filesnpaths as filesnpaths
import anvio.utils as utils

from prody.proteins.pdbfile import parsePDB
from prody.proteins.compare import matchChains
from prody.measure.transform import calcRMSD
from prody.measure.transform import calcTransformation

ap = argparse.ArgumentParser()
ap.add_argument('-s', '--structure')
args = ap.parse_args()

sdb = sops.StructureDatabase(args.structure)
templates = sdb.db.get_table_as_dataframe('templates')
models = sdb.db.get_table_as_dataframe('models')
goi = set([int(x.strip()) for x in open('goi').readlines()])
templates = templates[templates['corresponding_gene_call'].isin(goi)]

# downsize templates to only include the best template, where best template is defined as
# the template with the highest (alignment coverage * percent similarity = proper percent similarity)
templates = templates.groupby('corresponding_gene_call').apply(lambda df: df[df['proper_percent_similarity'] == df['proper_percent_similarity'].max()])
templates = templates.drop_duplicates(subset=['corresponding_gene_call', 'proper_percent_similarity'])

def get_RMSD(template_path, gene_path):
    x = parsePDB(template_path)
    y = parsePDB(gene_path)

    try:
        matches = matchChains(x, y, seqid=30, overlap=30)
    except:
        matches = None

    if matches is None:
        return matches

    calcTransformation(matches[0][0], matches[0][1]).apply(matches[0][0])
    return calcRMSD(matches[0][0], matches[0][1])

for i, row in templates.iterrows():
    template = row['pdb_id'] + row['chain_id']
    gene_id = row['corresponding_gene_call']

    # export structures
    try:
         template_path = utils.download_protein_structure(protein_code=row['pdb_id'], chain=row['chain_id'], output_path=filesnpaths.get_temp_file_path(), raise_if_fail=True)
    except:
        continue
    gene_path = sdb.export_pdb_content(gene_id, filesnpaths.get_temp_file_path(), ok_if_exists=True)

    templates.loc[i, 'RMSD'] = get_RMSD(template_path, gene_path)

models = models[models['picked_as_best'] == 1]
templates.reset_index(drop=True, inplace=True)
df = templates.merge(models, on='corresponding_gene_call', how='left')
df = df[(df.RMSD <= 7.5) & (df.GA341_score >= 0.95)]

with open('12_GENES_WITH_GOOD_STRUCTURES_MODELLER', 'w') as f:
    for gene_id in df['corresponding_gene_call']:
        f.write(f"{gene_id}\n")

#! /usr/bin/env python

import anvio.utils as utils
import anvio.terminal as terminal
import anvio.structureops as sops
import anvio.filesnpaths as filesnpaths

from prody import confProDy
from prody.measure.measure import calcGyradius, calcDistance
from prody.proteins.pdbfile import parsePDB
from prody.proteins.compare import matchChains
from prody.measure.transform import calcRMSD, calcTransformation

import numpy as np
import pandas as pd
import tmscoring
import matplotlib.pyplot as plt

from pathlib import Path

# -------------------------------------------------------
# Load up the structure databases and auxiliary info
# -------------------------------------------------------

confProDy(verbosity='none')
progress = terminal.Progress()
run = terminal.Run()

db_AF = sops.StructureDatabase('09_STRUCTURE.db')
db_MOD = sops.StructureDatabase('09_STRUCTURE_MOD.db')

plddt = pd.read_csv('09_STRUCTURES_AF/pLDDT_gene.txt', sep='\t').set_index('gene_callers_id')
plddt_residue = pd.read_csv('09_STRUCTURES_AF/pLDDT_residue.txt', sep='\t')

# -------------------------------------------------------
# Get intersection considered trustworthy for both methods
# -------------------------------------------------------

good_AF = set([int(x.strip()) for x in open('12_GENES_WITH_GOOD_STRUCTURES').readlines()])
good_MOD = set([int(x.strip()) for x in open('12_GENES_WITH_GOOD_STRUCTURES_MODELLER').readlines()])
good_intersect = good_AF.intersection(good_MOD)

run.info('Number of trustworthy AlphaFold predictions', len(good_AF))
run.info('Number of trustworthy MODELLER predictions', len(good_MOD))
run.info('Number of shared trustworthy predictions', len(good_intersect))

# -------------------------------------------------------
# Define metrics for assessing ALIGNMENTS
# -------------------------------------------------------

def get_RMSD(path1, path2):
    struct1 = parsePDB(path1)
    struct2 = parsePDB(path2)
    matches = matchChains(struct1, struct2, seqid=30, overlap=30)
    calcTransformation(matches[0][0], matches[0][1]).apply(matches[0][0])
    return calcRMSD(matches[0][0], matches[0][1])

def get_TM_score(path1, path2):
    """Get TM score using tmscoring.

    For whatever reason, there is a bit of an issue with tmscoring. In a rare number of cases, the
    TM score returned depends upon the order of filepaths passed to tmscoring.TMscoring, and by
    trial and error, I have found that taking the max of these values consistently yields the same
    results as the original TM web server https://zhanggroup.org/TM-score/. As such, I return the
    maximum of both values here. See issue: https://github.com/Dapid/tmscoring/issues/6
    """
    aln = tmscoring.TMscoring(path1, path2)
    aln.optimise()
    TM_score = aln.tmscore(**aln.get_current_values())

    aln = tmscoring.TMscoring(path2, path1)
    aln.optimise()
    TM_score_2 = aln.tmscore(**aln.get_current_values())

    return max(TM_score, TM_score_2)

def get_contact_map_mean_absolute_error(gene_id):
    structure_AF = db_AF.get_structure(gene_id)
    contact_map_AF = get_contact_map(structure_AF)
    structure_MOD = db_MOD.get_structure(gene_id)
    contact_map_MOD = get_contact_map(structure_MOD)
    MAE = np.abs(contact_map_MOD - contact_map_AF).sum() / np.product(contact_map_MOD.shape)
    return MAE

def get_contact_map(structure):
    d = {}
    protein_length = len(structure.structure)
    contact_map = np.zeros((protein_length, protein_length))
    for i, residue1 in enumerate(structure.structure):
        if i not in d:
            d[i] = structure.get_residue_center_of_mass(residue1)
        COM1 = d[i]
        for j, residue2 in enumerate(structure.structure):
            if i > j:
                contact_map[i, j] = contact_map[j, i]
            else:
                if j not in d:
                    d[j] = structure.get_residue_center_of_mass(residue2)
                COM2 = d[j]
                contact_map[i, j] = np.sqrt(np.sum((COM1 - COM2)**2))
    return contact_map

# -------------------------------------------------------
# Define metrics for assessing STRUCTURES
# -------------------------------------------------------

def get_Rg(path):
    """Get the radius of gyration"""
    struct = parsePDB(path).select('not hydrogen')
    return calcGyradius(struct)

def get_end_to_end(path):
    """Get the radius of gyration"""
    struct = parsePDB(path).select('not hydrogen')
    nter, cter = struct[0], struct[-1]
    return calcDistance(nter, cter)

def get_mean_RSA(gene_id, db):
    return db.get_residue_info_for_gene(gene_id)['rel_solvent_acc'].mean()

def get_frac_sec_struct(gene_id, db):
    res_info = db.get_residue_info_for_gene(gene_id)
    alpha = res_info['sec_struct'].isin(['G','H','I']).sum()/len(res_info)
    beta = res_info['sec_struct'].isin(['E','B']).sum()/len(res_info)
    loop = res_info['sec_struct'].isin(['S','T','C']).sum()/len(res_info)
    return alpha, beta, loop

def get_length(gene_id, db1, db2):
    try:
        return len(db1.get_structure(gene_id).get_sequence())
    except:
        return len(db2.get_structure(gene_id).get_sequence())

# -------------------------------------------------------
# Run the metrics for each gene
# -------------------------------------------------------

d = dict(
    gene_callers_id = [],
    length = [],
    has_AF = [],
    has_MOD = [],
    RMSD = [],
    TM_score = [],
    contact_MAE = [],
    mean_RSA_AF = [],
    mean_RSA_MOD = [],
    gyradius_AF = [],
    gyradius_MOD = [],
    endtoend_AF = [],
    endtoend_MOD = [],
    frac_alpha_AF = [],
    frac_alpha_MOD = [],
    frac_beta_AF = [],
    frac_beta_MOD = [],
    frac_loop_AF = [],
    frac_loop_MOD = [],
    template_similarity = [],
    mean_pLDDT = [],
    mean_trunc_pLDDT = [],
)

progress.new('Calculating metrics', progress_total_items=len(good_AF.union(good_MOD)))

for gene_id in good_AF.union(good_MOD):
    has_AF = True if gene_id in good_AF else False
    has_MOD = True if gene_id in good_MOD else False

    path_AF = db_AF.export_pdb_content(gene_id, filesnpaths.get_temp_file_path()) if has_AF else np.nan
    # Calculate structure metrics
    mean_RSA_AF = get_mean_RSA(gene_id, db_AF) if has_AF else np.nan
    gyradius_AF = get_Rg(path_AF) if has_AF else np.nan
    endtoend_AF = get_end_to_end(path_AF) if has_AF else np.nan
    frac_alpha_AF, frac_beta_AF, frac_loop_AF = get_frac_sec_struct(gene_id, db_AF) if has_AF else (np.nan, np.nan, np.nan)
    # AF-specific metrics
    mean_pLDDT = plddt.loc[gene_id, 'plddt'] if has_AF else np.nan
    mean_trunc_pLDDT = plddt_residue.loc[plddt_residue['gene_callers_id']==gene_id, 'plddt'][10:-10].mean() if has_AF else np.nan

    path_MOD = db_MOD.export_pdb_content(gene_id, filesnpaths.get_temp_file_path()) if has_MOD else np.nan
    # Calculate structure metrics
    mean_RSA_MOD = get_mean_RSA(gene_id, db_MOD) if has_MOD else np.nan
    gyradius_MOD = get_Rg(path_MOD) if has_MOD else np.nan
    endtoend_MOD = get_end_to_end(path_MOD) if has_MOD else np.nan
    frac_alpha_MOD, frac_beta_MOD, frac_loop_MOD = get_frac_sec_struct(gene_id, db_MOD) if has_MOD else (np.nan, np.nan, np.nan)
    # Calculate MOD-specific metrics
    template_similarity = db_MOD.get_template_info_for_gene(gene_id)['percent_similarity'].mean() if has_MOD else np.nan

    # Calculate alignment metrics
    RMSD = get_RMSD(path_AF, path_MOD) if (has_MOD and has_AF) else np.nan
    TM_score = get_TM_score(path_AF, path_MOD) if (has_MOD and has_AF) else np.nan
    MAE = get_contact_map_mean_absolute_error(gene_id) if (has_MOD and has_AF) else np.nan

    progress.update(f"Gene ID {gene_id}, RMSD {RMSD:.1f}, contact MAE {MAE:.1f}, TM score {TM_score:.2f}")
    progress.increment()

    d['gene_callers_id'].append(gene_id)
    d['length'].append(get_length(gene_id, db_AF, db_MOD))
    d['has_AF'].append(has_AF)
    d['has_MOD'].append(has_MOD)
    d['RMSD'].append(RMSD)
    d['TM_score'].append(TM_score)
    d['contact_MAE'].append(MAE)
    d['mean_RSA_AF'].append(mean_RSA_AF)
    d['mean_RSA_MOD'].append(mean_RSA_MOD)
    d['gyradius_AF'].append(gyradius_AF)
    d['gyradius_MOD'].append(gyradius_MOD)
    d['endtoend_AF'].append(gyradius_AF)
    d['endtoend_MOD'].append(gyradius_MOD)
    d['frac_alpha_AF'].append(frac_alpha_AF)
    d['frac_alpha_MOD'].append(frac_alpha_MOD)
    d['frac_beta_AF'].append(frac_beta_AF)
    d['frac_beta_MOD'].append(frac_beta_MOD)
    d['frac_loop_AF'].append(frac_loop_AF)
    d['frac_loop_MOD'].append(frac_loop_MOD)
    d['template_similarity'].append(template_similarity)
    d['mean_pLDDT'].append(mean_pLDDT)
    d['mean_trunc_pLDDT'].append(mean_trunc_pLDDT)

progress.end()

df = pd.DataFrame(d)
df.to_csv('09_STRUCTURE_comparison.txt', sep='\t', index=False)

tables_dir = Path('WW_TABLES')
tables_dir.mkdir(exist_ok=True)
with pd.ExcelWriter(tables_dir/'STRUCT_COMP.xlsx') as writer:
    df.to_excel(writer, sheet_name='Structure comparison')

#! /usr/bin/env python

import anvio
import pandas as pd
import anvio.structureops as sops

from pathlib import Path

tables_dir = Path('WW_TABLES')
tables_dir.mkdir(exist_ok=True)

# ------------------------------------------------------------
# Load up all of the dataframes
# ------------------------------------------------------------

sheet1 = pd.read_csv("08_INTERACDOME-binding_frequencies.txt", sep='\t')
sheet2 = pd.read_csv("08_INTERACDOME-domain_hits.txt", sep='\t')
sheet3 = pd.read_csv("08_INTERACDOME-match_state_contributors.txt", sep='\t')
sheet4 = pd.read_csv(Path(anvio.__file__).parent / 'data/misc/Interacdome/representable_interactions.txt', sep='\t', comment='#')

# ------------------------------------------------------------
# Save the excel sheets
# ------------------------------------------------------------

with pd.ExcelWriter(tables_dir/'LIG.xlsx') as writer:
    sheet1.to_excel(writer, sheet_name='Ligand binding frequencies')
    sheet2.to_excel(writer, sheet_name='Domain hits')
    sheet3.to_excel(writer, sheet_name='Match state contributions')
    sheet4.to_excel(writer, sheet_name='Representable interactions')