Discussion :

[karaudrey88]

Bonjour,
Je doit comparer plusieurs algorithmes d'alignement. Je viens d'obtenir tous les fichiers et j'aimerais comparer ces fichiers SAM a un fichier SAM de référence, pour savoir si le même read est aligné au même endroit et possède la même CIGAR LINE. Quelqu'un aurait il une idée pour faire cette comparaison assez rapidement?
Merci

[trex7g2]

Ce que je verrais de plus simple serait de parser ta référence pour créer une hash contenant le nom de ton read, le CIGAR et les coordonnées. ensuite tu n'auras plus qu'à vérifier si, pour le même read, tu as les mêmes coordonnées et la même CIGAR.

Par contre attention je sais d'expérience que certains logiciels traitent différemment leurs CIGAR donc les résultats peuvent être surprenant.

[karaudrey88]

Merci. j'ai fait le hash de la manière suivante :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
 
open(FILE2,"$samBWA") || die ("Erreur d'ouverture de samBWA") ;
while (my $ligne2 = <FILE2>){
	@tmp2 = split("\t", $_);
	$qname2=$tmp2[0];
	$flag2=$tmp2[1];
	$rname2=$tmp2[2];
	$pos2=$tmp2[3];
	$mapQ2=$tmp2[4];
	$cigar2=$tmp2[5];
	${ $hash2{$qname2} }{$pos2} = $cigar2;
 
}
close (FILE2);

Et j'ai fait le 2e fichier de la même manière.Pour tester les clés j'ai vu qu'il existait une fonction existe mais je ne voit pas comment comparer les valeurs $pos2 et $cigar2 à celles de l'autre fichier.

[trex7g2]

Ce que je te propose c'est de changer très légèrement ta hash :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
 
${ $hash2{$qname2} }{0} = $pos2;
${ $hash2{$qname2} }{1} = $cigar2;

De cette manière, tu auras moins de problèmes je pense (toutes tes hash existeront ç partir du moment ou ta référence est complète). C'est aussi plus facile de faire des comparaisons : par exemple il m'est déjà arrivé d'avoir un logiciel qui considérait les extrémités des reads comme "potentiellement incorrecte" et n'hésitait pas a couper 2/3 nucléotides. Si tu tombes sur un cas comme ça tu pourras facilement tester la distance entre la position obtenue et la position de référence.

Ensuite tu fais la même chose avec ton fichier de test, et tu n'auras plus qu'a comparer une à une chaque hash.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
open(FILE1,"$samBWA") || die ("Erreur d'ouverture de samBWA") ;
while (my $ligne1 = <FILE1>){
	@tmp1 = split("\t", $_);
	$qname1=$tmp1[0];
	$flag1=$tmp1[1];
	$rname1=$tmp1[2];
	$pos1=$tmp1[3];
	$mapQ1=$tmp1[4];
	$cigar1=$tmp1[5];
	if ( $pos1 == ${ $hash2{$qname1} }{0} ){
		 # la position est bonne
	}

edit : des balises c'est bien