<jupyter_start><jupyter_text>Alignment handling and analysis using the Alignment class Contents:1. Alignment input and output2. Selecting rows and columns3. Calculating statistics of alignments4. Learning exercises to test your skillsAuthors: <NAME>, <NAME><jupyter_code>import numpy as np
from evcouplings.align import Alignment, map_matrix<jupyter_output><empty_output><jupyter_text>1. Alignment input and output<jupyter_code>## Read in a sequence alignment from a fasta file
with open("example/example_aln.a2m", "r") as infile:
    aln = Alignment.from_file(infile, format="fasta")
# Sequence length and number of sequences
print(f"alignment is of length {aln.L} and has {aln.N} sequences")
# Save the alignment to a file
with open("example/_temporary_file.fasta", "w") as outfile:
    aln.write(fileobj=outfile, format="fasta")<jupyter_output><empty_output><jupyter_text>You can also define a sequence alignment using a dictionary of identifiers that point to sequences, as follows<jupyter_code>d = {
    "seq1": "HELLO",
    "seq2": "WORLD"
    }
aln2 = Alignment.from_dict(d)<jupyter_output><empty_output><jupyter_text>2. Selecting rows and columns<jupyter_code># aln.matrix returns a numpy matrix of the sequences in the alignment (rows = sequences, columns=positions)
# This matrix can be indexed into as any numpy matrix, with indexing beginning at 0

#sequence 10, positions 1-10 
aln.matrix[9,0:10]

# The matrix can be indexed by looking up the index of an identifier of interest

# view the sequence identifiers
print("all the ids:", aln.ids, "\n")

# Find the index of the sequence PCDBI_MOUSE
index = aln.id_to_index["PCDBI_MOUSE"]
print(f"sequence PCDBI_MOUSE has index {index}\n")

# select the sequence of interest
sequence = aln.matrix[index,:]
print("the sequence is", sequence)
# aln.matrix can also be indexed with getitem-style access
aln["PCDBI_MOUSE"]
# Create a new sequence alignment with just the sequences of interest
ids_to_select = ['Q641K6_MOUSE','B9EJ92_MOUSE','Q91Y16_MOUSE']
selection_index = [aln.id_to_index[x] for x in ids_to_select]
aln_subsection = aln.select(sequences=selection_index)
print(f"the new alignment has {aln_subsection.N} sequences")
# Say the first 50 positions contain a domain of interest. 
# You can also select only certain columns

aln_subsection = aln.select(columns=range(0,50))
print(f"the new alignment has {aln_subsection.L} positions")<jupyter_output>the new alignment has 50 positions<jupyter_text>3. Calculating statistics of alignments<jupyter_code># The sequence alignment can be converted into a matrix of integers for faster handling
aln.matrix_mapped = map_matrix(aln.matrix, aln.alphabet_map)
print("mapping indices: ",aln.alphabet_map)
# Display the frequencies of each amino acid in position 21
frequency_21 = aln.frequencies[20,:]
print(f"There are {frequency_21.size} amino acids in the alignment")

for idx, amino_acid in enumerate(aln.alphabet):
    print(f"{amino_acid}, {np.round(frequency_21[idx], 2)}")

#calculates conservation of each positions using the column entropy
cons = aln.conservation()[20]
print(f"position 21 is {np.round(cons,2)} percent conserved")
# calculate the percent identity of every sequence in the alignment to the first sequence
aln.identities_to(aln.matrix[0])<jupyter_output><empty_output><jupyter_text>4. Learning exercises and solutions<jupyter_code>## Some learning exercises

# Which sequence in the example alignment contains the most gaps?

# Calculate the percent identity of all sequences to all other sequences and plot as a heatmap

# Create a new alignment of 100 random sequences, where each position is drawn from the frequency distribution of the natural sequences in the example alignment at that site. 

----Scroll down for solutions----<jupyter_output><empty_output><jupyter_text>Solutions<jupyter_code># Which sequence contains the most gaps?

#use the "count" method of the class
maximum1 = aln.count(axis="seq",char="-").argmax()

#or, use the numpy sum
maximum2 = (aln.matrix=="-").sum(axis=1).argmax()

assert maximum1==maximum2

print(aln.ids[maximum1])
# Calculate the percent identity of all sequences to all other sequences and plot as a heatmap

import seaborn as sns
import matplotlib.pyplot as plt

identity_matrix = np.zeros((aln.N, aln.N))
identity_matrix.fill(np.nan)

#loop through each sequence
for idx in range(identity_matrix.shape[0]):

    #calculate all the percent identitites to that sequence
    idents = aln.identities_to(aln.matrix[idx,:])
    identity_matrix[idx] = idents

#plot the heatmap
fig = plt.figure(figsize=(11,10))
ax = fig.gca
sns.heatmap(identity_matrix,cmap="Blues", xticklabels=aln.ids, yticklabels=aln.ids)
# Create a new alignment of 100 random sequences, where each position is drawn from the frequency distribution of the natural sequences at that site.  
from collections import defaultdict

#Slow, loopy way

aln_dict = defaultdict(str)

for i in range(0,100):
    #name for fake sequence
    seq_name = "rand_"+str(i)
    
    sequence = np.zeros((1,aln.L),dtype=str)
    # for each position
    for col_index in range(aln.L):
        
        # choose a letter 
        letter = np.random.choice(
            list(aln.alphabet),
            1, 
            replace=False,
            p = aln.frequencies[col_index,:]
        )[0]
        #save to new sequence
        sequence[0,col_index] = letter
        
    # save sequence to alignment dictionary
    aln_dict[seq_name] = "".join(list(sequence.reshape(-1)))

rand_aln = Alignment.from_dict(aln_dict)

# verify that the two alignments have similar residue frequencies

print("highest difference in frequency", abs(rand_aln.frequencies - aln.frequencies).max())
    
# faster way!
n_seqs = 10000
sequence_matrix = np.zeros((n_seqs,aln.L),dtype=str)

#loop through columns
for i in range(aln.L):
    
    # make n_seqs random choices
    letters = np.random.choice(
            list(aln.alphabet),
            n_seqs, 
            replace=True,
            p = aln.frequencies[i,:]
        )

    #save the new columns
    sequence_matrix[:,i] = letters
    
aln_dict = {f"rand_{i}":"".join(list(x)) for i,x in enumerate(sequence_matrix)}
rand_aln = Alignment.from_dict(aln_dict)

#verify similar frequencies
print("highest difference in frequency", abs(rand_aln.frequencies - aln.frequencies).max())<jupyter_output>highest difference in frequency 0.01198113207547169
{'fraction_non_alphanumeric': 0.09643638637074085, 'fraction_numerical': 0.016567677399187245, 'mean_word_length': 4.770063119927863, 'pattern_counts': {'":': 3, '<': 29, '<?xml version=': 0, '>': 29, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '286432', 'license': 'MIT', 'n_tokens_mistral': 2106, 'n_tokens_neox': 2007, 'n_words': 638, 'path': 'notebooks/alignment_analysis.ipynb', 'repo_name': 'berkalpay/EVcouplings', 'size': '6413'}