Module

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.

Installer un package R

Il vous est possible d’installer un package R dans votre répertoire personnel. Pour cela, 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), et enfin:

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:

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

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

Lectures

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 = ' ' )

freq 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)

freq 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)

freq 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 neoud

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 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)

Requette 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

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) : qsub -V test_snowfall_mpi.R