Commit 5b182bf3 authored by Theo Serralta's avatar Theo Serralta

Add some comments and new code for all seq

parent 5c598d30
import pysam
def calcul_depth(bamfile, seq, window_size=100):
def calcul_depth_seq(bamfile, seq, window_size=100):
samfile = pysam.AlignmentFile(bamfile, "rb")
seq_length = samfile.lengths[samfile.references.index(seq)]
seq_length = samfile.lengths[samfile.references.index(seq)] #Permet d'obtenir la longueur de la sequence
for pos_start in range(0, seq_length, window_size):
pos_end = min(pos_start + window_size, seq_length)
window_depths = []
for pos_start in range(0, seq_length, window_size): #Itère sur la sequence par "fenêtre" en fonction de la taille
pos_end = min(pos_start + window_size, seq_length) #pos_start + window_size donne la position de début de la fenêtre, puis grâce à min() et seq_length je peux définir position de fin. le min me permet d'être certain que pos_end < seq_lenght. C-a-d que si la fenêtre dépasse la taille de la sequence, j'aurais quand même pos_end
window_depths = [] #liste pour stocker la profondeur de chaque fenêtre
for pileupcolumn in samfile.pileup(seq, start=pos_start, stop=pos_end):
window_depths.append(pileupcolumn.n)
for pileupcolumn in samfile.pileup(seq, start=pos_start, stop=pos_end): #itère sur chaque colone de pileup
window_depths.append(pileupcolumn.n)#ajoute le nb de lecture couvrant cette position à la liste window_depths
if window_depths:
if window_depths:#calcul la profondeur moyenne pour chaque sequence
avg_depth = sum(window_depths) / len(window_depths)
print(f"Profondeur moyenne pour la fenêtre {seq}:{pos_start}-{pos_end}: {avg_depth}")
samfile.close()
def calcul_depth_all(bamfile, window_size=100):
samfile=pysam.AlignmentFile(bamfile, "rb")
for seq in samfile.references:
calcul_depth_seq(bamfile, seq, window_size)
samfile.close()
bamfile_path = "/home/theo/dev/git/cnvCallerGPU/align.bam"
sequence_of_interest = "chr21"
calcul_depth(bamfile_path, sequence_of_interest, window_size=1000)
#sequence_of_interest = "chr21"
calcul_depth_all(bamfile_path, window_size=100)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment