import pysam

from pysam import AlignmentFile

"""
Explanation: Esercizio1 - Pysam
Prendere in input un BAM file che contiene allineamenti ad un'unica reference (cromosoma X) e:

controllare se sono presenti paired-end reads
determinare il numero e le lunghezze di tutti gli introni supportati dagli allineamenti del BAM file
determinare la base della reference a copertura massima e produrre un file SAM contenente i soli allineamenti che coprono tale base

Importare il modulo pysam
End of explanation
"""


pysam.index('./sample.bam')
bamfile = AlignmentFile('./sample.bam', 'rb')

"""
Explanation: Caricare il file BAM in input
End of explanation
"""


alignment_list = list(bamfile.fetch())

any(alignment.is_paired for alignment in alignment_list)

"""
Explanation: Controllare se ci sono dei reads paired-end.
End of explanation
"""


spliced_cigars = [alignment.cigarstring for alignment in alignment_list if 'N' in alignment.cigarstring]

import re

intron_lengths = [re.findall('(\d+)N', cigar)[0] for cigar in spliced_cigars]

"""
Explanation: Determinare il numero e le lunghezze di tutti gli introni supportati dagli allineamenti
Determinare la lista delle lunghezze degli introni supportati con duplicati.
End of explanation
"""


set(intron_lengths)

"""
Explanation: Determinare il set delle lunghezze degli introni supportati.
End of explanation
"""


pileup_iter = bamfile.pileup()

pileup_columns = list(pileup_iter)

"""
Explanation: Trovare la base a copertura massima e produrre il SAM file contenente tutti gli allineamenti che coprono questa base.
Determinare la lista delle pileup columns.
End of explanation
"""


max_col = max(pileup_col.nsegments for pileup_col in pileup_columns)

max_pileup_col = [pileup_col for pileup_col in pileup_columns if pileup_col.nsegments == max_col][0]

max_pileup_col.nsegments

"""
Explanation: Estrarre la colonna di pileup con il maggior numero di reads.
End of explanation
"""


max_pileup_col.set_min_base_quality(0)

"""
Explanation: Settare a 0 la qualità minima della colonna di pileup (altrimenti le reads non vengono tenute perché hanno una qualità bassa).
End of explanation
"""


max_pileup_col.pileups

"""
Explanation: Estrarre la lista dei read di pileup relativi alla colonna di altezza massima.
End of explanation
"""


pileup_slignments = [pileup_read.alignment for pileup_read in max_pileup_col.pileups]

"""
Explanation: Estrarre la lista degli allineamenti relativi alla colonna di altezza massima.
End of explanation
"""


outputfile = pysam.AlignmentFile('./prova.sam', "w", template=bamfile)

for read in pileup_slignments:
    outputfile.write(read)
    
outputfile.close()

"""
Explanation: Produrre gli allineamenti in un SAM file utilizzando la stessa Header Section del BAM file.
End of explanation
"""


