<jupyter_start><jupyter_text>Statistics of the sequencing data:- counts: samples, species- amplicon recovery- allele counts per sample- extent of gaps in the alignmentManual steps: none<jupyter_code>%run common.ipynb<jupyter_output><empty_output><jupyter_text>Read files<jupyter_code>seq = pd.read_csv(COMB_SEQ)
display(seq.shape)
seq.sample()
meta = pd.read_csv(COMB_META)
display(meta.shape)
meta.sample()<jupyter_output><empty_output><jupyter_text>General stats<jupyter_code># samples per source (run/reference)
meta.groupby('Source')['s_Sample'].nunique()
# samples per sample set
meta.groupby('Sample Set')['s_Sample'].nunique()
# species total
meta.Species.nunique()
# species per sample set
meta.groupby('Sample Set')['Species'].nunique()
# samples per species
meta.groupby('Species')['s_Sample'].nunique().describe()<jupyter_output><empty_output><jupyter_text>Overall missingness<jupyter_code># for some overall statistics, exclude reference-derived data
miseq = seq[~seq.s_Sample.str.startswith('anopheles')]
# amplicon recovery per sample
aps = miseq.groupby('s_Sample')['target'].nunique()
display(aps.describe())
plt.hist(aps, bins=range(62))
plt.axvline(35, color='r', linestyle='solid')
plt.xlabel('Amplicons recovered');
# same, combined with reference data
aps = seq.groupby('s_Sample')['target'].nunique()
display(aps.describe())
plt.hist(aps, bins=range(62))
plt.axvline(35, color='r', linestyle='solid')
plt.xlabel('Amplicons recovered');<jupyter_output><empty_output><jupyter_text>Missingness vs species<jupyter_code># assign taxonomy
def get_taxonomy(species_series, taxonomy):
    tax_df = pd.DataFrame(species_series)
    ampl_counts['subgenus'] =  ampl_counts.species.replace(taxonomy['subgenus'].to_dict())
    ampl_counts['series'] =  ampl_counts.species.replace(taxonomy['series'].to_dict())
# recovered amplicon counts 
ampl_counts = seq.groupby(['species','s_Sample']) \
                    ['target'].nunique() \
                    .reset_index() \
                    .sort_values(ascending=False, 
                                 by=['target','species']) \
                    .rename(columns={'target':'targets_recovered'})
                    
taxonomy = pd.read_csv(TAXONOMY_FILE, index_col=0)
# set subgenera species mapping
ampl_counts['subgenus'] =  ampl_counts.species.replace(taxonomy['subgenus'].to_dict())
ampl_counts['series'] =  ampl_counts.species.replace(taxonomy['series'].to_dict())

ampl_counts.subgenus.value_counts(dropna=False)
ampl_counts.groupby('subgenus')['species'].nunique()

# series resolution for cellia
ampl_counts['group'] = ampl_counts['subgenus'] + ' subgenus'
ampl_counts.loc[ampl_counts.subgenus=='Cellia', 'group'] = ampl_counts.loc[ampl_counts.subgenus=='Cellia', 'series'] + ' series'
ampl_counts.group.value_counts(dropna=False)
ampl_counts[ampl_counts.group.isna()].species.unique()
# shorten genus name in species
ampl_counts['species'] = ampl_counts['species'].str.replace('Anopheles_','An. ')
# numeric representation of species
all_spp = ampl_counts.species.drop_duplicates()
spp_no = dict()
for i, sp in enumerate(all_spp):
    spp_no[sp] = i
ampl_counts['species_no'] = ampl_counts['species'].replace(spp_no)
# jitter in case of equal targets recovered
j = .4
for spn in ampl_counts.species_no.unique():
    acs = ampl_counts[ampl_counts.species_no == spn]
    for tr in acs.targets_recovered.unique():
        nseq = acs[acs.targets_recovered == tr].shape[0]
        upd_spn = [spn + i * j - j * (nseq - 1) / 2 for i in range(nseq)] 
        ampl_counts.loc[(ampl_counts.species_no == spn) & 
                        (ampl_counts.targets_recovered == tr), 'species_no'] = upd_spn
            
ampl_counts['sample_type'] = ampl_counts.s_Sample.str.startswith('anopheles').replace({True:'ref',False:'seq'})

fig, ax = plt.subplots(1,1,figsize=(10,6))
ax.grid(b=True, which='major', color='lightgrey', linewidth=1.0)
ax.set_axisbelow(True)
sns.scatterplot(data = ampl_counts,
                x='species_no',
                y='targets_recovered',
                hue='group',
                style='sample_type',
#                 x_jitter=1, 
                hue_order=('Pyretophorus series','Neocellia series',
                         'Myzomyia series','Neomyzomyia series',
                         'Anopheles subgenus','Nyssorhynchus subgenus',
                         'Kerteszia subgenus'),
                ax=ax)

ax.set_xticks(range(len(all_spp)))
ax.set_xticklabels(all_spp, rotation=90);
# sns.swarmplot(x='species', 
#               y='targets_recovered',
#               hue='group',
#               marker='v',
#               hue_order=('Pyretophorus series','Neocellia series',
#                          'Myzomyia series','Neomyzomyia series',
#                          'Anopheles subgenus','Nyssorhynchus subgenus',
#                          'Kerteszia subgenus'),
#               data=ampl_counts[ampl_counts.s_Sample.str.startswith('anopheles')], ax=ax)
# sns.swarmplot(x='species', 
#               y='targets_recovered',
#               hue='group',
#               marker='v',
#               hue_order=('Pyretophorus series','Neocellia series',
#                          'Myzomyia series','Neomyzomyia series',
#                          'Anopheles subgenus','Nyssorhynchus subgenus',
#                          'Kerteszia subgenus'),
#               data=ampl_counts[~ampl_counts.s_Sample.str.startswith('anopheles')], ax=ax)
# plt.xticks();<jupyter_output><empty_output><jupyter_text>Heterozygosity<jupyter_code># alleles per sample-target combination
ac = miseq.groupby(['s_Sample','target'])['consensus'].nunique()
# heterozygosity as more than one allele
(ac > 1).sum() / ac.shape[0]
(ac == 0).sum(), (ac == 1).sum(), (ac == 2).sum(), (ac > 2).sum(), ac.shape[0]
# alleles per sample-target combination
ac = miseq.groupby(['s_Sample','target'])['consensus'].nunique()
# allele counts
ac = ac.unstack()
mean_ac = ac.mean()
std_ac = ac.std()
# combine all multiallelic calls
ac[ac > 2] = 3
# count values
ac = ac.apply(pd.Series.value_counts).T
# normalise
# ac = ac.apply(lambda x: x / x.sum(), axis=1)
# plot
ac.columns='Homozygous Heterozygous Multiallelic'.split()
ac.plot.bar(stacked=True, figsize=(12,4), width=0.9);
plt.title('Allelism across amplicons')
plt.legend(bbox_to_anchor=(1.1, 1.05))

plt.show()
# set up exon-intron plotting settings
exon_c = 'orangered'
intron_c = 'gold'
inter_c = 'royalblue'
# exon vs intergenic
annot = (panel_data.intergenic == 'Yes').replace({True:inter_c,False:exon_c})
annot.loc[panel_data.intron == 'Yes'] = intron_c
annot = annot.to_list()
import matplotlib.patches as mpatches
red_patch = mpatches.Patch(color=exon_c, label='Exon')
green_patch = mpatches.Patch(color=intron_c, label='Intron')
blue_patch = mpatches.Patch(color=inter_c, label='Intergenic')
fig, ax = plt.subplots(1,1, figsize=(12,4))
b = ax.bar(range(62), height=mean_ac, yerr=std_ac, color=annot);
ax.set_xlabel('Amplicon')
ax.set_ylabel('Allele count')
ax.legend(handles=[red_patch,green_patch,blue_patch], loc=2)
ax.grid()
plt.xticks(range(62), rotation=45);<jupyter_output><empty_output><jupyter_text>Explore alignments<jupyter_code>import Bio
# read alignments and distance matrices
aln_data = dict()
for ampl in AMPLS:
    aln = AlignIO.read(open(ALN_ALL.format(ampl)), 'fasta', alphabet=Bio.Alphabet.generic_dna)
    for record in aln:
        record.seq = record.seq.upper()
    aln_data[ampl] = aln<jupyter_output><empty_output><jupyter_text>Gaps in alignment vs number of sequences<jupyter_code># how gappy are alignments?
ngaps = list()
for ampl in AMPLS:
    a = aln_data[ampl]
    ngaps.append([str(s.seq).count('-') for s in a])
# estimate means and stds
mean_gaps = [np.mean(g) for g in ngaps]
std_gaps = [np.std(g) for g in ngaps]
nseq = [len(g) for g in ngaps]
fig, ax = plt.subplots(1,1, figsize=(12,4))
b = ax.bar(range(62), height=[-g for g in mean_gaps], yerr=std_gaps, color=annot);
ns = ax.bar(range(62), height=nseq, color=annot);
ax.set_xlabel('Amplicon')
ax.set_ylabel('# gaps{}# seq'.format(' '*30))
ax.axhline(0,c='k')
ax.set_yticklabels([200,100,0,100,200])
# ax.grid()
ax.legend(handles=[red_patch,green_patch,blue_patch], loc=2)
plt.xticks(range(62), rotation=45);<jupyter_output><empty_output><jupyter_text>Sequence diversity<jupyter_code># collect dists
dists = list()
for ampl in AMPLS:
    a = aln_data[ampl]
    # need for blastn and trans DNA models to work
#     for record in a:
#         record.seq = str(record.seq).replace('N','-')
    l = a.get_alignment_length()
    # get distance matrix
    dm = DistanceCalculator('identity').get_distance(a).matrix
    # flatten
    d = [item for sublist in dm for item in sublist]
    # remove zeroes, multiply by lengths
    d = [i * l for i in d if i > 0]
    dists.append(d)
# get mean and sd
mean_dists = [np.mean(d) for d in dists]
std_dists = [np.std(d) for d in dists]
fig, ax = plt.subplots(1,1, figsize=(12,4))
b = ax.bar(range(62), height=mean_dists, yerr=std_dists, color=annot);
ax.set_xlabel('Amplicon')
ax.set_ylabel('Pairwise dists')<jupyter_output><empty_output><jupyter_text>combined plot<jupyter_code>from matplotlib import lines

fig, (ax1, ax2, ax3, ax4) = plt.subplots(4,1, 
                                    figsize=(12,10), 
                                    constrained_layout=True,
                                    gridspec_kw={'height_ratios':(2,1,1,1)})

w = 1
# split chromosomes
for ax in (ax2, ax3, ax4):
    ax.axvline(11.5, c='k', lw=2, alpha=.4)
    ax.axvline(28.5, c='k', lw=2, alpha=.4)
    ax.axvline(38.5, c='k', lw=2, alpha=.4)
    ax.axvline(55.5, c='k', lw=2, alpha=.4)
# subplot 1
sns.scatterplot(data = ampl_counts,
                x='species_no',
                y='targets_recovered',
                hue='group',
                style='sample_type',
                hue_order=('Pyretophorus series','Neocellia series',
                         'Myzomyia series','Neomyzomyia series',
                         'Anopheles subgenus','Nyssorhynchus subgenus',
                         'Kerteszia subgenus'),
                ax=ax1)

ax1.set_xticks(range(len(all_spp)))
ax1.set_xticklabels(all_spp, rotation=90)
ax1.grid(b=True, which='major', color='lightgrey', linewidth=1.0)
ax1.set_axisbelow(True)
ax1.set_title('Amplicon recovery across taxonomy', loc='left')
ax1.set_ylabel('Amplicons recovered')
ax1.legend(bbox_to_anchor=(1, 0.9), frameon=False)
# subplot 2
ac.plot.bar(stacked=True, width=w, ax=ax2);
ax2.set_ylabel('Number of samples')
ax2.legend(bbox_to_anchor=(1, 0.9), frameon=False)
ax2.set_xlabel('')
ax2.set_xticklabels(range(62), rotation=90)
ax2.set_title('Allelism and sample recall across amplicons', loc='left')
# subplot 3
ax3.bar(range(62), width=w, height=nseq, color=annot);
ax3.set_xlabel('Amplicon')
ax3.set_xticks(range(62))
ax3.set_xticklabels(range(62), rotation=90)
ax3.set_xlabel('')
ax3.set_ylabel('Sequences')
ax3.legend(handles=[red_patch,green_patch,blue_patch], bbox_to_anchor=(1, 0.9), frameon=False)
ax3.set_title('Unique sequences across amplicons', loc='left')
# subplot 4
ax4.bar(range(62), height=[g+10 for g in mean_gaps], bottom=-10, yerr=std_gaps, color=annot, width=w);
ax4.set_xlabel('Amplicon')
ax4.set_xticks(range(62))
ax4.set_xticklabels(range(62), rotation=90)
ax4.set_ylabel('Gaps in alignment')
ax4.set_xlabel('')
ax4.legend(handles=[red_patch,green_patch,blue_patch], bbox_to_anchor=(1, 0.9), frameon=False)
ax4.set_title('Alignment gaps across amplicons', loc='left')
# shared axes
ax2.get_shared_x_axes().join(ax2, ax3, ax4)
ax2.autoscale(tight=True)
# labels
for label, ax, yloc in zip('A B C D'.split(), (ax1, ax2, ax3, ax4), (1.08, 1.2, 1.2, 1.2)):
    ax.text(-0.02, yloc, label, transform=ax.transAxes,
          fontsize=16, fontweight='bold', va='top', ha='right')

    
plt.savefig('data/fig2.pdf')
plt.show()<jupyter_output><empty_output><jupyter_text>Combine sample metadata with sequencing stats<jupyter_code>comb_meta = meta.set_index('s_Sample')
comb_meta['Total Reads'] = seq.groupby('s_Sample')['reads'].sum()
comb_meta['Amplicons Recovered'] = seq.groupby('s_Sample')['target'].nunique()
# using count instead of nunique to account for perfect repetitions in reference genomes
ac = seq.groupby(['s_Sample','target'])['consensus'].count().unstack()
comb_meta['Homozygous Calls'] = (ac == 1).sum(axis=1)
comb_meta['Heterozygous Calls'] = (ac == 2).sum(axis=1)
comb_meta['Multiallelic Calls'] = (ac > 2).sum(axis=1)
comb_meta.to_csv(COMB_SEQ_META)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15594495119219076, 'fraction_numerical': 0.018803008481357016, 'mean_word_length': 4.258308792595709, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 50, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '175771', 'license': 'MIT', 'n_tokens_mistral': 4820, 'n_tokens_neox': 4530, 'n_words': 807, 'path': 'work/7_species_id/1_explore_data.ipynb', 'repo_name': 'malariagen/ano-spp', 'size': '12498'}