Module

Warning Depuis l'été 2015, le cluster MBB et celui de l'ISE-M sont passés sous module pour gérer R (voir lien). Il faut donc penser à précéder votre ligne de commande par module load R-... dans votre script de soumission SGE (points de suspension correspondants à la version désirée). De même, les packages doivent être installés dans votre environnement utilisateur (home).

Utilisation générale de R

Après avoir chargé la version souhaitée avec module, il faudra lancer votre script de soumission avec l'option SGE -V pour exporter les variables d'environnement actuelles dans le job. Par ailleurs, dans votre script de soumission il faut préférer l'utilisation de R CMD BATCH my_script.R ou R --vanilla --slave --args [your args here] < my_script.R en lieu et place de Rscript my_script.R ref En effet, Rscript contient en dur le chemin d'installation du binaire R qui correspond à l'installation système.

Bonnes pratiques

Après avoir chargé R par module, il peut d'avérer intéressant d'utiliser le package renv. Ce dernier permet d'isoler par projet chaque environnement R (un peu de la même manière que pour pyenv/virtualenv en python). Utilisation :

qrsh
module load R/3.5.3
mkdir my_project
R

install.packages('renv')
library('renv')
renv::init()
q()

R

Exemple d'un nouveau projet :

install.packages("BiocManager")
BiocManager::install("GenomeInfoDb")
renv::snapshot()

Le cache de renv se trouve dans ~/.local/share/renv. Tous vos packages seront installés par défaut dans ~/R/x86_64-pc-linux-gnu-library/3.5; avec renv, ça sera dans ./renv/library/R-3.5/x86_64-pc-linux-gnu/ dans le cadre de R-3.5.3, d'où l'intérêt de créer un nouveau dossier par projet.

Vous pouvez ensuite recharger votre environnement ainsi (après avoir chargé R par module) :

cd my_project
module load R/3.5.3
R

renv::restore()

Plus d'informations ici et ici.

Installer un package R

Il vous est possible d'installer un package R dans votre répertoire personnel.

D'abord, chargez la version de R appropriée :

module load R/3.4.3

Dans la plupart des cas, il sufira ensuite de lancer R, puis :

# "mon_package" à remplacer par votre package
install.packages("mon_package")

Répondez aux questions, en choisissant des mirroirs proches (France), puis acceptez d'installer le package dans le dossier proposé dans votre home.

Sinon, si nécessaire, téléchargez les sources (éventuellement sur le site cran : http://cran.r-project.org/web/packages/ ), puis créez vous un répertoire rlibrary (par exemple).

# pour un utilisateur toto
mkdir -p /home/toto/nas1/rlibrary

Puis :

# utilisateur toto
R CMD INSTALL mon_package.tar.gz -l /home/toto/nas1/rlibrary

Si le package est disponible sur le site CRAN, vous pouvez même l'installer directement à partir de R :

install.packages("mon_package",lib="/home/toto/nas1/rlibrary").

Pour la charger dans R :

library(mon_package,lib.loc="/home/toto/nas1/rlibrary")

Pour lister les packages déjà présents, depuis votre terminal :

echo 'installed.packages()' | R --slave

Info Pour installer le package rgdal, la procédure est la suivante :

qrsh
module load R-3.3.1 gdal-2.1.2 proj-4.9.3
R

install.packages('rgdal', type = "source", configure.args=c('--with-proj-include=/share/apps/bin/proj/4.9.3/include','--with-proj-lib=/share/apps/bin/proj/4.9.3/lib/'))
## on peut ensuite sortir de R et du noeud de calcul

Installation du package lme4

lme4 dépend du package R nloptr qui dépend de la librairie nlopt. nlopt est présent dans /share/apps/lib/nlopt/2.4.2; pour la charger et installer lme4 :

# cette ligne devra aussi être utilisée avant de lancer les jobs
# jobs qui devront être soumis avec qsub, avec l'option '-V'
module load gcc4.9 nlopt-2.4.2 R-3.2.0
mkdir rlibrary
cd $_
wget https://cran.r-project.org/src/contrib/nloptr_1.0.4.tar.gz
wget https://cran.r-project.org/src/contrib/lme4_1.1-14.tar.gz
cd
# l'option '-l' permet de définir un répertoire d'installation des libraries. La valeur dépend de la version de R est de l'architecture d'installation. Ces informations peuvent être retrouvées avec "R --version"
R CMD INSTALL -l ~/R/x86_64-unknown-linux-gnu-library/3.2/ --configure-vars="LIBS=-L/share/apps/lib/nlopt/2.4.2/lib CXXFLAGS=-I/share/apps/lib/nlopt/2.4.2/include" rlibrary/nloptr_1.0.4.tar.gz
R CMD INSTALL -l ~/R/x86_64-unknown-linux-gnu-library/3.2/ rlibrary/lme4_1.1-14.tar.gz

Installer un package avec dependances

Des dépendances importantes nécessiteront de passer par un conteneur singularity. En effet, nous ne pouvons pas surcharger les disques de tous les noeuds en installations nombreuses et parfois exotiques (avec parfois des problèmes de stabilité, de non compatibilité entre packages ou de dépendances croisées)...

Lectures

• Introduction à l'application des méthodes probabilistes de reconstruction phylogénétique (maximum de vraisemblance et approche Bayésienne)

• Quelques liens utiles:
• Quelques doc
• Bioinformatics with R
• Spatial Statistics with R
• Graphical User Interfaces for R
• Different graphics created with R
• Available Bundles and Packages for R

Analyser plusieurs fichiers sous R

Vous voulez lire tous les fichiers d'un dossiers pour faire des calculs dessus.

Exemple de fichiers .txt dans le dossier data de votre home directory :

> files <- list.files(« ~/data »)
> txtfiles <- files[grep(glob2rx(« *.txt »), files)]
> txtfiles
> for(fictxt in txtfiles) {
DF <- read.table(fictxt, …)
…# traitement
}

Accès à Ensembl avec R et biomaRt

On va récupérer les séquences et d'autres infos pour 2 identifiants (Ensembl) de gènes chez homo sapiens

library('biomaRt')

homo_ens_mart = useMart('ensembl', dataset = 'hsapiens_gene_ensembl')

monfiltre = 'ensembl_gene_id'
filtreval = c('ENSG00000000460','ENSG00000003989')

mesattrib = c('description', 'ensembl_gene_id', 'chromosome_name', 'start_position','end_position', 'gene_exon_intron') #unspliced gene
getBM(attributes = mesattrib, filters = monfiltre, values = filtreval, mart = homo_ens_mart)

Simulation simple forward d’une pop avec un locus bi allélique sous selction et derive

simu <-function(taillepop = 1000, max_generation = 1000, freq_init = 0.5, s0 = 0, s1 = 0)
{
male <- numeric(taillepop)
femelle <- numeric(taillepop)
freqs <- numeric(max_generation)
freqs[1] <- freq_init
for (i in 2: max_generation) {
p = freqs[i – 1]
q = 1 – p;
femelle <- sample(x = c('A', 'B'), size = taillepop, replace = T, prob = c(p, q))
male <-sample(x = c('A', 'B'), size = taillepop, replace = T, prob = c(p,q) )
descendants = paste(femelle, male, sep = '-' )
nhomoA <-sum(descendants == 'A-A')
nhomoB <-sum(descendants == 'B-B')
nhetero <-taillepop – nhomoA – nhomoB
freqs[i] <-(nhomoA * (1 + s0) + 0.5 * nhetero)/(nhomoA *
(1 + s0) + nhetero + nhomoB * (1 + s1))
}
return(freqs)
}

Essai avec une pop de taille 20 :

t = 20
plot(simu(taillepop = t, max_generation = 500), type = 'l', ylab = 'fréquence A', xlab = 'temps', main =
paste('n=', t, sep ="), ylim = c(0, 1))

Génét des pops basique et R

Construisons quelques genotypes :

alleles = c('A','B','C', 'D')
freqs = rep(0.25,4)

Tirage de 20 gamettes mâles pour un locus à partir d'une urne gamétique à 4 allèles equi frequents

m = sample(alleles, size = 20, replace = T, prob = freqs )

Idem pour les femelles :

f = sample(alleles, size = 20, replace = T, prob = freqs )

construire les genotypes des descendants :

F1 = paste(f, m, sep = ' ' )

Frequence allèliques :

gregexpr('A', F1) #rechercher l'allèle A
unlist(gregexpr('A', F1))
freqA = sum( unlist(gregexpr('A', F1)) != 1) / (2*length(F1))

Appliquer à tous les allèles à l'aide de la fonction sapply (evite le recours à une boucle for ...) :

f1_freqs = sapply(alleles , function(al) sum( unlist(gregexpr(al, F1)) != 1 ) )/ (2*length(F1))

Genotype counts :

al1 = substr(F1, 1, 1)

al2 = substr(F1, 3, 3)
obs1 < table(factor(al1), factor(al2) )
obs2 < table(factor(al2), factor(al1) )
obs = (obs1 + obs2) / 2

test du chi 2 :

test < chisq.test(obs)
k = length(alleles)
test$parameter < k*(k 1)/2 #changer le nombre de degré de libertés au lieu de (k 1) * (k 1)
test$p.value < pchisq(test$statistic, test$parameter, lower = FALSE)
names(test$statistic) < 'X squared'
names(test$parameter) < 'df'

to compute p values by Monte Carlo simulation :

test = chisq.test(obs, simulate.p.value = T)

Frequence homozygotes :

hA = sum( unlist(gregexpr('A A', F1)) != 1) / (20)
hf1_freqs = sapply(alleles , function(al) sum( unlist(gregexpr(paste(al,al,sep = ' '), F1)) != 1 ) ) / length(F1)

frequence hetero :

het_gen = sapply( 1:4, function(x) paste(alleles[x], alleles[ x], sep = ' ') ) # generer tous les genotypes héterozygotes
hetf1_freqs = sapply(het_gen , function(gen) sum( unlist(gregexpr(gen, F1)) != 1 ) ) / length(F1)
....

Phylogénie basique avec R

Aligner des seq. à l'aide d'un programme externe

library(ape)

Quelques accession numbers de Ramphocelus (Passereaux) as used in Paradis (1997)

ref <- c('U15717', 'U15718', 'U15719', 'U15720','U15721', 'U15722', 'U15723', 'U15724')
meseq <- read.GenBank(ref)

Ecriture dans un fichier au format = "interleaved", "sequential", or "fasta" :

write.dna(meseq, format = 'fasta', 'Ramphocelus.fas')
system('clustalw Ramphocelus.fas') # lancer clustal en local
t = read.tree('Ramphocelus.dnd') #lire l'arbre produit par clustalw
plot(t) #le dessiner

au fait les seq sont déjà alignées

calcul des distances génétiques par paires et phylogenie

ape

Différents modèles de substitution "JC69", "K80" (the default), "F81", "K81", "F84", "BH87", "T92", "TN93", "GG95", "logdet","paralin".

d_kim = dist.dna(meseq) #model Kimura 2 p (defaut)
print(d_kim)
njtree_kim = nj(d_kim) # Un arbre nj à partir de ces distances
write.tree(njtree_kim) #Affichage au format Newick
plot( njtree_kim, 'u' ) #On le dessine avec l'option unrooted

Changer la couleur du nom d'un taxon :

couleurs = rep('black',8)
couleurs[1] = 'red'
plot( njtree_kim, 'u', tip.color = couleurs )

Ajout d'un histogramme de données sur les 8 taxons :

mesure = rnorm(8, 0.25, sd=0.1)/10
max = max(mesure)
names(mesure) = njtree_kim$tip.label
plot(njtree_kim, x.lim=0.14, font=1, cex = 0.8)
axisPhylo()
deb= 0.08
segments(rep(deb,8),1:8,rep(deb,8) + mesure, 1:8, lwd=3, col='red')
axis(1,at=c(deb, deb+(max/2), deb+max), labels= format(c(0, (max/2) ,max ),digits=3) )
mtext('Mesure / espèces', at=deb+(max/2), side=1, line=2)

Illustrer l'incertitude sur un noeud :

nodelabels(node=14,'?',adj=1, bg='red')

Arbre de maximum de vraisemblence :

mod = DNAmodel('F84') #modèle d'evolution des sequences selon Felsenstein 1984 (tr tv et freq variables)
tree_F84 = mlphylo(model = mod, meseq, phy = njtree_kim) #pas sûr que cela marche !

Manipulation de séquences à l'aide de R

Reverse et/ou Complement Function

rev_comp < function(seq=mystr, rev=T, comp=T)
{
    if(rev==T) {
        seq < as.vector(unlist(strsplit(seq, split="")))
            seq < rev(seq)
            seq < paste(seq, collapse="")
    }
    if(comp==T) {
        seq < gsub("A", "1", seq, ignore.case = T)
            seq < gsub("T", "2", seq, ignore.case = T)
            seq < gsub("C", "3", seq, ignore.case = T)
            seq < gsub("G", "4", seq, ignore.case = T)
            seq < gsub("1", "T", seq, ignore.case = T)
            seq < gsub("2", "A", seq, ignore.case = T)
            seq < gsub("3", "G", seq, ignore.case = T)
            seq < gsub("4", "C", seq, ignore.case = T)
    }
    seq

}

exemple :

seq='tcgatcgtacgttcagcttactacgttcgttc'
rev_comp(seq, rev=T, comp=T)

ex avec seqinr :

v_seq = s2c(seq) # pour transformer en vecteur de caractères
comp(v_seq ) # complémenter
comp( rev( v_seq ) ) #reverser et complémenter

Taux de GC

ape :

GC.content(as.DNAbin(v_seq) )

seqinr :

GC( v_seq)

frequences en base

ape :

base.freq(as.DNAbin(v_seq)) # multiplier par length(v_seq) pour obtenir les comptes

seqinr :

obs = count(v_seq, 1) # 2, 3 pour les dimères trimères ...

Test d'equi fréquence de la composition en bases :

chisq.test(obs,p=c(0.25,0.25,0.25, 0.25))

Recherche de motifs

find_pattern < function(seq, motif, rev=F, comp=F){
    seq = rev_comp(seq, rev, comp);
    pos < gregexpr(motif, as.character(seq));
    #liste des positions où le motif est trouvé
    unlist(pos)
}

exemple :

find_pattern(seq='TCGATCGTACGTTCAGCT', motif='CGT', rev=F, comp=F)

Lecture de séquences à l'aide de R

On a le fichier en local et on veut le lire à l'aide du package ape :

library(ape)
read.dna('seq1.txt')

A l'aide de seqinr :

library(seqinr)
read.fasta('seq1.txt') # voir également write.fasta() read.alignment

Depuis GeneBank avec des accessions numbers :

ref< c('NM_005368') #réference du gene de la myoglobine
myoglobin< read.GenBank(ref)

Depuis une des bases de données du système ACNUC à l'aide du package seqinr :

choosebank('genbank')
query('MitCatsCDS', 'sp=felis catus AND t=cds AND o=mitochondrion') #tous les cds de la mito du chat
MitCatsCDS # Info sur l'objet de la requette
MitCatsCDS$req # liste avec détails sur les résultats de la requette
MitCatsCDS$req[[1]] # la première séquence de la liste
getSequence(MitCatsCDS$req[[1]]) # en extraire la séquence

Depuis une des bases de données structurée selon le système BioMart (biomart.org) et le package biomaRt :

library(biomaRt)
listMarts() # lister toutes les bases de données disponibes
ens_mart = useMart('ensembl') # choisir les données d'Ensembl
listDatasets(ens_mart) # lister tous les jeux de données d'Ensembl
mus_ens_mart = useDataset('mmusculus_gene_ensembl', ens_mart) # choisir le jeux de données des genes de la souris

# Définir les nom des filtres, les valeurs prises par ces filtres et les attributs à récuperer
# Filtres = les genes du chromosome 10 qui codent pour des proteines
monfiltre = c('chromosome_name', 'biotype'); # nom des filtres
filtreval = list(chromosome_name ='10', biotype='protein_coding'); #valeurs des filtres
# Attributs = récuperer les identifiants Ensembl (gene et transcrits) et MGI, la description ainsi que la position debut et position de fin sur le chr
mesattrib = c('description', 'ensembl_gene_id', 'ensembl_transcript_id', 'mgi_symbol', 'start_position','end_position');
mus_chr10_genes = getBM(attributes = mesattrib, filters = monfiltre, values = filtreval, mart = mus_ens_mart)

Requête sur les bases de données de genomes du site d'ensembl à l'aide du package RMySQL (et DBI) :

library('RMySQL')
m < dbDriver('MySQL') #specifier le type de base de données relationnelle
# Connexion en anonyme au site d'Ensembl et à la base Mus musculus core
con < dbConnect(m, host='ensembldb.ensembl.org', user='anonymous',password='', dbname='mus_musculus_core_45_36f')
dbListTables(con) #lister les tables
dbListFields(con, 'repeat_feature') # lister les champs de la table des élements repétés
dbListFields(con, 'gene') # idem pour la table des genes
#Quelques requettes sql qui retournent des data.frame
dbGetQuery(con, 'select count(*) from gene') #compter tous les genes
dbGetQuery(con, 'select * from gene LIMIT 5') #recup. des infos sur les 5 premiers genes
dbGetQuery(con, 'select * from gene where description LIKE '%meiosis%'') #tous les genes avec le mot clef "meiosis"
dbGetQuery(con, 'select * from seq_region LIMIT 10') #recup des infos sur les dix premières regions de sequences

Parallelisation avec R sur le cluster

La suite de cet article n'est valable que jusqu'à une version de R < 3.4. Pour les versions ultérieures, nous avons choisi d'utiliser des conteneurs. Vous pouvez donc vous reporter à cet article pour faire fonctionner Rmpi.

R permet l'utilisation d'un certain nombre de packages pour paralléliser du code (Parallel, Snow, snowfall, foreach, Rmpi ...).

Un certain nombre de ces packages sont installés sur le cluster, d'autres peuvent l'être sur le dossier local des utilisateurs.

Ces différents packages offrent au moins trois types de parallélisation : en local (multithread), par socket ou via mpi sur plusieurs noeuds de calcul.

La parallélisation multithread (souvent utilisée par défaut) a un intérêt limité vu le nombre faible de coeurs dans la plupart des noeuds (le code ne pourra être dispaché efficacement sur plus que 8 coeurs en général).

En revanche la parallélisation multi-noeuds par sockets ou par mpi offre plus de capacité si l'on utilise plusieurs coeurs sur plusieurs noeuds disponibles.

La contrainte majeure pour l'utilisation de ces techniques sur le cluster, est celle qui oblige à faire appel au gestionnaire de job (ici SGE) pour allouer les ressources disponibles pour un meilleur equilibrage des charges.

Il faut donc lors de la soumission d'un script R, utilisant les packages parallèles, pouvoir faire appel aux mécanismes offerts par SGE pour mieux paralléliser tout en se conformant aux règles de gesion des queues de calcul offert par le cluster.

En effet SGE fourni les moyens d'exécuter les jobs parallèles en utilisant un des environnements pré-configurés et adaptés pour la parallélisation à passage de message (mpi) ou à mémoire partagée (openmp ...).

L'option à utiliser dans ce cas est –pe suivi du nom de l'environnement parallèle (voir plus bas) et du nombre de slots désiré.

En ligne de commande:

qsub -pe robin 20

Dans un script de soumission :

#$ -pe robin 20

Quand le job est soumis avec cette option, SGE lui attribue une liste de noeuds qu'il stocke dans un fichier accessible via la variable d'environnement $PE_HOSTFILE. Il suffit donc de transmettre ce fichier à votre code R et d'utiliser cette liste dans la fonction de création du "cluster" de votre package.

Exemple de fichier de soumission test_parallel.sge:

#!/bin/bash
# nom du job:
#$ -N mon_test_sge
#
# utilise le repertoire courant pour lancer le job
#$ -cwd
#
# utiliser l'environnement parallèle robin avec 20 slots
#$ -pe robin 20

R CMD BATCH test_foreach_sock.R $PE_HOSTFILE

A soumettre avec : qsub test_parallel.sge

Votre code R doit contenir du code pour lire le fichier contenant la liste des machines qui sont allouées au job.

Exemple avec le package foreach test_foreach_sock.R :

library(doSNOW)
library(foreach)

table <- data.frame(a=rnorm(1000),b=rnorm(1000))
process <- function(table)
{for (loop in (1:nrow(table)))
{table[loop,"c"] <- with(table[loop,], a*b)
assign("table",table,envir=.GlobalEnv)
}
}

#lecture du fichier contenant la liste des machines
args <- commandArgs(TRUE)
peFile=args[1]
liste_Nodes=read.table(peFile, sep=" ",header=F, stringsAsFactors=F)

#construction de la liste des slots pour le "cluster"

node_Names=liste_Nodes[,1]
nb_slots=liste_Nodes[,2]
workers=rep(node_Names, nb_slots)

nbworkers <- length(workers)

#We will Run in parallel mode (socket) with
cl <- makeSOCKcluster(workers)
registerDoSNOW(cl)

#Pour avoir une idée de ce qui ce passe
cat(sprintf('%s backend is registered\n', if(getDoParRegistered()) 'A' else 'No'))
cat(sprintf('Running with %d worker(s)\n', getDoParWorkers()))
(name <- getDoParName())
(ver <- getDoParVersion())
if (getDoParRegistered())
cat(sprintf('Currently using %s [%s]\n', name, ver))

#pour avoir la liste des noeuds alloués
#foreach(j=1:nbworkers) %dopar% system("hostname")

#Do the job !
system.time(foreach(j=1:nbworkers ) %dopar% process(table))
stopCluster(cl)

Exemple avec le package snowfall test_snowFall_sock.R :

#Ici les packages sont installés en local sur le home dir
library(snow,lib.loc="/home/khalid/R/library")
library(snowfall,lib.loc="/home/khalid/R/library")

workerFunc <- function(n) { return(n^2) }

values <- 1:100

#lecture du fichier contenant la liste des machines
args <- commandArgs(TRUE)
peFile=args[1]

#construction de la liste des slots pour le "cluster"
liste_Nodes=read.table(peFile, sep=" ",header=F, stringsAsFactors=F)
node_Names=liste_Nodes[,1]
nb_slots=liste_Nodes[,2]

workers=rep(node_Names, nb_slots)

numWorkers <- length(workers)

#We will Run parallel mode (socket) with
sfInit( parallel=TRUE, cpus=numWorkers, type="SOCK",socketHosts=workers )
stopifnot( sfCpus() == numWorkers )
stopifnot( sfParallel() == TRUE )

cl <- sfGetCluster()

#Do the job in //
res <- parLapply(cl, values, workerFunc)

sfStop()

print(unlist(res))

Un autre test avec SNOW :

library(snow,lib.loc="/home/khalid/R/library")

workerFunc <- function(n) { return(n^2) }

values <- 1:100

#lecture du fichier contenant la liste des machines
args <- commandArgs(TRUE)
peFile=args[1]

#construction de la liste des slots pour le "cluster"
liste_Nodes=read.table(peFile, sep=" ",header=F, stringsAsFactors=F)
node_Names=liste_Nodes[,1]
nb_slots=liste_Nodes[,2]

workers=rep(node_Names, nb_slots)

#We will run in parallel mode (socket) with
cl <- makeSOCKcluster(workers)

res <- parLapply(cl, values, workerFunc)
stopCluster(cl)

print(unlist(res))

Dernier exemple avec snowfall en utilisant le mode MPI :

library(snow,lib.loc="/home/khalid/R/library")
library(snowfall,lib.loc="/home/khalid/R/library")

workerFunc <- function(n) { return(n^2) }
values <- 1:100

#lecture du fichier contenant la liste des machines
args <- commandArgs(TRUE)
peFile=args[1]

#construction de la liste des slots pour le "cluster"
liste_Nodes=read.table(peFile, sep=" ",header=F, stringsAsFactors=F)
node_Names=liste_Nodes[,1]
nb_slots=liste_Nodes[,2]
workers=rep(node_Names, nb_slots)
numWorkers <- length(workers)

#We will Run parallel mode (MPI) with
sfInit( parallel=TRUE, cpus=numWorkers, type="MPI",socketHosts=workers )
stopifnot( sfCpus() == numWorkers )
stopifnot( sfParallel() == TRUE )
cl <- sfGetCluster() #makePSOCKcluster(workers)

#Do the job
res <- parLapply(cl, values, workerFunc)

sfStop()

print(unlist(res))

Il faut lancer ce dernier exemple avec l'option -V (transmission de toutes les variables d'environnement) : sbatch --export=ALL test_snowfall_mpi.R