<jupyter_start><jupyter_text>30 April 2015Reassortment is more highly represented when there is a phylogenetic difference between the host species than when there isn't.Phylogenetic difference is defined by the 18S rRNA. What if we used cytochrome C oxidase?<jupyter_code>np.shape(bins)
def permute_18s(G):
    """
    Be sure to pass in a copy of G, and not the original!
    Takes in a graph G, shuffles the labels, and returns G with shuffled labels.
    """
    from random import shuffle
    rrna_labels = [d['18s_rRNA_sequence'] for n, d in G.nodes(data=True)]
    shuffle(rrna_labels)
    
    for i, n in enumerate(G.nodes()):
        G.node[n]['18s_rRNA_sequence'] = rrna_labels[i]
        
    return G

shuffledG = G.copy()

expected_reassort = []

from time import time

data = np.zeros(shape=(100, 10))
for i in range(100): # changed to 100 in script format.
    start_time = time()
    shuffledG = permute_18s(shuffledG)
    reassortant_shuffled_phyloD, full_complement_shuffled_phyloD = count_delta_phyloD(shuffledG)
    
    reassortant_shuffled_binned = np.histogram(reassortant_shuffled_phyloD, bins=bins)[0].astype(float)
    full_complement_shuffled_binned = np.histogram(full_complement_shuffled_phyloD, bins=bins)[0].astype(float)
    
    expected_proportion_reassortant = reassortant_shuffled_binned / (full_complement_shuffled_binned + reassortant_shuffled_binned)
    
    end_time = time()
    
    data[i] = expected_proportion_reassortant
    
    print('Round {0}, time: {1} s'.format(i, end_time - start_time))
# pd.DataFrame(data).fillna(0).to_csv('phylod_expected_reassortant_fraction.csv')

# Read in the data from running the above cells on Rous.
expected_shuffled = pd.read_csv('phylod_expected_reassortant_fraction.csv', index_col=0)
expected_shuffled = expected_shuffled[['0', '1', '2', '3', '4']]
fig = plt.figure(figsize=(3,2))
ax = fig.add_subplot(111)

bp = ax.boxplot(expected_shuffled.values, showfliers=False)
plt.setp(bp['boxes'], color='red')
plt.setp(bp['whiskers'], color='red')
plt.setp(bp['fliers'], color='red')
plt.setp(bp['caps'], color='red')
plt.setp(bp['medians'], color='red')

xpoints = np.arange(1, 11, 1)
ax.scatter(xpoints[0:5], proportion_reassortant[0:5], color='blue')
ax.plot(xpoints[0:5], proportion_reassortant[0:5], color='blue')
ax.set_ylabel('Proportion Reassortant')
ax.set_ylim(0, 1)
ax.set_xlabel('Host Dissimilarity')
ax.set_xlim(0.5,5)
ax.set_xticks(np.arange(0,6,1)+0.5)
ax.set_xticklabels(np.arange(0,0.3,0.05))

lines, labels = ax.get_legend_handles_labels()
# lines2, labels2 = ax2.get_legend_handles_labels()
ax.legend() 
ax.annotate('a.', xy=(0,1), xycoords='figure fraction', ha='left', va='top')

plt.subplots_adjust(left=0.15, bottom=0.18, right=0.95, top=0.95)
plt.savefig('Phylogenetic Susceptibility and Reassortment.pdf')
# # Plot proportion non-identical vs. identical ratio
# import numpy as np
# simulated = pd.read_csv('ratio_identical_nonidentical_expected.csv', index_col=0)
# actual = pd.read_csv('ratio_identical_nonidentical.csv', index_col=0)

# mean = simulated.mean().values
# ci99 = simulated.std() * 3
# ci99 = ci99.values

# barwidth = 0.3

# fig = plt.figure(figsize=(1, 1.5))
# ax = fig.add_subplot(111)
# ax.bar(1, mean, width=barwidth, yerr=ci99, color='blue', alpha=0.5, label='Shuffled')
# ax.bar(1-barwidth, actual.values, width=barwidth, color='blue', label='Data')
# ax.set_ylabel('Fold Representation')
# ax.set_xticks([])
# ax.set_yticks(np.arange(0,5,1))
# ax.annotate('b.', xy=(0,1), xycoords='figure fraction', ha='left', va='top')
# plt.subplots_adjust(left=0.3, bottom=0.05)
# plt.savefig('PhyloD Reassortant Fold Over-Representation.pdf')
# Plot viral PWD against host PWD
viral_pwd_reassortant = []
viral_pwd_full = []
host_pwd_reassortant = []
host_pwd_full = []
# edge_reassortant_status = []
for sc, sk, d in G.edges(data=True):
    sc_18s = G.node[sc]['18s_rRNA_sequence'].seq.ungap('-')
    sk_18s = G.node[sk]['18s_rRNA_sequence'].seq.ungap('-')

    if len(sc_18s) != 0 and len(sk_18s) != 0:
        dist = host_phylogenetic_distance(G, sc, sk)
        
    #viral_pwd.append(d['pwi'])
    #host_pwd.append(dist)
    if d['edge_type'] == 'full_complement':
        viral_pwd_full.append(d['pwi'])
        host_pwd_full.append(dist)
    if d['edge_type'] == 'reassortant':
        viral_pwd_reassortant.append(d['pwi'])
        host_pwd_reassortant.append(dist)
print(8 - np.mean(viral_pwd_reassortant))
print(8 - np.mean(viral_pwd_full))
print(np.mean(host_pwd_reassortant))
print(np.mean(host_pwd_full))
(sum(Counter(host_pwd_full).values()) - Counter(host_pwd_full)[0]) / float(Counter(host_pwd_full)[0])
(sum(Counter(host_pwd_reassortant).values()) - Counter(host_pwd_reassortant)[0]) / float(Counter(host_pwd_reassortant)[0])
# plt.hist2d
from matplotlib.colors import LogNorm 

fig = plt.figure()#figsize=(3,2))
ax = fig.add_subplot(111)
ax.scatter(host_pwd_reassortant, 8 - np.array(viral_pwd_reassortant), color='blue', label='Reassortant', alpha=0.1)
ax.scatter(host_pwd_full, 8 - np.array(viral_pwd_full), color='green', label='Full Complement', alpha=0.1)
ax.legend()
from scipy.stats import linregress

slope_r, intercept_r, r_value_r, p_value_r, std_error_r = linregress(8 - np.array(viral_pwd_reassortant), host_pwd_reassortant)
slope_f, intercept_f, r_value_f, p_value_f, std_error_f = linregress(8 - np.array(viral_pwd_full), host_pwd_full)
print(slope_r, slope_f)
print(r_value_r**2, r_value_f**2)
print(p_value_r, p_value_f)
print(std_error_r, std_error_f)
# plt.scatter([0] * len(host_pwd_full), host_pwd_full, alpha=0.1, color='green')
# plt.scatter([1] * len(host_pwd_reassortant), host_pwd_reassortant, alpha=0.1, color='blue')
plt.hist(host_pwd_full)
plt.hist(host_pwd_reassortant)
# Add in ecotype information

# Load the data
df = pd.read_csv('Country-Species_to_be_Coded_by_Nichola.csv', index_col=0)
df['Habitat/setting'] = df['Habitat/setting'].replace('Domestic ', 'Domestic')
# Make it into a dicitonary where each entry is (country, species):habitat/setting
ecotype_map = dict()

for row, data in df.iterrows():
    country = data['Country']
    species = data['Species']
    ecotype = data['Habitat/setting']

    ecotype_map[(country, species)] = ecotype

# Add the ecotype into each node.
for n, d in G.nodes(data=True):
    country = d['country']
    species = d['host_species']
    if '/' in species:
        species = species.split('/')[0]
    if 'null' in species:
        species = 'Unknown'


    
    G.node[n]['ecotype'] = ecotype_map[(country, species)]
def correct_metadata(G, node):
    """
    Helper function for providing the correct subtype graph metadata.    
    """
    d = G.node[node]
    subtype = d['subtype']
    ecotype = d['ecotype']
    host = d['host_species']
    if ecotype == 'Unknown':
        ecotype = 'W'
    if host == 'Human':
        ecotype = 'H'
    if ecotype == 'Domestic':
        ecotype = 'D'
    if ecotype == 'Wild':
        ecotype = 'W'
    subtype = subtype + '::' + ecotype
    
    return ecotype, subtype, host
    
from collections import Counter, defaultdict
transition_hostpd_reassortant = defaultdict(list)
transition_hostpd_full = defaultdict(list)
# Plot the number of ecotype jumps detected.
for sc, sk, d in G.edges(data=True):
    sc_ecotype, sc_subtype, sc_host = correct_metadata(G, sc)
    sk_ecotype, sk_subtype, sk_host = correct_metadata(G, sk)

    sc_18s = G.node[sc]['18s_rRNA_sequence'].seq.ungap('-')
    sk_18s = G.node[sk]['18s_rRNA_sequence'].seq.ungap('-')


    transition = sc_ecotype + '::' + sk_ecotype
    if len(sc_18s) != 0 and len(sk_18s) != 0:
        dist = host_phylogenetic_distance(G, sc, sk)
        if d['edge_type'] == 'reassortant':
            transition_hostpd_reassortant[transition].append(dist)
        else:
            transition_hostpd_full[transition].append(dist)
transition_hostpd_means = defaultdict(list)
for k, v in transition_hostpd.items():
    transition_hostpd_means[k] = np.mean(v)

fig = plt.figure()
ax = fig.add_subplot(111)
# ax.bar(np.arange(len(transition_hostpd_means.keys()))-0.5, transition_hostpd_means.values())
# ax.set_xticks(np.arange(len(transition_hostpd_means.keys())))
# ax.set_xticklabels(transition_hostpd_means.keys())
dataset = sorted(transition_hostpd_reassortant.items(), key=lambda x:x[0])
compiled = []
for i, (transition, hostpds) in enumerate(dataset):
    #np.array(hostpds)
    compiled.append(np.array(hostpds))
    print(len(np.array(hostpds)))
    ax.scatter([i] * len(hostpds), hostpds, alpha=0.1, color='blue')
    
dataset = sorted(transition_hostpd_full.items(), key=lambda x:x[0])
compiled = []
for i, (transition, hostpds) in enumerate(dataset):
    #np.array(hostpds)
    compiled.append(np.array(hostpds))
    print(len(np.array(hostpds)))
    ax.scatter([i] * len(hostpds), hostpds, alpha=0.1, color='green')
    


# Try tree patristic distance rather than hamming distance.
import dendropy as dp

tree = dp.Tree().get_from_path('RAxML_bestTree.host_COI5P_trimmed_alignment.nexus', 'nexus')

pds = dp.treecalc.PatristicDistanceMatrix(tree).distances()
taxon_dict = {str(i).split(': ')[-1]:i for i in tree.taxon_set}
taxon_dict
network_pds = dict()
network_pds['reassortant'] = list()
network_pds['full_complement'] = list()

for sc, sk, d in G.edges(data=True):
    sc_host = G.node[sc]['host_species']
    sk_host = G.node[sk]['host_species']
    if sc_host in taxon_dict.keys() and sk_host in taxon_dict.keys():
        patr_dist = dp.treecalc.patristic_distance(tree, taxon_dict[sc_host], taxon_dict[sk_host])
        network_pds[d['edge_type']].append(patr_dist)
plt.hist(network_pds['reassortant'])
plt.hist(network_pds['full_complement'])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14859147662779898, 'fraction_numerical': 0.019812196883706532, 'mean_word_length': 5.0879396984924625, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '105294', 'license': 'MIT', 'n_tokens_mistral': 3859, 'n_tokens_neox': 3693, 'n_words': 669, 'path': '16 Does phylogenetic distance of host impact reassortment?.ipynb', 'repo_name': 'ericmjl/influenza-reassortment-analysis', 'size': '9691'}