<figure><img src="image_1.png" /></figure>

HAL
open science

# Caractérisation de gènes et d'ARN longs non codants régulés par BMP15 chez la truie et la vache.

Gabryelle Agoutin

► To cite this version:

Gabryelle Agoutin. Caractérisation de gènes et d'ARN longs non codants régulés par BMP15 chez la truie et la vache.. Sciences du Vivant [q-bio]. 2022. hal-04225079

HAL Id: hal-04225079

https://hal.inrae.fr/hal-04225079

Submitted on 2 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Caractérisation de gènes et d'ARN longs non codants régulés par BMP15 chez la truie et la vache.

Gabryelle AGOUTIN

M2 Bioinformatique et Biologie des systèmes

2021-2022

## Encadrantes :

- Carine Genet (GenROC, Génomique des Ruminants, Ovins et Caprins)

- Sarah Maman (SIGENAE, Système d'Information et Génétique des Animaux d'Elevage)
Résumé

Pour l'amélioration des systèmes de production animale, des études se sont centrées sur les fonctions de reproduction. Parmi ces fonctions se trouvent la folliculogenèse, processus de développement du follicule. Plusieurs gènes appartenant à la superfamille du TGF-β sont impliqués dans le contrôle de celle-ci. Un de ces gènes est celui codant le facteur BMP15 (bone morphogenetic protein 15), sécrété de façon paracrine par l'ovocyte vers les cellules de la granulosa. Il a été montré que ce facteur a des effets sur la prolificité chez la brebis et qu'il existe un déterminisme génétique. L'équipe s'intéresse à la différence d'effet du BMP15 sur l'expression des gènes entre une espèce poly-ovulante (la truie) et une espèce mono-ovulante (la vache). En plus de caractériser les gènes différentiellement exprimés, nous nous intéressons aussi aux lncRNAs (long non-coding RNA) en tant qu'acteurs potentiels de la régulation.

Dans le cadre d'un projet ANR, une analyse transcriptomique des cellules de granulosa issues de 10 vaches et 10 truies, traitées ou non avec du BMP15 pendant 48h a été réalisée. L'objectif de mon stage était de reprendre l'analyse de ces données en partant des séquences jusqu'à la caractérisation des gènes différentiellement exprimés.

Pour cela, j'ai utilisé un pipeline récent, le nf-core/rnaseq. Il est lancé avec Nextflow ce qui le rend reproductible, facile à utiliser et accessible. Ce workflow permet de générer la matrice de comptage et de reconstruire les transcrits servant pour la prédiction des lncRNAs. C'est FEELnc qui a été utilisé pour cette prédiction. La matrice de comptage annotée avec les lncRNAs a servi pour l'analyse différentielle ce qui nous a permis d'obtenir une liste de gènes différentiellement exprimés (DEGs). J'ai réalisé ensuite une analyse de sur-représentation des termes GO pour caractériser les fonctions moléculaires et les processus biologiques impliqués dans la réponse à BMP15 chez nos deux espèces d'études.

Concernant l'obtention d'un répertoire de lncRNAs dans les cellules de granulosa, nous avons identifié 2285 nouveaux lncRNAs chez la vache et 2817 chez la truie. L'analyse différentielle a permis d'identifier 436 gènes différentiellement exprimés chez la vache et 1008 chez la truie (padj<0.05 et |LFC|>1). L'analyse de sur-représentation des termes GO à montrer des différences dans l'action du BMP15 entre la truie et la vache notamment au niveau des processus liés au cycle cellulaire.
Abstract

For the improvement of animal production systems, studies have focused on reproductive functions. Among these functions, folliculogenesis is the process of follicle development. Several genes, belonging to the TGF-β superfamily, are involved in the control of folliculogenesis. One of these genes, the BMP15 factor (bone morphogenetic protein 15), is secreted in a paracrine way by the oocyte towards the granulosa cells. It has been shown that BMP15 modulates prolificity in ewes and that these effects are mediated by a genetic determinism. The team is interested in the difference in effect of BMP15 on gene expression between a poly-ovulating species (the sow) and a mono-ovulating species (the cow). In addition to characterizing differentially expressed genes, we are also interested in lncRNAs (long-non-coding RNAs) as potential regulatory players.

In the framework of an ANR project, a transcriptomic analysis of granulosa cells from 10 cows and 10 sows, treated or not with BMP15 for 48h, was undertaken. The objective of my internship was to start over the analysis of these data from the sequences to the characterisation of differentially expressed genes.

For that, I launched a recent pipeline, the nf-core/rnaseq, launched with Nextflow which makes it reproducible, easy to use and to access. This workflow allows to generate the counting matrix and to reconstruct the transcripts used for the prediction of lncRNAs. FEELnc was used for this prediction. The count matrix annotated with the lncRNAs was used for differential analysis which allowed us to extract a list of differential expressed genes (DEGs). I then performed an over-representation analysis of GO terms to characterize the molecular functions and biological processes involved in the response to BMP15 in cow and sow.

The differential analysis revealed a total of 436 differentially expressed genes in cows and 1008 in sow. Prediction of lncRNAs gave us 2285 novel lncRNAs in cows and 2817 in sows. The over-representation analysis of GO terms showed differences in the action of BMP15 between sow and cow, especially in processes related to cell cycle.
Remerciements

Je remercie en premier lieu mes encadrantes, Carine Genêt et Sarah Maman, notamment pour leur pédagogie et leur présence. Elles se sont montrées patientes et à l'écoute à chaque instant de ce stage. Apprendre à leurs côtés a été une chance incroyable et je les remercie pour cette superbe opportunité.

Je remercie chaleureusement Stéphane Fabre, Laurence Drouilhet, Gwenola Tosser-Klopp, Florence Plisson-Petit, Julien Sarry, Florent Woloszyn et Maxime Ben Braiek de l'équipe GenROC, qui m'ont accueilli et surtout bien intégré, ce qui a rendu ce stage encore plus agréable.

Je remercie aussi l'équipe Sigenae pour m'avoir offert la chance de participer à l'une de leur réunion pour présenter mon stage. Cela m'a surtout permis d'avoir un petit aperçu du fonctionnement d'une plateforme bioinformatique, ce qui est pour moi une immense chance.

Enfin, je remercie Alexia Rivero, étudiante de la promotion 2019-2021 du Master Bioinformatique et Biologie des systèmes, d'avoir toujours pris le temps de me donner de précieux conseils durant ces deux années.

We are grateful to the genotoul bioinformatics platform Toulouse Occitanie (Bioinfo Genotoul, https://doi.org/10.15454/1.5572369328961167E12) for providing help and/or computing and/or storage resources.
Table des matières

<table><tr><td>Résumé</td><td>1</td></tr><tr><td>Abstract</td><td>2</td></tr><tr><td>Remerciements</td><td>3</td></tr><tr><td>1. Introduction</td><td>6</td></tr><tr><td>2. Matériels et méthodes</td><td>10</td></tr><tr><td>2.1 Obtention des échantillons biologiques</td><td>10</td></tr><tr><td>2.2 Production des séquences et analyse de leurs qualités</td><td>10</td></tr><tr><td>2.3 Génération de la matrice de comptage</td><td>12</td></tr><tr><td>2.4 Reconstruction et quantification des transcrits "de novo"</td><td>12</td></tr><tr><td>2.5 Prédiction des lncRNAs</td><td>13</td></tr><tr><td>2.6 Analyse statistique des résultats de FEELnc</td><td>15</td></tr><tr><td>2.7 Valorisation de la matrice de comptage</td><td>15</td></tr><tr><td>2.7.1 Annotation des gènes : obtention des noms de gènes</td><td>15</td></tr><tr><td>2.7.2 Annotation des lncRNAs dans la matrice</td><td>16</td></tr><tr><td>2.8 Analyse différentielle</td><td>17</td></tr><tr><td>2.9 Vérification d'un MSTRG : IGV</td><td>18</td></tr><tr><td>2.10 Analyse de la Gene Ontology</td><td>18</td></tr><tr><td>3. Résultats</td><td>22</td></tr><tr><td>3.1 Contrôle qualité</td><td>22</td></tr><tr><td>3.2 Matrice de comptage</td><td>22</td></tr><tr><td>3.4 FEELnc : Analyse statistique des résultats</td><td>24</td></tr><tr><td>3.4.1 Résultat de la prédiction des lncRNAs</td><td>24</td></tr><tr><td>3.4.2 Proportion de lncRNAs Ensembl vs ceux prédits par FEELnc</td><td>24</td></tr><tr><td>3.5 Analyse différentielle</td><td>26</td></tr><tr><td>3.6 Exemple de vérification de transcrits "de novo" d'intérêts</td><td>28</td></tr><tr><td>3.7 Analyse de la Gene Ontology</td><td>30</td></tr><tr><td>4. Discussion et perspectives</td><td>31</td></tr><tr><td>4.1 Impact sur le stage du délai de récupération des FASTQ au NCBI</td><td>31</td></tr><tr><td>4.2 Difficulté pour la quantification des transcrits "de novo"</td><td>31</td></tr><tr><td>4.3 Confiance apportée aux prédictions de FEELnc</td><td>33</td></tr><tr><td>4.4 Réseaux de gènes</td><td>33</td></tr><tr><td>Références</td><td>35</td></tr><tr><td>Liste des abréviations</td><td>39</td></tr></table>
<figure><img src="image_2.png" /><figcaption>Figure 1 : Schéma des différentes étapes de la folliculogénèse au niveau de l'ovaire.</figcaption></figure>
# 1. Introduction

Dans le cadre du défi de l'agriculture et de l'alimentation de demain, un des objectifs de la recherche agronomique est l'amélioration des techniques d'élevage. Entre la pression démographique et celle de l'environnement, il est devenu primordial d'optimiser la rentabilité et la productivité dans l'industrie des ruminants tout en limitant les émissions de gaz à effet de serre. Un élément central pour tout système de production animale est la maîtrise des fonctions de reproduction. Elle contribue à la diffusion du progrès génétique et à une alimentation plus durable en contrôlant, par exemple, le nombre de descendants. Une des composantes majeures de l'efficacité de la reproduction est la prolificité (Juengel et al, 2021), c'est-à-dire le nombre de descendants obtenus par une mère en un an.

La reproduction sexuée débute au moment de la fécondation, lorsqu'un gamète mâle (spermatozoïde) fusionne avec un gamète femelle (ovule). L'ovogenèse, c'est à dire la croissance et l'évolution du gamète femelle, s'effectue au niveau de l'ovaire à l'intérieur des follicules ovariens au cours de leurs processus de développement appelé folliculogenèse (cf. Figure 1).

La folliculogenèse est un ensemble de processus continus par lesquels un follicule primordial évolue et accompagne l'ovocyte jusqu'à l'ovulation ou dans la majorité des cas, régresse par atrésie. Les follicules primordiaux constituent un pool de “réserve” élaboré principalement pendant la vie intra-utérine. A partir d'une sélection de follicules primordiaux composés d'une couche unique de cellules aplaties autour de l'ovocyte, une phase de croissance débute et génère un follicule primaire. Les cellules aplaties deviennent cubiques et prennent le nom de cellules de la granulosa. Ces follicules primaires continuent de croître pour passer au stade folliculaire secondaire dans lequel les cellules de la granulosa forment plusieurs couches autour de l'ovocyte, et les cellules de la thèque se forment. L'apparition d'une cavité folliculaire (ou antrum) dans la granulosa caractérise le passage du follicule secondaire au follicule antral. Enfin, le dernier stade avant l'ovulation est celui du follicule de Graaf (follicule mûr ou pré-ovulatoire) qui va éclater pour libérer l'ovocyte mûr (ovule). Le follicule, en absence de fécondation deviendra un corps jaune qui entrera dans un processus de dégénérescence. Ainsi, le follicule fournit à l'ovocyte l'environnement adéquat pour son développement pendant lequel le follicule est sous l'action d'hormones gonadotropes (FSH et LH), ou d'interactions paracrines entre l'ovocyte et les cellules folliculaires (granulosa et thèque). Plusieurs gènes appartenant à la superfamille du TGF-β sont impliqués dans le contrôle de la folliculogenèse et du taux d'ovulation (Fabre et al, 2006).
<figure><img src="image_3.png" /></figure>

**Figure 2 :** Facteurs de croissance spécifiques à l'ovocyte dans la folliculogenèse (d'après Jozkowiak et al, 2020).
Ainsi, l'importance de certains membres de la famille des protéines
morphogénétiques osseuses (BMPs) et des facteurs de croissance et de différenciation
(GDFs) ou d'un de leurs récepteurs (BMPR1B) ont été mis en évidence. Chez la brebis, il
existe un déterminisme génétique de la variabilité de la prolificité résultant de mutations dans
BMPR1B (mutations dite "FecB"; Mulsant et al, 2001) dans GDF9 (mutations dite "FecG";
Hanrahan et al, 2004) ou dans BMP15 (mutations dite "FecX"; Galloway et al, 2000; Demars
et al, 2013). Selon si la mutation est hétérozygote ou homozygote, des phénotypes allant de
l'hyperprolificité jusqu'à la stérilité sont observés (Galloway et al, 2000; Hanrahan et al,
2004; Silva et al, 2011).

BPM15 est un facteur codé par le gène BMP15 situé sur le chromosome X chez tous
les mammifères. Il appartient à la superfamille des TGF-β (transforming growth
factor) (Moore *et al*, 2003). Il est sécrété de façon paracrine par l'ovocyte vers les cellules de
la granulosa. Son expression par l'ovocyte induit l'expression de gènes spécifiques dans les
cellules de la granulosa (Feuerstein *et al*, 2006). Cette protéine s'homodimérise avec
elle-même ou s'hétérodimérise avec GDF9 puis se fixe à son récepteur BMPR2 (Bone
morphogenetic protein receptor 2), et active ensuite la voie de signalisation des SMAD
(Fabre et al. 2006)(cf. Figure 2). Il a alors plusieurs actions sur la folliculogenèse dont la
régulation de la prolifération et de la différenciation des cellules de granulosa (Otsuka *et al*,
2000; Moore *et al*, 2003) et l'inhibition des effets de la FSH en diminuant le niveau
d'expression de l'ARNm du récepteur à la FSH (Otsuka *et al*, 2001).

Dans le contexte de l'étude du déterminisme génétique de la prolificité, l'équipe
GenROC s'intéresse non seulement aux gènes différentiellement exprimés induits par
l'expression de BMP15, mais aussi aux ARN longs non codants (lncRNAs). L'intérêt pour le
"non-codant" s'est fortement intensifié ces dernières années notamment depuis la découverte
de leurs implications dans différentes maladies et cancers (Bhan et al, 2017). Les lncRNAs
appartiennent à la plus grande classe de gènes non codants (Gudenas et al, 2019). Ce sont des
transcrits d'une longueur supérieure à 200 nucléotides sans capacité de codage des protéines
(Ponting et al, 2009). Ils sont épissés de manière prédominante et sont en général faiblement
exprimés avec une expression tissu-spécifique. La séquence de leurs régions exoniques est
faiblement conservée entre les espèces mais les régions promotrices des lncRNAs sont
conservées (Gudenas et al, 2019). Leurs fonctions biologiques sont variées et sont classées
dans quatre thèmes mécanistiques : ils agissent comme un signal (d'initiation de la transcription), un leurre (empêcher la transcription), un guide (liaison des facteurs de transcription, recrutement de protéines pour l'initiation à la transcription) ou comme un support (médiateur d'interaction physique entre protéines et ARN non codant)(Gudenas et al, 2019).

Les objectifs de mon stage s'inscrivent dans ce contexte et portent en particulier sur l'identification des gènes différentiellement régulés par BMP15 dans les cellules de la granulosa de deux espèces divergentes par leur nombre de descendants (la vache, une espèce mono-ovulante et la truie une espèce poly-ovulante). J'ai également travaillé sur l'obtention d'un inventaire des lncRNAs exprimés dans les cellules de granulosa porcine et bovine répondant ou non à l'action de BMP15. Enfin, j'ai réalisé les analyses d'ontologies des gènes différentiellement exprimés.

A mon arrivée au laboratoire, j'ai repris l'analyse des données de RNAseq générées dans le cadre du projet ANR (Agence Nationale de la Recherche) Monopoly (mono ou poly ovulation, rôle évolutif du gène BMP15 chez les mammifères, ANR 2010, porteur S. Fabre). Les échantillons biologiques et leur séquençage avaient déjà été analysés par ma co-encadrante Sarah Maman en 2020 avec des outils différents de ceux utilisés pour mon stage. Les FASTQ ont été publiés sur le NCBI en mars 2021. Leur analyse s'étend du traitement des lectures de séquences de RNAseq à leur interprétation biologique.
## 2. Matériels et méthodes

### 2.1 Obtention des échantillons biologiques

Les ovaires de 10 vaches et 10 truies ont été prélevées à l'abattoir sur des animaux tout-venant, jeunes et impubères. Après dissection, les cellules de granulosa ont été mises en culture pendant 48h *in vitro* avec (traitées) ou sans (contrôles) la protéine BMP15 humaine (R&D systems) à une concentration de 200 ng/µl. L'ARN total a été extrait avec le kit Nucleospin RNA II de Macherey-Nagel selon le protocole fourni. La qualité des ARNs et leurs concentrations ont été déterminés par dosage sur nanodrop ND-8000.

### 2.2 Production des séquences et analyse de leurs qualités

Les 20 ARNs par espèce (soit 10 traités et 10 contrôles) ont été purifiés, rétro-transcrits et séquencés en paire ("paired-end") sur un séquenceur Illumina HiSeq2000 par la plateforme GetPlaGe (https://get.genotoul.fr/). Chaque échantillon a été déposé sur 2 lanes différentes de séquenceur car lorsque le séquençage a été réalisé en 2013, se posait la question d'un possible biais lié à l'effet "lane" du séquenceur. Après vérification d'une forte corrélation (r>0.9) au sein d'un même échantillon biologique, les données de séquences d'un même individu ont été additionnées.

Les FASTQ (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA713946) ont été téléchargés à l'aide de l'outil fastq-dump contenu dans **sratookit** (version 3.0.0) (https://github.com/ncbi/sra-tools) à partir de leurs identifiants SRR.

Le contrôle de qualité des FASTQ a été réalisé avec **FastQC** (version 0.11.9) au sein du pipeline Nextflow nf-core/rnaseq, qui fournit un rapport permettant de vérifier la qualité des séquences ainsi que le bon déroulement des outils bioinformatiques. Certaines informations dans ce rapport sont importantes à prendre en compte comme la taille des reads ou encore le pourcentage du contenu en base GC. Cela permet aussi de vérifier que les adaptateurs de séquençage Illumina ont bien été enlevés.
<figure><img src="image_4.png" /><figcaption>**Figure 3 :** Schéma du pipeline nf-core rnaseq au format "ligne de métro". Il représente les différents outils utilisés une fois la ligne de commande lancée (from : https://github.com/nf-core/rnaseq).</figcaption></figure>
## 2.3 Génération de la matrice de comptage

La matrice de comptage est générée à l'aide du pipeline **nf-core/rnaseq** (version 3.5; Harshil Patel *et al*, 2021) (cf. Figure 3) qui utilise **Nextflow** (version 21.10.6). **Nextflow** est codé en DSL (domain-specific language) et permet d'exécuter une suite de tâches (workflow) sur n'importe quelle infrastructure. Grâce à des conteneurs logiciels, un pipeline créé avec **Nextflow** sera reproductible. Cela confère à ce pipeline un caractère FAIR (Findable, accessible, interoperable and reusable) très attendu dans les publications bioinformatiques et les projets ANR. Les données doivent être faciles d'accès, réutilisables et un pipeline lancé plusieurs fois doit toujours générer le même résultat.

L'étape d'alignement des lectures (ou reads) sur le génome de référence est réalisée avec l'outil **STAR** (version 2.7.6a; Dobin *et al*, 2013) qui est l'un des aligneurs les plus utilisés actuellement. Les génomes de références sont récupérés depuis Ensembl, un navigateur de génomes (www.ensembl.org, v101, August 2020). Les génomes utilisés sont ceux de Bos Taurus (ARS-UCD1.2) pour l'espèce bovine et Sus scrofa (Sscrofa11.1) pour l'espèce porcine. Une fois l'alignement réalisé, la deuxième étape est la compression des fichiers SAM (sequence alignment map) en BAM (Binary alignment map) avec la suite de programme **Samtools** (version 1.14; Li *et al*, 2009), qui permet d'interagir avec les données de séquençage à haut débit. Ces BAM sont triés avec **Samtools**. Le logiciel **RSEM** (version 1.3.1; Li & Dewey, 2011) permet de quantifier les niveaux d'expression des transcrits et des gènes à partir de données RNA-Seq.

Le pipeline **nf-core rnaseq** produit un rapport html avec **MultiQC** (version 1.12). Il contient l'interprétation des différents logs des outils utilisés dans le workflow. **MultiQC** produit, entre autres analyses, le contrôle qualité des séquences avec **FastQC** et l'analyse de l'alignement des lectures sur génome de référence avec la représentation graphique des logs de **STAR**.

## 2.4 Reconstruction et quantification des transcrits “de novo”

Une étape importante dans le pipeline **nf-core/rnaseq** pour notre analyse est celle de la reconstruction des transcrits. Pour cela, le pipeline utilise **StringTie** (version 2.1.7; Pertea *et al*, 2015), un assembleur rapide et très efficace d'alignements RNA-Seq en transcrits
potentiels. **StringTie** est lancé avec l'option --stringtie ignore_gtf qui permet un assemblage
*de novo* sans se limiter au fichier GTF (General feature format) de référence. Cela permet
d'obtenir de nouveaux transcrits potentiellement inconnus dans l'annotation (nommé
MSTRG par StringTie). Pour chaque échantillon, un fichier d'annotation au format GTF est
généré. Ce fichier contient des informations telles que les identifiants, les biotypes ou encore
les positions des transcrits sur le génome de référence. **StringTie Merge** permet ensuite de
fusionner les différents GTF en un seul fichier GTF mergé.

Cependant, ces transcrits nouvellement identifiés ne sont pas intégrés à la matrice de
comptage par le pipeline nf-core. Afin de disposer de leurs comptages bruts, il est nécessaire
de relancer une quantification. Nous avons choisi de garder **StringTie** pour cette
quantification car la reconstruction des transcrits est faite avec cet outil. Nous avons mis
comme annotation de référence la concaténation entre le GTF de référence Ensembl récupéré
précédemment et le fichier GTF des transcrits reconstruits par **StringTie**. Ce fichier GTF
"mergé" contient bien les annotations des MSTRG, ceux-ci ont donc pu être quantifiés. Nous
disposons alors de leurs valeurs TPM (transcripts *per* million). Le package R
**IsoformSwitchAnalyzeR** (version 1.16.0; Vitting-Seerup, 2017) qui fait aussi partie de la
suite **StringTie**, permet de convertir les comptages TPM en comptage bruts et de filtrer les
transcrits reconstruits. Il supprime les MSTRG non exprimés et ceux qui chevauchent un
exon. Il corrige aussi les problèmes d'annotation de **StringTie** car, parfois, **StringTie** appelle
un transcrit MSTRG alors qu'il dispose déjà d'un Ensembl ID de référence.

2.5 Prédiction des lncRNAs

Pour la prédiction des lncRNAs, nous utilisons **FEELnc** (FlExible Extraction of LncRNAs; version 0.2; Wucher *et al*, 2017). **FEELnc** est un programme d'annotation sans alignement qui prédit les lncRNAs en se basant sur un modèle de forêts aléatoires à partir de transcrits assemblés. Il utilise une "Random Forest" ou Forêt d'arbres de décision pour déterminer si un transcrit est un lncRNA ou mRNA en fonction de différents calculs. C'est une technique de "Machine learning" qui peut être utilisée pour des tâches de classification. Une donnée passe un ensemble de tests dont la réponse est oui/non jusqu'à la décision finale en bas de l'arbre. Chaque arbre est entraîné sur un sous-ensemble de notre jeu de données (jeu d'entraînement) puis les résultats de tous les arbres de décision sont combinés.
La prédiction des lncRNAs se déroule en trois étapes : une étape de filtrage, une étape de calcul du potentiel codant des transcrits et une étape de classification des lncRNAs nouvellement annotés.

La première étape *FEELnc Filter* permet la suppression des transcrits qui sont prédits comme codants, comme ceux qui chevauchent un exon dans le même sens et les transcrits inférieurs à 200 nucléotides. **FEELnc** permet de garder ou non les modèles intergéniques mono-exon. Nous avons dû choisir entre ces deux possibilités et pour cela j’ai testé les deux options.

La seconde étape *FEELnc codpot* permet le calcul du potentiel codant (CPS: coding potential score). Ce potentiel codant est déterminé à l'aide de trois "scores prédicteurs" qui sont incorporés aux forêts aléatoires. Le premier score annote tous les cadres ouverts de lectures ORF (pour open reading frame). Il y a des ORF de type 0 (strict) jusqu'au type 4 en fonction de la couverture de l'ORF sur la taille du transcrit. Le second score produit les fréquences multi-Kmer entre les ARNm et les lncRNAs. Un compteur de Kmer (KIS) spécifique au mRNA ou au lncRNA permet de trouver le total de Kmer de taille k sur le transcrit. Un "score de k-mer" est alors associé à chaque séquence. Le troisième score concerne la longueur totale de la séquence ARN.

Le seuil optimal de CPS permet de maximiser la spécificité et la sensibilité. **FEELnc** calcule 2 seuils : un pour les lncRNAs et un pour les mRNAs. Ils seront de "haute confiance" en dépassant leur seuil respectif. Il permet aussi de prendre en compte les transcrits proches du seuil (TUCp transcrit of unknown coding potential ou transcrit ambigu).

Dans **FEELnc**, il y a deux stratégies principales pour le jeu d'entraînement : la méthode shuffle et la méthode Intergenic. Cela est dû au fait que les lncRNAs ne sont généralement pas disponibles pour tous les organismes contrairement aux mRNA. Or, pour qu'un algorithme d'apprentissage automatique soit efficace, il faut un grand nombre de lncRNAs et mRNAs disponibles dans le jeu d'entraînement. C'est pour modéliser les lncRNAs que ces deux méthodes sont lancées. La méthode intergenic est une méthode naïve qui extrait des séquences intergéniques aléatoires de longueur L en tant qu'ARN non-codant. La méthode shuffle est basée sur l'idée que les lncRNAs sont dérivés de "débris" de gènes codants. Elle mélange au hasard les séquences d'ARNm en préservant une fréquence k-mer donnée.
La troisième étape *FEELnc classifier* est celle permettant de cataloguer les lncRNAs nouvellement annotés en fonction de leur gène codant le plus proche. **FEELnc** utilise une stratégie de fenêtre glissante pour parcourir le lncRNA et signaler les transcrits à proximité de cette fenêtre. Les règles font intervenir le sens et le type de chaque interaction.

Ce qu'il est important de noter c'est qu'à ce moment de l'analyse nous disposons de deux sources de lncRNAs: ceux déjà connus comme tels dans le GTF de référence Ensembl et ceux prédits par **FEELnc**. Pour ces derniers, ils sont soit des lncRNAs "de novo" trouvés à partir des nouveaux transcrits reconstruits par **StringTie** soit des lncRNAs qui sont déjà annotés dans le GTF de référence mais sans biotype "lncRNA" (par exemple un biotype "pseudogène").

## 2.6 Analyse statistique des résultats de FEELnc

L'objectif de cette analyse est d'explorer plus en détail la sortie du logiciel **FEELnc** afin de caractériser nos lncRNAs. Pour visualiser la proportion de lncRNAs connus dans Ensembl face à ceux prédits par **FEELnc**, nous avons utilisé un diagramme de Venn généré avec l'outil **jVenn** (Bardou *et al*, 2014). Pour cela, il faut récupérer les identifiants correspondants à des biotypes lncRNAs dans le fichier Ensembl de référence et les identifiants dans le fichier candidat des lncRNAs de **FEELnc** pour ceux qui proviennent d'Ensembl et qui n'avaient pas le biotype "lncRNA" avant.

## 2.7 Valorisation de la matrice de comptage

### 2.7.1 Annotation des gènes : obtention des noms de gènes

Notre matrice de comptage est indexée par des identifiants Ensembl ou Ensembl ID. Pour plus de visibilité, j'ai codé un script avec **R** (version 4.3.1) et en particulier le package **biomaRt** (version 2.50.3 ; Huber, 2017) pour l'annotation des gènes. Il permet de retrouver les noms des gènes à partir des identifiants Ensembl en fouillant dans plusieurs bases de données. Ce qu'on appelle ici nom de gène correspond au symbole attribué par le HGNC (HUGO Gene Nomenclature Committee), qui est le comité qui approuve un nom unique et
valable pour chaque gène humain connu. Nous avons choisi de confronter nos identifiants de *Bos taurus* et *Sus scrofa* aux identifiants *Homo sapiens* car le génome humain est nettement plus annoté.

Le package **biomaRt** permet l'extraction de données de la DataBase Biomart (i.e Ensembl). Il fait le lien entre les identifiants Ensembl de nos espèces d'intérêt et le symbole HGNC du gène dans plusieurs bases de données pour annoter le plus de gènes possible.

Nous avons cherché les symboles HGNC dans 3 bases de données différentes : NCBI (https://www.ncbi.nlm.nih.gov/), Wikigene (https://www.wikigenes.org/) et Uniprot (https://www.uniprot.org/). Le script consiste donc à créer une colonne "Gene name" et à y insérer la valeur correspondant à l'Ensembl ID. Nous commençons par chercher dans External gene name (= identifiant NCBI), nous ajoutons les symboles et s'il n'y a pas de valeur correspondante à notre ID alors nous mettons la valeur "NA" (Not available ie valeur manquante). Nous récupérons les ID de ces valeurs manquantes et nous les cherchons dans la base de données Wikigene. Enfin, nous prenons encore les derniers NA non annotés par Wikigene et nous cherchons dans Uniprot. Les gènes qui ne sont pas annotés par leur symbole HGNC garderons alors leur Ensembl ID.

## 2.7.2 Annotation des lncRNAs dans la matrice

L'annotation de notre matrice consiste à ajouter une colonne indiquant si le gène est un lncRNA ou non. Pour cela, nous utilisons un script python (https://forgemia.inra.fr/sarah.maman-haddad/matrix-lncrna-enrichment) développé par Sarah Maman que j'ai pu tester et déboguer à l'aide de fichiers tests légers que j'ai généré avant de l'utiliser sur nos données. Ce script python récupère tous les identifiants des lncRNAs déjà connus comme tels dans notre GTF de référence, ainsi que les identifiants des lncRNAs prédits par FEELnc. Il a donc les identifiants de ceux qui ont déjà un Ensembl ID car ils sont connus d'Ensembl mais pas comme lncRNAs et de ceux qui auront comme identifiant MSTRG.X, correspondant à des nouveaux transcrits lncRNAs. Le script récupère cette liste d'identifiants et les cherche dans la colonne des Ensembl ID de notre matrice de comptage. Quand il retrouve l'identifiant de la liste dans notre matrice, il ajoute la mention "lncRNA". Ceux qui n'auront pas la mention lncRNA auront un "-" dans la colonne.
## 2.8 Analyse différentielle

L'analyse différentielle est réalisée sur R en suivant le Workflow rnaseqGene (Love et al, 2016). Il permet de générer différentes étapes pour la vérification des données de séquences puis pour l'analyse d'expression différentielle avec le package DESeq2 (version 1.34.0. ; Love et al, 2014). Il se sert aussi du package apeglm (version 1.18.0) qui permet de modérer la variation d'expression des gènes à faible niveau de comptage ou avec des comptages variables. Avec un ajustement GLM (Generalized linear model), ce package élimine le problème de la surestimation des changements d'expression pour des faibles comptages (Love et al, 2014).

Ce workflow comprend une partie "Analyse exploratoire et visualisation" qui nécessite une transformation VST (Variance stabilizing transformation) des comptages pour visualiser les échantillons et leurs relations. Les données transformées permettent de générer une carte thermique (Heatmap) et/ou une ACP (analyse en composante principale) afin d'explorer la corrélation entre nos échantillons et l'effet du traitement par rapport au contrôle. C'est une première vérification de la qualité de notre expérience.

Le package R DESeq2 permet la détection de gènes différentiellement exprimés à l'aide de modèles linéaires généralisés binomiaux négatifs, des estimations de la dispersion et de la variation logarithmique en base 2 du changement d'expression (log2FoldChange : LFC). Ainsi, DESeq2 fournit en sortie le LFC et la p-value ajustée (valeur de la statistique du test corrigée pour la multiplicité des tests). Ce sont ces deux valeurs qui nous permettent de déterminer si un gène est différentiellement exprimé, dans quel sens, et si c'est significatif.

D'un point de vue technique, nous avons appliqué la fonction DESeqDataSetFromMatrix de DESeq2, qui nous a permis d'obtenir un objet DESeqDataSet. C'est un objet de classe personnalisé qui sert à stocker les données nécessaires et vérifie qu'elles sont cohérentes. Les paramètres pour le créer sont : la matrice, un "coldata" contenant les informations de l'expérience (individu, traitement) et le design de l'expérience. Ici le design est :

~~ Pig + Treatment ou ~ Cow + Treatment
Cela signifie que nous voulons tester l'effet du traitement en prenant en compte nos animaux.
Comme les échantillons contrôles et traités proviennent des cellules d'un même animal, les
échantillons sont dits appariés. Donc ce modèle permet de contrôler l'effet de l'appariement
de nos données.

Une fois que l'objet est créé, un deuxième filtre sur les gènes est appliqué. Nous
gardons les gènes pour lesquels la somme des comptages bruts de dix échantillons est
supérieure à 1. Nous obtenons les gènes exprimés dans notre tissu d'intérêt : 16 178 chez la
truie et 15 163 chez la vache. C'est sur ces résultats que le pipeline DESeq() pourra être lancé
afin d'obtenir les résultats de l'analyse différentielle.

Pour la visualisation des résultats de l'analyse différentielle, nous avons réalisé un
volcano-plot avec le package EnhancedVolcano (version 1.12.0). Il représente un nuage de
point avec le log2FoldChange en abscisse et la significativité (-Log10P) en ordonnée.

## 2.9 Vérification d'un MSTRG : IGV

Lorsqu'un transcrit de novo est identifié comme différentiellement exprimé, il est
important de vérifier s'il correspond à une réalité biologique. Pour cela, nous avons utilisé
IGV (Integrative Genome Viewer, Robinson *et al*, 2011) au travers d'un script JavaScript
html. Il permet, facilement, de visualiser les transcrits sous IGV sans avoir à installer l'outil
en local. Nous avons précisé à l'intérieur de ce script, le chemin d'accès aux fichiers suivants
: génome de référence (au format fasta et au format indexé fai), un merge de tous les fichiers
BAM (et le fichier merge des BAM indexés), le GTF de référence et enfin le GTF des
transcrits reconstruits par StringTie.

## 2.10 Analyse de la Gene Ontology

A partir d'une liste de gènes différentiellement exprimés, il faut réussir à tirer des
informations pertinentes biologiquement à partir de ces données transcriptomiques. Pour cela,
nous avons réalisé l'analyse ontologique des gènes (GO pour Gene Ontology) permettant de
définir les propriétés majoritaires et les relations entre gènes. La base de données Gene Ontology (The Gene Ontology Consortium et al, 2021) est la base de données de référence mondiale pour la fonction des gènes. Elle est élaborée comme un graphe orienté acyclique où chaque terme est en relation avec d'autres termes du même domaine biologique ou d'un autre domaine. Il y a trois grands domaines : les composants cellulaires, les fonctions moléculaires et les processus biologiques. Nous réalisons ici une analyse de sur-représentation des termes GO les plus significatifs parmi notre liste de gènes différentiellement exprimés.

Pour mener à bien notre analyse ontologique, nous utilisons **Cluster Profiler** (version 4.2.2; Wu et al, 2021), un package R qui implémente différentes méthodes pour l'analyse des gènes. Il prend notamment en charge l'**ORA** (Over Representation Analysis; Boyle et al, 2004) avec sa fonction enrichGO(). L'**ORA** est une approche qui nous permet spécifiquement de déterminer si des termes GO (et donc des fonctions moléculaires ou des processus biologiques) sont surreprésentés dans notre liste de gènes différentiellement exprimés.

Pour utiliser l'ORA, nous avons besoin d'une liste de gènes "background" (N) qui est constituée de l'ensemble des gènes exprimés que nous étudions. Parmi eux, n gènes sont différentiellement exprimés (c'est le nombre de DEG total). Parmi nos gènes "background" se trouvent des groupes de gènes (M) associés à un terme GO particulier (S). Dans cet ensemble de gènes (M) associé à (S), il y a aura k gènes différentiellement exprimés. Pour déterminer si ce terme GO est significativement surreprésenté et que cela n'est pas le fait du hasard, l'ORA procède au calcul de la pvalue p. Dans Cluster Profiler, p est calculé à l'aide de la distribution hypergéométrique de la manière suivante :

$$p = 1 - \sum_{i=0}^{k-1} \frac{\binom{M}{i} \binom{N-M}{n-i}}{\binom{N}{n}}$$

Un "gène ratio" pour un terme GO (S) est calculé correspondant à la formule k/n (nombre de gènes différentiellement exprimés se retrouvant dans un ensemble de gènes (M) dont leur fonction moléculaire est annotée par un même terme GO divisé par le nombre de total de gènes différentiellement exprimés). Plus ce ratio est grand, plus il y a de gènes différentiellement exprimés associé à ce terme GO.
Enfin, **Cluster Profiler** utilise des objets Orgdb ("Organism database" objets contenant les annotations à l'échelle du génome) comme référence pour retrouver nos gènes dans la gene ontology. Il faut alors déterminer quel Orgdb nous utilisons.
<figure><img src="image_5.png" /></figure>

**Figure 4 :** Score de qualité moyenne des lectures pour tous les échantillons chez la vache (à gauche) et la truie (à droite). Ces graphes nous montrent une bonne qualité moyenne des FASTQ;

<figure><img src="image_6.png" /></figure>

**Figure 5 :** Pourcentage d'alignement des lectures sur le génome de référence par échantillon chez la vache (à gauche) et la truie (à droite). Les pourcentages d'alignement sont bons quand ils dépassent 90%. Nous pouvons donc confirmer que notre séquençage est bien réalisé.
# 3. Résultats

## 3.1 Contrôle qualité

Grâce à **MULTIQC**, nous pouvons observer le score de qualité moyen pour chaque échantillon (cf. Figure 4). Nous confirmons que nous pouvons avoir confiance en notre séquençage et qu'il n'y a pas d'échantillon aberrant (outlier) à supprimer de nos analyses.

La seconde chose importante que nous vérifions est la qualité de l'alignement réalisé par **STAR**. Un alignement est considéré comme bon s'il dépasse 90% de lectures alignées sur le génome de référence. Dans notre analyse, nos lectures sont correctement alignées, nous trouvons une moyenne d'alignement de 96,4% chez la vache et de 94,6% chez la truie (cf. Figure 5).

## 3.2 Matrice de comptage

La matrice de comptage est une matrice contenant les gènes en ligne *i* et les noms des échantillons en colonnes *j*. Un élément *i,j* est un comptage du nombre de lectures alignées sur ce gène dans cet échantillon. Ce nombre de lectures est proportionnel à l'expression du gène dans la condition "traité" ou "contrôle".

La matrice de comptage pour les vaches contient 28 028 gènes. Les gènes qui sont non exprimés c'est-à-dire ceux pour lesquels la somme de leurs colonnes est égale à zéro sont retirés. Il en reste 20 190 pour la vache. La matrice de comptage pour les truies contient 32 178 gènes initialement et 23 476 après le filtre.
Méthode shuffle

<figure><img src="image_7.png" /></figure>

Méthode Intergenic

<figure><img src="image_8.png" /><figcaption>Figure 6: Courbe ROC représentant le score de potentiel codant en fonction des performances : exemple chez le bovin.</figcaption></figure>

La courbe à gauche est celle pour la méthode shuffle et la courbe à droite est celle pour la méthode intergenic.
La ligne horizontale (0,98 pour shuffle, 0,989 pour intergenic) représente le seuil de performance atteint.
Lorsque la courbe de la spécificité (en rouge) croise cette droite, nous obtenons le seuil de score de potentiel codant pour lequel les transcrits ayant un score supérieur à celui-ci sont des mRNA (0,479 chez pour, 0,554 pour intergenic). Lorsque la courbe de la sensibilité (en bleue) croise le seuil de performance, nous obtenons le seuil de score de potentiel codant pour lequel les transcrits ayant un score inférieur à celui-ci sont des lncRNA (0,31 pour shuffle, 0,554 pour intergenic). Au milieu des deux seuils, ce sont les transcrits ambigus (TUCp).

<figure><img src="image_9.png" /><figcaption>Figure 7 : Diagramme de Venn représentant la comparaison de lncRNAs prédits par FEELnc par rapport à ceux déjà connus dans Ensembl chez la vache (à gauche) et chez la truie (à droite). Nous comparons la liste des gènes id unique du fichier résultant de FEELnc et les gènes id avec un biotype lncRNA dans le fichier GTF de référence d'Ensembl.</figcaption></figure>
## 3.3 FEELnc : Choix de la méthode et des paramètres

Lors de l'étape 1 de FEELnc, nous pouvons choisir de laisser ou non les mono-exons. J'ai donc lancé avec et sans mono-exons et j'ai comparé les résultats. Enlever les mono-exons faisait perdre un nombre non négligeable de lncRNAs potentiels (plus de trois quarts), nous avons donc choisis de les garder.

FEELnc lance les prédictions avec les deux méthodes : shuffle et intergenic. Parmi les résultats se trouve des courbes ROC (receiver operating characteristic)(cf. Figure 6) qui nous permettent d'évaluer les performances et de faire un choix entre ces deux méthodes. Nous choisissons la méthode qui donnera le moins de transcrits ambigus et avec les meilleures performances.

Pour les bovins (cf. Figure 6) le seuil de performance est de 0,989 pour la méthode intergénique ce qui supérieur à 0,98 pour shuffle. Le seuil déterminant les lncRNAs pour le modèle intergenic est de 0,554, bien supérieur au modèle shuffle qui est de 0,31. De plus, il n'y a pas de transcrit ambigu pour la méthode intergénique. Pour les porcins, le seuil de performance est de 0,98 pour les deux méthodes. Le seuil déterminant les lncRNAs pour intergenic est de 0,42 ; lui aussi bien supérieur à shuffle qui est de 0,254. Il y a des transcrits ambigus entre 0,42 et 0,444 ce qui est inférieur à l'intervalle pour la méthode shuffle (0,245-0,447). Pour nos deux espèces, la méthode choisie est donc la méthode intergenic, qui donne des résultats moins ambigus pour nos deux espèces d'intérêt.

## 3.4 FEELnc : Analyse statistique des résultats

### 3.4.1 Résultat de la prédiction des lncRNAs

Dans le fichier gtf contenant les lncRNA candidats chez l'espèce porcine, il y a 25 537 transcrits de lncRNAs prédits représentant 7 573 lncRNAs uniques. Chez l'espèce bovine, il y a 8146 transcrits prédits pour 3 111 lncRNAs avec un identifiant de gène unique.

### 3.4.2 Proportion de lncRNAs Ensembl vs ceux prédits par FEELnc

Une autre statistique intéressante est de regarder ce que nous apporte FEELnc par rapport aux lncRNAs déjà annotés dans Ensembl. Le GTF de référence contient 1 480 lncRNAs annotés chez la vache et 6 790 annotés chez la truie (cf. Figure 7). Il y a 3 111 lncRNAs prédits par FEELnc chez la vache dont 2285 nouveaux et 8146 chez la truie dont 2817 potentiellement nouveaux.
<figure><img src="image_10.png" /><figcaption>**Figure 8 : Analyse en Composantes Principales pour la truie** avec les deux premières composantes principales. Il y a sous forme de triangle les échantillons traités par BMP15 et en rond ceux contrôles. Les couleurs correspondent aux individus. Nous avons, pour chaque axe, le pourcentage de variance phénotypique expliquée.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9 : Volcano plot des gènes chez la vache.** La significativité statistique (-Log10P) est exprimée en fonction de la variation d'expression (Log2FoldChange). Chaque point sur le graphe correspond à un gène. S'il est de couleur verte il a une valeur absolue du Log2FoldChange >1; s'il est de couleur bleue, a une p-value significative (p<0,05); s'il est de couleur rouge, cumule les deux critères de significativité (= DEGs). Enfin, en gris sont visualisés les gènes non significatifs.</figcaption></figure>
### 3.5 Analyse différentielle

J'ai réalisé une ACP (cf. Figure 8) avec les deux premières composantes principales. Le pourcentage de variance phénotypique expliquée par le traitement BMP15 sur le premier axe est respectivement de 48 et 30 % pour l'espèce porcine et bovine. Cette ACP nous montre clairement une différence entre les échantillons traités et les contrôles, traduisant l'effet de la molécule BMP15. Elle met aussi en évidence l'appariement de nos échantillons.

L'analyse différentielle révèle un total de 4 570 gènes significatifs chez la vache et 7 882 chez la truie (p-value ajustée < 0,05). Nous avons également ajouté un filtre sur le log2FoldChange : en prenant une valeur absolue du log2FoldChange >1 nous obtenons les gènes qui sont deux fois plus (LFC > 1) exprimés (dit up-régulés) ou deux fois moins (LFC < -1) exprimés (gènes dit down-régulés). Le sens de variation s'exprime par rapport à la condition contrôle. Ainsi, nous obtenons 436 DEGs dont 195 down-régulés et 241 up-régulés chez la vache (cf. Figure 9). Tandis que chez la truie, nous trouvons 1008 DEGs dont 517 down-régulés et 491 up-régulés. Nous pouvons déjà voir qu'il y a plus de gènes différentiellement exprimés dans la granulosa chez la truie après action du BMP15 par rapport à la vache.

Parmi ces gènes différentiellement exprimés, il y a 139 lncRNAs significatif chez la vache dont 29 avec une valeur absolue du LFC>1. Treize lncRNAs sont down-régulés et 16 sont up-régulés. Chez la truie, il y a 403 lncRNAs significatif dont 78 avec une valeur absolue du LFC>1. Parmi ceux-ci, 31 sont down-régulés et 47 sont up-régulés.
<figure><img src="image_12.png" /><figcaption>**Figure 10 :** *Captures d'écran de l'interface IGV réalisé avec un script JavaScript HTML. IGV permet la visualisation de la localisation génomique (10:13,666,712-13,807,335) du MSTRG.803 (en haut), différentiellement exprimé chez le bovin. Il englobe un gène annoté du GTF de référence ENSBTAG00000049288. Il est cependant plus long. Le MSTRG.12918 (en bas) est en position 4:79,250,035-79,329,906. Il se superpose à deux transcrits déjà annotés dans Ensembl et semble les relier.*</figcaption></figure>
3.6 Exemple de vérification de transcrits “de novo” d'intérêts

Nous avons vérifié la réalité biologique de deux MSTRG (les 2 plus up-régulés chez l'espèce bovine) en les explorant visuellement à l'aide d'IGV (cf. Figure 10). Nous avons entré la position du MSTRG.803 dans IGV et dans Ensembl. Sa position nous indique qu'il chevauche un ENSBTAG00000049288, un lncRNA et qu'il est encadré par deux gènes : SMAD6 et SMAD3. L'ENSBTAG00000049288 est un lncRNA qui est différentiellement exprimé dans les résultats de DESeq2, bien qu'il ne présente pas autant de comptage que le MSTRG.803. Il se superpose au lncRNA et continue encore pendant des dizaines de milliers de paires de bases.

Bien que nous ne pouvons pas infirmer ou confirmer la réalité biologique de ce MSTRG, il nous a permis de nous intéresser à l'ENSBTAG00000049288. Le fait qu'il se retrouve proche de SMAD6, gène lui aussi différentiellement exprimé et présent sur le volcano-plot, a attiré notre attention. Il a été montré que SMAD6 est un régulateur de l'activité de la famille des TGF-β et inhibe préférentiellement la voie des BMPs, en agissant comme un leurre de SMAD4 (Goto *et al*, 2007).

Le MSTRG.12918 est superposé à deux transcrits déjà annotés par Ensembl et il semble les relier. Il remet en question la réalité biologique de certains transcrits reconstruits par StringTie qui semblent exagérément longs et dont les comptages sont hauts.
<figure><img src="image_13.png" /><figcaption>**Figure 11:** DotPlot résultant de l'ORA montrant le top 7 des termes GO les plus représentés dans le domaine des processus biologiques chez les gènes communs entre la vache et la truie.</figcaption></figure>

Notre liste de gènes d'intérêt ici sont les gènes différentiellement exprimés en communs chez les deux espèces. Nous regardons si un même gène se comporte différemment en fonction de l'espèce sous l'action du BMP15. Les points correspondent au gène ratio. Plus ils sont gros, plus le ratio est important et plus les points sont rouges, plus ils sont significatifs.

<figure><img src="image_14.png" /><figcaption>**Figure 12:** DotPlot des 10 processus biologiques les plus significatifs résultant de l'ORA pour les DEGs truie et vache.</figcaption></figure>

Les points correspondent au gène ratio. Plus ils sont gros, plus le ratio est important et plus les points sont rouges, plus ils sont significatifs.
## 3.7 Analyse de la Gene Ontology

Tout d'abord pour configurer la fonction enrichGo() de clusterProfiler, j'ai dû choisir l'Orgdb le plus adapté. Nous avions le choix entre l'Orgdb de l'espèce (*Sus Scrofa* ou *Bos Taurus*) ou bien celui de l'homme (*Homo Sapiens*). Le nombre de gènes annotés en terme GO est de 43 894 pour l'*Homo Sapiens*, 31 175 pour *Sus Scrofa* et 25 566 pour Bos Taurus. Après avoir effectué les deux analyses c'est à dire en projetant sur l'*Homo Sapiens* ou en utilisant les Orgdb de l'espèce analysée, nous nous sommes rendus compte qu'il était plus judicieux d'utiliser l'OrgDB Humain.

Dans un premier temps, nous avons identifié les gènes communs entre la vache et la truie (cf. Figure 11) et nous avons déterminé quels étaient les processus biologiques (BP) communs partagés par les deux espèces, après un traitement au BMP15, et qui dessinerait ainsi la réponse “générique” à BMP15. A partir de la figure 11, nous observons que la GO “response to BMP” est présente, ce qui démontre que notre traitement a bien marché et que les gènes de réponse à BMP15 sont bien activés. Nous observons aussi trois GO en lien avec la matrice extracellulaire ou l'environnement de la cellule (“extracellular matrix organisation”; “extracellular structure organisation” et “external encapsulating structure organisation”) et un GO “regulation of system process” qui traduisent l'implication de BMP15 dans les relations cellules-cellules et leur interaction avec les composants de la matrice extra-cellulaire (Sengle et al., 2011). Nous retrouvons aussi le GO:“cardiac chamber development” impliqué dans la maturation d'une structure cardiaque, qui peut traduire ici l'effet de BMP15 sur la différentiation des cellules (Persani *et al*, 2014).

Dans un second temps, à partir des DEGs (LFC>1 et padj <0.05) chez les deux espèces, nous avons identifié les 10 termes GO BP les plus significatifs (cf. Figure 12). Le GO: “epithelial cell proliferation” est commun aux deux espèces et retranscrit l’effet prolifératif de BMP15 (Otsuka *et al*, 2000). Nous remarquons que 3 GO en lien avec la maturation ou le développement des cellules, ne sont enrichis que chez l'espèce bovine : “outflow tract morphogenesis”, “mesenchyme development” et “axon development”. Nous observons que 5 GO sont en lien avec la division cellulaire notamment au niveau de la ségrégation des chromosomes ou des chromatides soeurs : “mitotic sister chromatid segregation”, “mitotic nuclear division”, “chromosome segregation”, “sister chromatid
segregation" et "nuclear chromosome segregation". Ses 5 BP sont significativement enrichis que chez l'espèce porcine. Ainsi, cette première analyse des Gene Ontology, nous a permis d'identifier des processus biologiques communs aux deux espèces et surtout de mettre en évidence des spécificités à chaque espèce. Ces premiers résultats encourageants, nous incitent à aller plus loin dans les analyses et notamment d'aller regarder quels sont les réseaux de gènes impliqués dans la réponse commune à BMP15 et surtout les réponses spécifiques de chaque espèce.

# 4. Discussion et perspectives

## 4.1 Impact sur le stage du délai de récupération des FASTQ au NCBI

A mon arrivée au sein du laboratoire, je devais reprendre les données du projet MONOPOLY afin d'y appliquer le pipeline nf-core/rnaseq qui n'avait pas encore été utilisé pour ce projet. Il était de plus en plus utilisé par les plateformes bio-informatiques tel que SIGENAE pour d'autres projets. Cependant, les données MONOPOLY ont été inaccessibles pendant tout le début du stage (jusqu'en mars 2022). La première raison était la version obsolète de **sratoolkit** de notre côté. Il a fallu demander l'installation de la mise à jour au cluster de calcul Genotoul ce qui a été résolu en quelques minutes. Une fois la mise à jour prête, les fichiers n'étaient quand même pas accessibles sur les serveurs du NCBI. Le débogage s'est fait du côté NCBI après les avoir contactés. Cela a évidemment impacté le temps restant pour l'approfondissement biologique du sujet.

Durant ce temps de latence, j'ai pu lancer le **pipeline nf-core/rnaseq** et **FEELnc** sur les données d'un autre projet : AMHAROC (L'hormone anti-Müllérienne ovarienne : régulation, activité et implication dans le syndrome des ovaires polykystiques). J'ai aussi pris en main les outils sur R pour l'analyse différentielle et l'analyse de la Gene Ontology. Tout cela m'a permis de préparer mes scripts pour un lancement plus rapide sur les données MONOPOLY.

## 4.2 Difficulté pour la quantification des transcrits "de novo"

Le **pipeline nf-core/rnaseq** est un très grand progrès pour la bio-informatique. Pouvoir lancer un pipeline pour l'analyse RNA-seq en une seule commande est un avantage
considérable. Cependant, nous nous sommes rendues compte qu'il était très adapté pour la génération de la matrice de comptage et l'étude des gènes mais pas forcément pour les lncRNAs. Le **pipeline nf-core/rnaseq** implique de faire des traitements en dehors tel que la quantification des transcrits de novo. Elle est faite avec **StringTie** qui n'a pas le même algorithme de comptage que **RSEM**.

A ce sujet, nous avons contacté un des développeurs du pipeline, H.patel, qui nous a confirmé que c'était un choix de ne pas intégrer les transcrits de novo à la matrice de comptage et que leur comptage n'était pas encore codé. Les développeurs se sont montrés intéressés par les solutions apportées par mon stage.

Bien que nous ayons trouvé des solutions pour au final avoir une quantification, l'analyse des MSTRG reconstruits montrent qu'ils ne respectent pas toujours les réalités biologiques. Ils sont parfois très longs ou chevauchent plusieurs gènes déjà existants (cf Figure 10). Il faut donc s'attendre à une analyse très approfondie pour leurs vérifications. Cette vérification demande du temps, il faut donc l'envisager uniquement pour les lncRNAs différentiellement exprimés par BMP15.

Un autre problème que nous avons observé est la sur-expression des MSTRG dans la matrice de comptage avec l'utilisation de **StringTie**. Leurs comptages bruts sont nettement supérieurs à ce que nous attendions pour un lncRNA. Nous avons pourtant comparé les comptages de gène commun entre la matrice **StringTie** et la matrice **RSEM** et nous avions constaté que les comptages étaient suffisamment similaires pour être assemblés dans une même matrice. Pour essayer de résoudre cette problématique, une des pistes est de relancer une quantification des MSTRG mais cette fois avec **RSEM**. Nous relançons le **pipeline nf-core/rnaseq** avec le GTF de référence qui sert à la quantification des transcrits de novo. Ces résultats n'ont pas encore fini d'être étudiés mais montrent une très légère amélioration sur les comptages.

Une autre piste pourrait être l'utilisation du pipeline **TAGADA** (https://github.com/FAANG/analysis-TAGADA). C'est aussi un pipeline **Nextflow** qui traite les données RNA-seq. Il s'inscrit dans le cadre du projet GENE-Switch, financé par le programme de recherche et d'innovation Horizon 2020 de l'Union Européenne. Il est différent car, par exemple, ce pipeline utilise une quantification avec StringTie et lance **FEELnc** directement à la suite. Il merge les GTF des transcrits reconstruits avec **Tmerge**, les quantifie avec **StringTie** et permet d'avoir un seul GTF avec tous les nouveaux transcrits potentiels dont ceux provenant de **FEELnc**. Pour une étude plus approfondie sur les
lncRNAs, ce pipeline semble plus adapté. Nous avons lancé ce pipeline courant avril, mais il a fallu beaucoup de tests avant d'arriver à un résultat sans erreur. Nous avons donc disposé des résultats fin-mai. Il restera donc à comparer les résultats de ce pipeline avec ceux du **pipeline nf-core/rnaseq** pour confirmer ou infirmer si ces résultats sont plus intéressants pour notre étude.

## 4.3 Confiance apportée aux prédictions de FEELnc

Lorsque nous utilisons **FEELnc** nous avons confiance en sa prédiction. S'il prédit que c'est un lncRNA, nous ne remettons généralement pas en cause ce résultat. C'est important à noter car s'il retrouve un lncRNA qui n'avait pas ce biotype initialement dans Ensembl, il faut faire un choix pour déterminer qui a raison. Ce choix ici est de suivre **FEELnc**.

Il existe un nouveau logiciel : **IncEvo** qui propose des fonctionnalités qui seraient extrêmement intéressantes à essayer notamment pour la prédiction et surtout pour l'étude de l'évolution des lncRNAs entre espèces. C'est un pipeline qui se compose de trois "modules" : l'assemblage du transcriptome à partir des données RNA-Seq, la prédiction des lncRNAs puis l'étude de conservation entre espèces (Bryzghalov *et al*, 2021). Nous pourrons comparer ses résultats avec FEELnc et cela nous apporterait surtout une analyse comparative évolutive des lncRNAs entre l'espèce bovine et porcine.

## 4.4 Réseaux de gènes

Un des objectifs ultimes d'une analyse transcriptomique est d'aboutir à la mise en évidence des réseaux de gènes gouvernant le phénotype d'intérêt. Ces réseaux peuvent être mis en évidence par l'étude de la co-expression des gènes. L'hypothèse sous-jacente est que les gènes qui sont co-exprimés avec d'autres gènes appartiennent à la même fonction biologique ou voie de signalisation. Il est également possible d'établir la liste des facteurs de transcription impliqués chez nos deux espèces afin de les utiliser pour déduire les régulations existantes dans notre modèle et identifier les réseaux régulateurs de gènes spécifiques ou communs à chaque espèce.

Afin de répondre à ces problématiques, il est possible d'utiliser le package R **Corto** (version 1.1.10; Mercatelli *et al*, 2020), qui implémente une combinaison de corrélation de Pearson et d'inégalité de traitement des données (DPI)(Reverter & Chan, 2008) afin de construire des réseaux entre des centroïdes et des cibles. Ici, les centroïdes sont des facteurs de transcription et leurs cibles sont des gènes. Il peut mettre en evidence une relation, une
co-expression entre gènes et facteurs de transcription et donc en déduire des mécanismes de régulation ou de fonction géniques.

En intégrant les lncRNAs dans ces réseaux, nous pourrions aussi déterminer si certains ont un rôle de régulateur de la transcription. Ainsi, si un ou plusieurs mRNA sont co-exprimés avec des lncRNAs, cela peut s'expliquer par le fait que les mRNA constituent les gènes cibles en "trans" du lncRNA. Nous utiliserons ainsi les réseaux régulateurs de gènes pour identifier les cibles potentielles en "trans" des lncRNAs que nous validerons avec l'outil lncTar (Li et al. 2014) qui prédit et quantifie l'interaction entre lncRNA et mRNA, en calculant l'énergie libre de la structure secondaire se formant entre les deux molécules.

Un autre package qui semble intéressant est le package Cetf (version 1.6.0.; Oliveira de Biagi *et al*, 2021). Il contient les fonctions permettant d'exécuter deux algorithmes: L'algorithme du coefficient de corrélation partielle avec la théorie de l'information (PCIT; Reverter & Chan, 2008) et l'algorithme des facteurs d'impact des régulateurs (RIF; Reverter *et al*, 2010). Il permet, en résumé, l'identification des corrélations significatives pour définir des réseaux de co-expression génique ainsi que l'identification des facteurs de transcription.

Enfin, une fois les réseaux réalisés sur R, l'objectif sera de les relier avec Cytoscape (version 3.9.1; Shannon *et al*, 2003) un logiciel permettant la visualisation des réseaux et de les rendre interactifs.

Pour conclure ce stage a permis d'identifier les gènes différentiellement exprimés chez deux espèces d'intérêt agronomique suite à un traitement par BMP15. Une première caractérisation des gènes impliqués dans la réponse à BMP15 dans les cellules de la granulosa a été entreprise et a mis en évidence des effets différents de BMP15 selon que l'espèce soit mono-ovulante ou poly-ovulante. Il a permis de tester, d'identifier les défauts et d'apporter des perspectives et des solutions pour réussir à obtenir un premier catalogue de lncRNAs dans les cellules de granulosa chez la vache et la truie et d'identifier les lncRNAs différentiellement exprimés entre ces deux espèces.
Références

Bardou P, Mariette J, Escudié F, Djemiel C & Klopp C (2014) jvenn: an interactive Venn diagram viewer. *BMC Bioinformatics* 15: 293

Bhan A, Soleimani M & Mandal SS (2017) Long Noncoding RNA and Cancer: A New Paradigm. *Cancer Res* 77: 3965–3981

Boyle EI, Weng S, Gollub J, Jin H, Botstein D, Cherry JM & Sherlock G (2004)
GO::TermFinder--open source software for accessing Gene Ontology information and
finding significantly enriched Gene Ontology terms associated with a list of genes.
*Bioinforma Oxf Engl* 20: 3710–3715

Bryzghalov O, Makałowska I & Szcześniak MW (2021) lncEvo: automated identification and conservation study of long noncoding RNAs. *BMC Bioinformatics* 22: 59

Demars J, Fabre S, Sarry J, Rossetti R, Gilbert H, Persani L, Tosser-Klopp G, Mulsant P,
Nowak Z, Drobik W, et al (2013) Genome-wide association studies identify two novel
BMP15 mutations responsible for an atypical hyperprolificacy phenotype in sheep.
*PLoS Genet* 9: e1003482

Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M & Gingeras TR (2013) STAR: ultrafast universal RNA-seq aligner. *Bioinformatics* 29: 15–21

Fabre S, Pierre A, Mulsant P, Bodin L, Di Pasquale E, Persani L, Monget P & Monniaux D (2006) Regulation of ovulation rate in mammals: contribution of sheep genetic models. *Reprod Biol Endocrinol RBE* 4: 20

Feuerstein P, Cadoret V, Dalbies-Tran R, Guérif F & Royère D (2006) Le dialogue ovocyte–cumulus. *Gynécologie Obstétrique Fertil* 34: 793–800

Galloway SM, McNatty KP, Cambridge LM, Laitinen MPE, Juengel JL, Jokiranta TS,
McLaren RJ, Luiro K, Dodds KG, Montgomery GW, et al (2000) Mutations in an
oocyte-derived growth factor gene (BMP15) cause increased ovulation rate and
infertility in a dosage-sensitive manner. *Nat Genet* 25: 279–283

Goto K, Kamiya Y, Imamura T, Miyazono K & Miyazawa K (2007) Selective Inhibitory Effects of Smad6 on Bone Morphogenetic Protein Type I Receptors. *J Biol Chem* 282: 20603–20611

Gudenas BL, Wang J, Kuang S, Wei A, Cogill SB & Wang L (2019) Genomic data mining for functional annotation of human long noncoding RNAs. *J Zhejiang Univ-Sci B* 20: 476–487

Hanrahan JP, Gregan SM, Mulsant P, Mullen M, Davis GH, Powell R & Galloway SM
(2004) Mutations in the Genes for Oocyte-Derived Growth Factors GDF9 and BMP15 Are Associated with Both Increased Ovulation Rate and Sterility in Cambridge and Belclare Sheep (Ovis aries)1. *Biol Reprod* 70: 900-909

Harshil Patel, Ewels P, Peltzer A, Hammarén R, Botvinnik O, Sturm G, Moreno D, Pranathi Vemuri, Silviamorins, Pantano L, et al (2021) nf-core/rnaseq: nf-core/rnaseq v3.5 - Copper Chameleon Zenodo

Huber SD <Biomartdev@Gmail C Wolfgang (2017) biomaRt Bioconductor

Jozkowiak M, Hutchings G, Jankowski M, Kulcenty K, Mozdziak P, Kempisty B, Spaczynski RZ & Piotrowska-Kempisty H (2020) The Stemness of Human Ovarian Granulosa Cells and the Role of Resveratrol in the Differentiation of MSCs—A Review Based on Cellular and Molecular Knowledge. *Cells* 9: 1418

Juengel JL, Cushman RA, Dupont J, Fabre S, Lea RG, Martin GB, Mossa F, Pitman JL, Price CA & Smith P (2021) The ovarian follicle of ruminants: the path from conceptus to adult. *Reprod Fertil Dev* 33: 621

Li B & Dewey CN (2011) RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. *BMC Bioinformatics* 12: 323

Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R, & 1000 Genome Project Data Processing Subgroup (2009) The Sequence Alignment/Map format and SAMtools. *Bioinformatics* 25: 2078-2079

Love MI, Anders S, Kim V & Huber W (2016) RNA-Seq workflow: gene-level exploratory analysis and differential expression. *F1000Research* 4: 1070

Love MI, Huber W & Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol* 15: 550

Mercatelli D, Lopez-Garcia G & Giorgi FM (2020) corto: a lightweight R package for gene network inference and master regulator analysis. *Bioinformatics* 36: 3916-3917

Moore RK, Otsuka F & Shimasaki S (2003) Molecular Basis of Bone Morphogenetic Protein-15 Signaling in Granulosa Cells. *J Biol Chem* 278: 304-310

Mulsant P, Lecerf F, Fabre S, Schibler L, Monget P, Lanneluc I, Pisselet C, Riquet J, Monniaux D, Callebaut I, et al (2001) Mutation in bone morphogenetic protein receptor-IB is associated with increased ovulation rate in Booroola Mérino ewes. *Proc Natl Acad Sci* 98: 5104-5109

Oliveira de Biagi CA, Nociti RP, Brotto DB, Funicheli BO, Cássia Ruy P de, Bianchi Ximenez JP, Alves Figueiredo DL & Araújo Silva W (2021) CeTF: an R/Bioconductor package for transcription factor co-expression networks using regulatory impact factors (RIF) and partial correlation and information (PCIT)
analysis. *BMC Genomics* 22: 624

Otsuka F, Yamamoto S, Erickson GF & Shimasaki S (2001) Bone Morphogenetic Protein-15 Inhibits Follicle-stimulating Hormone (FSH) Action by Suppressing FSH Receptor Expression. *J Biol Chem* 276: 11387-11392

Otsuka F, Yao Z, Lee T, Yamamoto S, Erickson GF & Shimasaki S (2000) Bone Morphogenetic Protein-15. *J Biol Chem* 275: 39523-39528

Persani L, Rossetti R, Di Pasquale E, Cacciatore C & Fabre S (2014) The fundamental role of bone morphogenetic protein 15 in ovarian function and its involvement in female fertility disorders. *Hum Reprod Update* 20: 869-883

Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT & Salzberg SL (2015) StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. *Nat Biotechnol* 33: 290-295

Ponting CP, Oliver PL & Reik W (2009) Evolution and Functions of Long Noncoding RNAs. *Cell* 136: 629-641

Reverter A & Chan EKF (2008) Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks. *Bioinformatics* 24: 2491-2497

Reverter A, Hudson NJ, Nagaraj SH, Pérez-Enciso M & Dalrymple BP (2010) Regulatory impact factors: unraveling the transcriptional regulation of complex traits from expression data. *Bioinformatics* 26: 896-904

Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G & Mesirov JP (2011) Integrative genomics viewer. *Nat Biotechnol* 29: 24-26

Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B & Ideker T (2003) Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. *Genome Res* 13: 2498-2504

Silva BDM, Castro EA, Souza CJH, Paiva SR, Sartori R, Franco MM, Azevedo HC, Silva T a. SN, Vieira AMC, Neves JP, *et al* (2011) A new polymorphism in the Growth and Differentiation Factor 9 (GDF9) gene is associated with increased ovulation rate and prolificacy in homozygous sheep. *Anim Genet* 42: 89-92

The Gene Ontology Consortium, Carbon S, Douglass E, Good BM, Unni DR, Harris NL, Mungall CJ, Basu S, Chisholm RL, Dodson RJ, *et al* (2021) The Gene Ontology resource: enriching a GOld mine. *Nucleic Acids Res* 49: D325-D334

Vitting-Seerup K (2017) IsoformSwitchAnalyzeR Bioconductor

Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, Feng T, Zhou L, Tang W, Zhan L, *et al* (2021) clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. *The*
Wucher V, Legeai F, Hédan B, Rizk G, Lagoutte L, Leeb T, Jagannathan V, Cadieu E, David A, Lohi H, *et al* (2017) FEELnc: a tool for long non-coding RNA annotation and its application to the dog transcriptome. *Nucleic Acids Res*: gkw1306
## Liste des abréviations

ACP : Analyse en composante principale

ANR : Agence National de la Recherche

BMP15 : Protéine morphogénétique osseuse 15

BMPR1B : Récepteur de protéine morphogénétique osseuse de type 1B

BP : Biological process

CPS : Score de potentiel codant

DEGs : Gènes différentiellement exprimés.

FEELnc : FlExible Extraction of LncRNAs

FSH : Hormone de stimulation folliculaire

GDF9 : Facteur de croissance et de différenciation 9

GDF : Facteurs de croissance et de différenciation

GLM : Modèle linéaire généralisé

GO : Gene Ontology

GTF : Format d'élément général

LFC : Log2FoldChange

LH : Hormone Lutéinisante

lncRNA : ARN long non codant

MF : Molecular Function

mRNA : ARN messager

NCBI : Centre américain pour les informations biotechnologiques

ORA : Analyse de surreprésentation

ORF : Cadre ouvert de lecture

RNA-seq : Séquençage de l'ARN

RSEM : RNA-Seq by Expectation Maximization
STAR : Spliced Transcripts Alignment to a Reference

TGF-β : Facteur de croissance transformant Beta

TPM : transcrit par million

VST : Transformation stabilisant la variance