Commit 29e1ab0e authored by Theo Serralta's avatar Theo Serralta

Remove files from git

parent 2fae1bf6
import pysam
samfile = pysam.AlignmentFile("/home/theo/dev/git/cnvCallerGPU/3544.bam", "rb")
align = pysam.AlignmentFile("align.bam", "wb", header = samfile.header)
for read in samfile:
if read.is_mapped:
align.write(read)
samfile.close()
align.close()
total_sequences_alignées = 0
align_file = pysam.AlignmentFile("/home/theo/dev/git/cnvCallerGPU/align.bam", "rb")
reference = align_file.references
reference_dict = {reference: i for i, reference in enumerate(reference)}
print("Index\tChromosome")
for chromosome, index in reference_dict.items():
print(f"{index}\t{chromosome}")
for read in align_file:
print(read)
chr_name = reference[read.reference_id]
position = read.reference_start + 1
print(f"Chromosome: {chr_name}")
print(f"Position: {position}")
print(f"Non aligné: {read.is_unmapped}")
total_sequences_alignées += 1
print(f"Total de sequences alignées : {total_sequences_alignées}")
align_file.close()
import pysam
def calcul_depth_seq(bamfile, seq, window_size=100):
samfile = pysam.AlignmentFile(bamfile, "rb")
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): #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): #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:#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_all(bamfile_path, window_size=100)
align_file = "/home/theo/dev/git/cnvCallerGPU/test_depth.txt"
depth_sup_mean = "/home/theo/dev/git/cnvCallerGPU/depth_sup_mean.txt"
def calc_depth_mean(align_file):
with open(align_file, 'r') as file:
total_depth = 0
total_window = 0
for line in file:
depth_str_split = line.split()
depth_str_value = depth_str_split[6:]
total_depth += float(depth_str_value[0])
total_window += 1
mean = total_depth / total_window
return mean
with open(depth_sup_mean, 'w') as output:
with open(align_file, 'r') as file:
mean = calc_depth_mean(align_file)
for line in file:
depth_str_split = line.split()
depth_str_value = depth_str_split[6:]
if float(depth_str_value[0]) > mean:
output.write(f"{depth_str_split[5]}: {depth_str_value[0]}\n")
#import pysam
#seq_inf= {}#On definit le dictionnaire
#samfile = pysam.AlignmentFile("/home/theo/dev/git/cnvCallerGPU/3544.bam", "rb")#Importation du fichier bam
#for read in samfile.fetch(until_eof=True):#Boucle qui va lire le fichier bam en entier
#seq_id=read.query_name#Attribution de l'ID de la sequence
#chro=read.reference_name#Attribution du nom du chromosome
#start=read.reference_start#Attribution position départ
#lenght=read.query_length#Attribution longueur sequence
#end=read.reference_end#Attribution position fin
#seq_inf[chro]={}#Création d'un nouveau dictionnaire avec le nom du chromosome à chaque itération de boucle
#seq_inf[chro]=({"seq_id" : seq_id, "start" : start, "end" : end, "lenght" : lenght, "depth" : 0})#attribution de différentes informations pour chaque chromosome à chaque itération de boucle
#samfile.close()
#def depth_function(seq_inf):#Definition d'une fonction qui va calculer la profondeur
# for chro, data in seq_inf.items():
# depth = [0] * (data[-1]["end"] +1) #on crée une liste qui stock la profondeur de chaque chromosome. On prend la dernière position lue de chaque dictionnaire. On rajoute +1 parce que ça compte à partir de 0.
# for position in range(start, end +1):#on reprend les position start et end défini au début
# depth[position] += 1 #On calcule la profondeur en ajoutant +1 à chaque itération
#depth_function(seq_inf)#On fait appel à la fonction pour le dictionnaire de base
#import pyopencl as cl
#ctx = cl.create_some_context()
#queue = cl.CommandQueue(ctx)
#mf = cl.mem_flags
import pysam
seq_inf = {} # On définit le dictionnaire
samfile = pysam.AlignmentFile("/home/theo/dev/git/cnvCallerGPU/3544.bam", "rb") # Importation du fichier bam
for read in samfile.fetch(until_eof=True): # Boucle qui va lire le fichier bam en entier
seq_id = read.query_name # Attribution de l'ID de la séquence
chro = read.reference_name # Attribution du nom du chromosome
start = read.reference_start # Attribution position de départ
length = read.query_length # Attribution longueur de la séquence
end = read.reference_end # Attribution position de fin
if chro not in seq_inf:
seq_inf[chro] = [] # Création d'une liste vide pour chaque chromosome
seq_inf[chro].append({"seq_id": seq_id, "start": start, "end": end, "length": length, "depth": 0})
samfile.close()
def depth_function(seq_inf): # Définition d'une fonction qui va calculer la profondeur
for chro, data in seq_inf.items():
for record in data:
start = record.get("start") # Utiliser get() pour obtenir la valeur de "start"
end = record.get("end") # Utiliser get() pour obtenir la valeur de "end"
if start is not None and end is not None:#Comme j'ai des None on calcule la profondeur quand on peut...
depth = [0] * (end - start + 1) # On crée une liste qui stocke la profondeur pour chaque position du chromosome
for position in range(start, end + 1): #on reprend les position start et end défini au début
depth[position - start] += 1 # On calcule la profondeur en ajoutant +1 à chaque iteration
depth_function(seq_inf) # On fait appel à la fonction pour le dictionnaire de base
print(seq_inf[chro])
import pysam
samfile = pysam.AlignmentFile("/home/theo/dev/git/cnvCallerGPU/3544.bam", "rb")
non_align = pysam.AlignmentFile("non_align.bam", "wb", header = samfile.header)
for read in samfile:
if read.is_unmapped:
non_align.write(read)
samfile.close()
non_align.close()
total_sequence_non_alignées = 0
non_align_file = pysam.AlignmentFile("/home/theo/non_align.bam", "rb")
reference = non_align_file.references
reference_dict = {reference: i for i, reference in enumerate(reference)}
print("Index\tChromosome")
for chromosome, index in reference_dict.items():
print(f"{index}\t{chromosome}")
for read in non_align_file:
print(read)
chr_name = reference[read.reference_id]
position = read.reference_start + 1
print(f"Chromosome: {chr_name}")
print(f"Position: {position}")
print(f"Non aligné: {read.is_unmapped}")
total_sequence_non_alignées += 1
print(f"Total de sequences non alignées : {total_sequence_non_alignées}")
non_align_file.close()
import pysam
nb_seq=0
samfile = pysam.AlignmentFile("/home/theo/Documents/alternance/Projet Alternance/dijex8000.bam", "rb")
for read in samfile.fetch("chrY"):
print(read)
nb_seq+=1
print(nb_seq)
samfile.close()
#for pileupcolumn in samfile.pileup("chr1", 1):
# print("\ncoverage at base %s = %s" % (pileupcolumn.pos, pileupcolumn.n))
# for pileupread in pileupcolumn.pileups:
# if not pileupread.is_del and not pileupread.is_refskip:
# query position is None if is_del or is_refskip is set.
# print('\tbase in read %s = %s' %
# (pileupread.alignment.query_name, #pileupread.alignment.query_sequence[pileupread.query_position]))
#samfile.close()
import numpy as np
import pyopencl as cl
a_np = np.random.rand(50000).astype(np.float32)
b_np = np.random.rand(50000).astype(np.float32)
ctx = cl.create_some_context()#C'est la "table de jeu". Permet aux devices de recevoir les kernels et de transferer les données
queue = cl.CommandQueue(ctx)#Permet de créer la queue de commandes dans le context
mf = cl.mem_flags #specifie comment la mémoire doit être gérée et utilisée
#print(mf.to_string(1))
a_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=a_np)#espace de mémoire temporaire. Ce sera une allocation de mémoire sur un periph (gpu) pour stocker les données qui seront utilisées dans des calculs parallèles
b_g = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=b_np)
#creation du programme
prg = cl.Program(ctx, """
__kernel void sum(
__global const float *a_g, __global const float *b_g, __global float *res_g)
{
int gid = get_global_id(0);
res_g[gid] = a_g[gid] + b_g[gid];
}
""").build()
res_g = cl.Buffer(ctx, mf.WRITE_ONLY, a_np.nbytes)#Crée un buffer dans le contexte allouant de la mémoire sur le gpu pour stocker les resultats du calcul. Comme "WRITE-ONLY", il est destiné à stocker les données de sortie du noyau OpenCl
knl = prg.sum#On defini knl afin d'appeler le noyau sum du programme qu'on a créé précédement
knl(queue, a_np.shape, None, a_g, b_g, res_g)#execute le noyau sum sur le GPU à l'aide de la queue
res_np = np.empty_like(a_np)#Crée un tableau numpy de même dimension que a_np mais vide afin de stocker les résultats du calcul
cl.enqueue_copy(queue, res_np, res_g)#Copie les données depuis le buffer sur le petipherique (GPU) vers le tableau res_np
# Check on CPU with Numpy:
print(res_np - (a_np + b_np))#Permet de voir si les donées ont été exportées correctement dans le tableau res_np en checkant la difference entre le tableau et les deux variables du départ. Ca doit etre egal à 0.
print(np.linalg.norm(res_np - (a_np + b_np)))#Calcul la norme L2 de la difference entre res_np et (a_np + b_np). Si les deux vecteurs sont identique, ça devrait etre egal à 0. Ca nous donne une mesure de l'erreur entre les resultats du noyau OpenCl et les resultats attendus
assert np.allclose(res_np, a_np + b_np)#Compare les tableaux et verifie s'ils sont suffisament proches les uns des autres. S'ils le sont, le programme tourne normalement sinon, l'assertion renvoie un message d'echec.
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