R Session 1 - DU Bii 2021 - module 3 "R and stats"

Practical on dataframes on the study case "Mouse fibrotic kidney"

Wednesday 3rd of March, 2021

teachers: Claire Vandiedonck & Anne Badel; helpers: Antoine Bridier-Nahmias, Bruno Toupace, Clémence Réda, Jacques van Helden

Application proposed by Olivier Sand and Jacques van Helden; adapted by Anne Badel and Claire Vandiedonck

-> In this practical, we will explore the UUO transcriptome data: "Raw" expression data were obtained after pseudocounting with Salmon algorithm (cf modules 4/5). In the paper, they further processed the data.

  • Objectives

    • reading of a first real data set
    • manipulation of the dataset
    • information retrieval
    • graphic representations

  • Data sources

Doc URL
Total RNA for the experiment on Unilateral ureter obstruction (UUO) model https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE118339
Supplementary material of the article with all the data tables (zip archive) https://zenodo.org/record/2592516

0. Avant toute chose

Faites une copie de ce notebook avec votre nom.

Définissez votre répertertoire de travail dans votre /shared/home/login/RSession1.

Tip : fonction recommandée : setwd()

1. The input data

Vous partirez des deux fichiers de données téléchargés suivants présents dans le repertoire /shared/projects/dubii2021/trainers/module3/data/pavkovic/1.datainit:

Ouvrez les d'abord dans l'édituer TSV du JupyterLab pour identifier s'ils contiennet des entêtes, des noms de lignes. Regardez également comment est remplie la colonne "color' du fichier metadata.

Ci-dessous, vous les importerez dans R et les reformaterez.

A. Le data.frame metadata

Importez-le dans R et affichez-le.

Tip : fonction recommandée : read.table() et les options header, sep et comment.char appropriées.

Tip : fonction recommandée : str() et dim.

Tip : fonction recommandée : order()

B. Le data.frame rawData

Importez-le dans R.

Tip : fonction recommandée : read.table() et les options header, sep et row.names appropriées.

Tip : fonctions recommandées : head(), tail et median.

Tip : fonctions recommandées : names(), ==.

Tip : fonctions recommandées : ls(), rm().

2. Quelques manipulations de base des dataframes

Tip : fonctions recommandées : subset() avec l'argument select, grep(), which(), identical().

Supprimez ensuite les sous-dataframes intermédiaires.

3. Generation de statistiques descriptives sur les données brutes d'expression

A. Pour chaque échantillon

fonction recommandées :

data.frame(), as.data.frame(), matrix(), cbind()

min(), max()ou range()

mean(), sd(), median() : ajoutez na.rm=T pour ne pas inclure les données manquantes

quantile() avec l'option prob=

apply() pour appliquer une fonction sur les vecteurs d’un dataframe en ligne ou en colonne

B. Pour chaque gène

fonction recommandées :

length(), unique(), ==

Pour chaque gène:

fonction recommandées :

data.frame(), as.data.frame(), matrix(), cbind()

min(), max()ou range()

mean(), sd(), median() : ajoutez na.rm=T pour ne pas inclure les données manquantes

quantile() avec l'option prob=

apply() pour appliquer une fonction sur les vecteurs d’un dataframe en ligne ou en colonne

colSums(), colMeans(), rowSums() ou rowMeans()

4. Quelques questions sur les données d'expression

A. Sur les gènes:

Tip : vous pouvez regarder la moyenne est égale à 0 dans ce cas! Ou alors vous faites la somme et comptez les sommes égales à 0.

Tip : vous pouvez utiliser le dataframe de de stats sur les gènes et la fonction which.max()

Tip : les cooronnées sur l'axe des X correspondant au jour (0,3, 7 ou 14), donc générez ces positions avec rep().

Tip : après avoir généré le plot(), utilisez la fonction lines() dans la même cellule du notebook.

Tip : après avoir généré le plot(), utilisez la fonction segments() dans la même cellule du notebook.

Est-ce toujours le même?

- extraire les 10 gènes en question

Tip : triez les données pour sélectionner les 10 plus fortes

- extraire les 10 gènes en questionreprésentation graphique (expression en fonction du jour du traitement)
   - quelle est la caractéristique de ces gènes ayant un CV très élevé ?

Dans la suite de l'analyse, nous retirerons ces gènes exprimés dans peu d'échantillons.

B. Sur les échantillons

Tip : utilisez le dataframe sample_raw_stats

Tip : utilisez la fonction colSums()

=> On voit donc que les librarires sont de taille variable. On ne peut donc pas comparer les échantillons entre eux directement.

5. Transformation des data pour les analyses ultérieures

A. transformation en variable discrète

Beaucoup d'outils d'analyse de RNASeq ont besoin de variables entières (données de comptage). Arrondissez toutes les valeurs à l'entier le plus proche avec la fonction round(). Sauvergardez ce dataframe sous le nom rawCount.

B. Distributions des comptes

AVant de transformer les données, il est nécessaire de regarder leur distribution globale et entre échantillons.

Affichez sous forme d'un seul histogramme (fonction hist()) l'ensemble des valeurs (fonction unlist()). Pensez à mettre un titre (option main =), un nom à l'axe des abcisses (option xlab =) et à l'axe des ordonnées (option ylab =).

1. distribution of raw values

2. distribution of raw counts

The distribution of raw counts is not very informative, because the range is defined by some outlier, i.e. a feature having a huge values respective to the rest of the distribution.

This problem is particularly sensitive for RNA-seq data: even with a strong zoom on the abcsissa range from 0 to 500, the histogram shows a steep drop in the first bins.

3. distribution of raw counts - right-side trimming

In order to see the relevant part of the data it is convenient to truncate the histogram to the 95th percentile.

4. distribution of raw counts - two-sides trimming

For visualization, it is sometimes convenient to truncate both tails of the distribution, in order to avoid both the outliers (right tail) and the sometimes large number of zero values. This is particularly useful for single-cell RNA-seq data, which can have a "zero-inflated" distribution, but it is also useful for other data types.

5. log2-transformed values of raw counts

A typical normalising approach is apply a log2 transformation .

For RNA-seq data, this might however create a problem when the counts of a given gene in a given sample is 0. To circumvent this, we can add an epsilon parameter (epsilon = 0.1) before the log2 transformation (function log2()).

We can now inspect the distribution of counts per sample with the boxplot() function.

- boxplot of the raw counts
- boxplot of the log2-transformed values of raw counts

C. feature filtering

We filter out the features having very low values across all the samples, which we consider as undetected.

Remember to share your graphical window using the command (par(mfrow=c(1,2)) and then return to the original parameters using the command (par(mfrow=c(1,1)).

We use the remaining genes for the subsequent analyses.

6. Représentation graphique des 10 gènes ayant les plus grandes valeurs de CV (optionnel)

Que se passe t'il après ces opérations de filtrage ?

Qu'en pensez-vous?

Pensez à sauvegardez vos données

à l'aide de la fonction save() ou de la fonction write.table()