Il y a habituellement plus d’une façon d’accomplir une tâche donnée. Certaines sont meilleures que d’autres, plusieurs sont équivalentes. Décider laquelle utiliser dépend bien souvent du temps de calcul, de la facilité d’utilisation et/ou de nos préférences et compétences personnelles.

Supposons que j’aie une matrice contenant des milliers de positions pour un chromosome donné avec les colonnes suivantes : Nom_de_l’élément, Début, Fin. Toutes les positions se rapportent à un même chromosome et la taille des éléments est variable. Pour une liste de SNVs donnée, je voudrais savoir lesquels chevauchent mes éléments chromosomiques.

Pour ce faire, je pourrais construire deux fichiers bed/gtf et utiliser dans ma console la commande intersectBed de la suite d’outils bedtools. Cependant, comme je suis une grande fan de R, c’est en R que j’ai élaboré mes trois façons.

1. Les boucles for imbriquées

Voici bien souvent le choix du débutant. On parcourt tous les SNVs (première boucle for) et on compare chacun à tous les éléments chromosomiques pour voir si l’un de ces éléments chevauchent le SNV. Si au moins un élément le chevauche, on retourne true. Il n’y a rien de mal avec cette méthode. Excepté que R n’est pas le champion des boucles imbriquées. Ce sera très long à rouler pour un gros jeu de données.
C’est la solution facile, mais vous deviendrez rapidement impatient si vous avez besoin de visualiser vos résultats rapidement.


for (i in nrow(snvs)) {
bool <- c()
 for (j in nrow(features)){
    b <- snvs[i, "position"] > features[j, "Start"] & snvs[i, "position"] < features[j, "Stop"]	 	 
    bool <- c(bool, b)  	 
 }
 snvs[i,'found'] <- any(bool)
}

2. L'utilisation de apply de apply

En R, vous pouvez toujours remplacer une boucle for par la fonction apply qui est beaucoup plus rapide. Pour notre exemple, nous pouvons parcourir les SNVs dans un premier apply et faire la comparaison dans un second. Nous avons aussi des apply imbriqués. Sur Linux ou Mac, vous pouvez même utiliser le package parallel pour rouler sur plusieurs processeurs. Cette solution est beaucoup plus rapide, mais beaucoup moins évidente en terme de code à écrire.



compare <- function(p, x) {return (x > as.numeric(p['start']) & x < as.numeric(p['end']))}	 	
search <- function(x, features){	 	
 b <- apply(features, 1, compare, x=x)	 	 
 return (any(b))	 	
}	
 	
rres <- cbind (as.matrix(snvs), found=apply(snvs, 1, search, features=features))

3. Le package IRanges de Bioconductor

La dernière façon de procéder que je veux présenter profite du talent et des efforts des autres développeurs de la communauté R. Ce n'est pas déplaisant de laisser les autres travailler pour nous!
Le package IRanges implémente des fonctions pour manipuler des intervalles génomiques (ou seulement des intervalles). Les fonctions sont relativement simples à comprendre et à utiliser et elles sont rapides. Donc, avant de réinventer la roue, il est utile de faire une petite recherche web (google!) pour vérifier s'il n'existerait pas déjà un outil pouvant nous aider.


ir = IRanges (start = starts, end=stops)
sr = IRanges (start = snvs, end = snvs)
findOverlaps(sr, ir)