# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # Biopython - Esercizio4

# [MAFFT](https://www.ebi.ac.uk/Tools/msa/mafft/) è un tool di allineamento multiplo sviluppato da EMBL-EBI (European Bioinformatics Institute - European Molecular Biology Laboratory) per sequenze di DNA.
#
# Usare MAFFT (scegliendo ClustalW come formato di output) per allineare i 14 genomi completi di SARS-CoV-2 presenti nel file `covid-sequences.fasta` sequenziati nel novembre 2021 e scaricati dal sito di [NCBI](https://www.ncbi.nlm.nih.gov/sars-cov-2/). Il primo, con identificatore `NC_045512.2`, è il genoma di riferimento.
#
# Trovare tutte le variazioni rispetto al genoma di riferimento.
#
# ---
#
# **Variazione**: posizione della colonna di allineamento in cui esiste almeno un genoma che ha mismatch con quello di riferimento.
#
# Esempio di allineamento con variazioni in posizione 8 e 13:
#
#     REF   AAGCTGATTGCACGC-T
#     G1    --GCAGAGTGCAGGCCT
#     G2    --GCCGAGTGCACGCCT
#
# **Variazione 5**: `T` nel reference e `A` in G1 e `C` in G2.
#
# **Variazione 8**: `T` nel reference e `G` sia in G1 e G2.
#
# **Variazione 13**: `C` nel reference e `G` in G1.
#
# **Variazione 16**: `-` nel reference e `C` sia in G1 che in G2.
#
# ---
#
# Si richiede di:
# - costruire il data frame delle variazioni in cui le colonne sono tutte le posizioni 1-based delle variazioni e le righe sono indicizzate con l'identificatore del genoma.
# - estrarre il genoma con più variazioni e quello con meno variazioni
# - ottenere il data frame delle variazioni "complete", cioè in cui tutti i genomi variano rispetto al riferimento.
# - produrre il data frame delle variazioni "stabili" in cui tutti i genomi variano allo stesso modo rispetto al riferimento. 
# - ottenere la lista delle posizioni in cui c'è un gap nel genoma di riferimento.
# - ottenere la lista delle posizioni in cui c'è un gap in almeno uno dei genomi (diversi dal riferimento)

# Installare il package `Bio` di Biopython.

# + active=""
# !conda install -y -c conda-forge biopython
# -

# Importare il package `Bio`.

import Bio

# Importare il package `AlignIO` che è il package per manipolare file contenenti allineamenti multipli in diversi formati (tra cui `clustal` che è quello del file di input).

from Bio import AlignIO

# #### Leggere l'allineamento in input

# Il package `AlignIO` mette a disposizione la funzione `read` per leggete un allineamento:
#
#        AligIO.read(input_file_name, format)
#        
# e restituisce un oggetto di tipo `MultipleSeqAlignment` che è un oggetto iterabile contenente oggetti `SeqRecord`, uno per ognuna delle righe dell'allineamento letto.

alignment = AlignIO.read("mafft-alignments.clustalw", "clustal")

# La lunghezza dell'allineamento in input (numero di colonne della matrice di allineamento) è:

alignment.get_alignment_length()

# Trasformare l'oggetto in una lista di oggetti `SeqRecord`.

alignment = list(alignment)
alignment

# #### Eliminare i gap iniziali.
#
# Trovare il più lungo prefisso di soli simboli `-` delle righe dell'allineamento. Supponendo che tale prefisso sia lungo `g`, eliminare da ogni riga dell'allinemento il prefisso di lunghezza `g`.
#
# Ad esempio il seguente allineamento composto da tre righe:
#
#     GTATGTGTCATGTTTTTGCTA
#     --ATGTGTCATG-TTT-----
#     ----GTGTCATGTTTTTG---
#     
# presenta un più lungo prefisso di soli simboli `-` di lunghezza `g=4` (terza riga). Eliminando da tutte le righe un prefisso di lunghezza 4 si ottiene:
#
#         GTGTCATGTTTTTGCTA
#         GTGTCATG-TTT-----
#         GTGTCATGTTTTTG---

# +
import re

gap_list = [re.findall('^-+', str(row.seq)) for row in alignment]
gap_size_list = [len(gap[0]) for gap in gap_list if gap]
gap_size_list[:0] = [0]
leading_gaps = max(gap_size_list)
alignment = [row[leading_gaps:] for row in alignment]
# -

alignment

# #### Eliminare i gap finali.
#
# Trovare il più lungo suffisso di soli simboli `-` delle righe dell'allineamento. Supponendo che tale suffisso sia lungo `g`, eliminare da ogni riga il suffisso di lunghezza `g`.
#
# Ad esempio il seguente allineamento composto da tre righe:
#
#         GTGTCATGTTTTTGCTA
#         GTGTCATG-TTT-----
#         GTGTCATGTTTTTG---
#         
# presenta un più lungo suffisso di soli simboli `-` di lunghezza `g=5` (seconda riga). Eliminando da tutte le righe un suffisso di lunghezza 5 si ottiene:
#
#         GTGTCATGTTTT
#         GTGTCATG-TTT
#         GTGTCATGTTTT

gap_list = [re.findall('-+$', str(row.seq)) for row in alignment]
gap_size_list = [len(gap[0]) for gap in gap_list if gap]
gap_size_list[:0] = [0]
trailing_gaps = max(gap_size_list)
alignment = [row[:len(row)-trailing_gaps] for row in alignment]

alignment

# #### Creare la lista degli identificatori dei genomi

index_list = [row.id for row in alignment]

index_list

# #### Creare il dizionario contenente i dati per costruire il data frame
#
# - `key`: posizione 1-based della variazione (posizione della colonna nell'allineamento in input)
#
# - `value`: lista dei simboli allineati coinvolti nella variazione (il primo simbolo deve essere quello del reference, mentre se un genoma non presenta una differenza con il reference si deve inserire la stringa vuota)

# +
df_data = {}

reference = alignment.pop(0)

for (i,c) in enumerate(reference):
    variant_list = []
    is_variant = False
    for row in alignment:
        variant = ''
        if row[i] != c and row[i] in {'A', 'C', 'G', 'T'}:
            is_variant = True
            variant = row[i]
            
        variant_list.append(variant)
        
    if is_variant:
        variant_list[:0] = [c]
        df_data[str(i+leading_gaps+1)] = variant_list
# -

df_data

# #### Creare il data frame
#
#     df = pd.DataFrame(df_data, index = index_list)

# +
import pandas as pd

df = pd.DataFrame(df_data, index = index_list)
# -

df

# #### Estrarre il genoma con più variazioni e quello con meno variazioni

# Determinare la lista del numero di variazioni per genoma (per tutti i genomi tranne quello di riferimento).

variants_per_genome = [len(list(filter(lambda x: x!='', list(row)))) for row in df.values]

variants_per_genome.pop(0)
variants_per_genome

# In alternativa:

variants_per_genome = [df.shape[1]-list(df.loc[index]).count('') for index in index_list[1:]]

variants_per_genome

# Estrarre il genoma con più variazioni.

index_list[variants_per_genome.index(max(variants_per_genome))+1]

# Estrarre il genoma con meno variazioni.

index_list[variants_per_genome.index(min(variants_per_genome))+1]

# In alternativa, per estrarre il genoma con meno variazioni:

null_df = pd.DataFrame((df == '').sum(axis=1), columns=['difference'])
null_df[1:][null_df[1:]['difference'] == null_df[1:]['difference'].max()]

# In alternativa, per estrarre il genoma con più variazioni:

null_df[1:][null_df[1:]['difference'] == null_df[1:]['difference'].min()]

# #### Determinare il data frame delle variazioni "complete"
#
# Selezionare dal data frame precedente le sole colonne relative a variazioni "complete".

df_complete = df[[col for col in df.columns if all(df[col] != '')]]

df_complete

# #### Determinare il data frame delle variazioni "stabili"
#
# Selezionare dal data frame precedente le sole colonne relative a variazioni "stabili".

df_stable = df_complete[[col for col in df_complete.columns if len(df_complete[col][1:].unique()) == 1]]

df_stable

# #### Ottenere la lista delle posizioni in cui c'è un gap nel genoma di riferimento.

ref_gaps = [col for col in df.columns if df[col][0] == '-']

ref_gaps

# #### Ottenere la lista delle posizioni in cui c'è un gap in almeno uno dei genomi (diversi dal riferimento).

other_gaps = [col for col in df.columns if any(df[col][1:] == '-')]

other_gaps

StarcoderdataJupyterScriptsDedupFiltered