Obtenir une matrice de valeurs d'expression avec easyRNASeq

Introduction

Plus à venir…

Protocole

Note: le protocole qui suit fonctionne avec les dernières versions de R et Bioconductor (3.4 et 3.5 respectivement, en date de juin 2017).

Il vous faut les fichiers suivants:
- Un répertoire contenant TOUS vos fichiers d'alignement, filtrés. Note: Si vous n'avez pas vos fichiers tous au même endroit (très probables car ils seraient répertorié selon type, condition,temps,etc.), il est possible de créer des liens symboliques vers les fichiers sources qui eux seront au même endroit.
- Le fichier d'annotations vous ayant servi à faire les alignements, en format GTF ou GFF3.
Démarrez R et chargez la librairie easyRNASeq:

> library("easyRNASeq")

Il faut se placer dans le répertoire où se trouve les liens symboliques:

> setwd("/chemin/vers/les/liens/symboliques")

On utilise les méthodes de manipulation de fichiers de R ainsi que la méthode getBamFileList de easyRNASeq:

# Permet de créer un vecteur contenant la liste des fichiers désirés
# Si nécessaire, utilisez un patron commun pour filtrer les fichiers désirés
> bamFiles<-list.files(
    path=".",
    pattern="*.properly_paired.*.bam$"
  )
# Transformer le vector pour extraire les vrais chemins des fichiers
> bamFiles<-Sys.readlink(bamFiles)
# Créer l'objet BamFilesList nécessaire
> bamFiles<-getBamFileList(bamFiles)

Il faut ensuite traiter le fichier d'annotation pour le rendre utilisable par easyRNASeq:

# Si on utilise un fichier GTF
# Pour GFF3 simplement remplacé gtf par gff3
> annPar<-AnnotParam(
    datasource="/un/chemin/vers/fichier/GTF",
    type="gtf"
  )
 
# easyRNASeq doit par la suite créer un table d'annotations
# via l'analyse des exons pour créer une carte virtuelle des
# transcrits
> annPar<-createSyntheticTranscripts(annPar)
 
# Avec le fichier genes.gtf de Illumina-iGenome pour GRCh38 v.2015 
# Ça donne ça:
> getAnnotation(annPar)
No validation performed at that stage
GRanges object with 276766 ranges and 3 metadata columns:
           seqnames               ranges strand |           exon transcript
              <Rle>            <IRanges>  <Rle> |          <Rle>      <Rle>
       [1]     chr1       [11874, 12227]      + | DDX11L1.exon.1  DDX11L1.0
       [2]     chr1       [12613, 12721]      + | DDX11L1.exon.2  DDX11L1.0
       [3]     chr1       [13221, 14409]      + | DDX11L1.exon.3  DDX11L1.0
       [4]     chr1       [14362, 14829]      - | WASH7P.exon.11   WASH7P.0
       [5]     chr1       [14970, 15038]      - | WASH7P.exon.10   WASH7P.0
       ...      ...                  ...    ... .            ...        ...
  [276762]     chrY [57190940, 57191085]      + |  IL9R-2.exon.6   IL9R-2.0
  [276763]     chrY [57191798, 57191999]      + |  IL9R-2.exon.7   IL9R-2.0
  [276764]     chrY [57192603, 57192708]      + |  IL9R-2.exon.8   IL9R-2.0
  [276765]     chrY [57194043, 57194127]      + |  IL9R-2.exon.9   IL9R-2.0
  [276766]     chrY [57196336, 57197337]      + | IL9R-2.exon.10   IL9R-2.0
              gene
             <Rle>
       [1] DDX11L1
       [2] DDX11L1
       [3] DDX11L1
       [4]  WASH7P
       [5]  WASH7P
       ...     ...
  [276762]  IL9R-2
  [276763]  IL9R-2
  [276764]  IL9R-2
  [276765]  IL9R-2
  [276766]  IL9R-2
  -------
  seqinfo: 282 sequences from an unspecified genome; no seqlengths

Il reste à créer un objet RnaSeqParam pour contenir les paramètres du décompte:

# Selon les paramètres du rnaseq effectué
> bamPar<- BamParam(
    paired = T,
    stranded = F
  )
 
# Selon les informations souhaitées: 
# countBy = c("exons", "features", "genes","transcripts")
# precision = c("read", "bp")
# Dans la plupart des cas, c'est "read" qui sera utilisé
> rnaSeqPar<-RnaSeqParam(
    annotParam = annPar,
    bamParam = bamPar,
    countBy = "genes",
    precision = "read"
  )

Dernière étape: faire le décompte Sur un cluster, on peut utiliser le paramètre nnodes pour spécifier le nombre de CPU à utiliser pour accélérer le calcul; la valeur par défaut est nnodes=1.

> exp <- simpleRNASeq(
    bamFiles=bamFiles,
    param=rnaSeqPar,
    verbose=TRUE
  )

Ça va prendre du temps…
Plus à venir…