Croisement des données¶

Claire Toffano-Nioche, I2BC, Paris-Saclay¶
Pauline Francois, Anses, Lyon¶
Selon les ressources de Matthias Zytnicki, INRAE, Toulouse¶
et toute la team Roscoff¶

Objectifs du cours¶

  • Permettre de croiser des données génomiques de types différents
  • Renforcer les notions vues précedemment
  • Se familiariser avec bedtools

Ce que ce cours n'est pas :

  • Une réponse à une vraie question biologique
  • Une méthode statistiquement valide

Croisement de données¶

Qu’est-ce c'est ?¶

  • Une comparaison des positions ou intervalles génomiques:
    • un variant par rapport à des gènes d’intérêt,
    • des pics de ChIP-Seq par rapport à des gènes différentiellement exprimés, etc.

À quel moment est-ce valide ?¶

  • Lorsque vous cherchez des co-occurrences
  • Lorsque vous donnez des distributions (de distance)

À quel moment est-ce douteux ?¶

  • Lorsque les résultats sont présentés comme significatifs.

Avant de commencer¶

Vérifier que vous êtes dans votre espace projet, sous-dossier TP_croisement

In [1]:
pwd

/shared/ifbstor1/projects/2336_nanopore/TP_croisement

Accès aux données du TP¶

  1. Rendez vous dans le dossier data/ et listez en le contenu
In [2]:
cd data
ls

chrs.len  DEgenes.txt  hsGRCh38.genes.gff3  picChip.bed  variant.vcf
  1. Afin de se familiariser avec le contenu des données, on affiche les 2 premières lignes de l'ensemble des fichiers :
In [3]:
head -n 2 *

==> chrs.len <==
1	248956422
10	133797422

==> DEgenes.txt <==
ENSG00000004846
ENSG00000005981

==> hsGRCh38.genes.gff3 <==
1	ensembl_havana	gene	65419	71585	.	+	.	ID=gene:ENSG00000186092;Name=OR4F5;biotype=protein_coding;description=olfactory receptor family 4 subfamily F member 5 [Source:HGNC Symbol%3BAcc:HGNC:14825];gene_id=ENSG00000186092;logic_name=ensembl_havana_gene;version=6
1	ensembl_havana	gene	450703	451697	.	-	.	ID=gene:ENSG00000284733;Name=OR4F29;biotype=protein_coding;description=olfactory receptor family 4 subfamily F member 29 [Source:HGNC Symbol%3BAcc:HGNC:31275];gene_id=ENSG00000284733;logic_name=ensembl_havana_gene;version=1

==> picChip.bed <==
22	16192349	16192565	region_1
22	16846630	16870710	region_2

==> variant.vcf <==
##fileformat=VCFv4.0
##fileDate=20180418
  1. Pour le fichier variant.vcf on ne voit que des commentaires (souvent des lignes qui commencent par "#") ; afficher les deux dernières lignes de ce fichier :
In [4]:
tail -n 2 variant.vcf

Y	26614668	rs376925794	TC	T	.	.	RS=376925794;RSPOS=26614671;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9968,0.003244;COMMON=1
Y	26624486	rs771929773	T	TA	.	.	RS=771929773;RSPOS=26624486;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9968,0.003244;COMMON=1
  1. Activer les outils bedtools (v2.30.0) dont nous aurons besoin plus loin
In [2]:
module load bedtools/2.30.0

Problème¶

Question scientifique¶

Mes gènes différentiellement exprimés contiennent-ils plus de pics de chip qu'attendu ? (oui/non, pas de quantification)

En d'autres termes, y a-t-il plus ou moins de pics de chip-seq dans les gènes différentiellement exprimés par rapport aux gènes non différentiellement exprimés ?

Données¶
  • DEgenes.txt
  • picChip.bed
  • hsGRCh38.genes.gff3
Données¶
  • DEgenes.txt
  • picChip.bed
  • hsGRCh38.genes.gff3

Afficher les 2 premières lignes des fichiers DEgenes.txt, picChip.bed et hsGRCh38.genes.gff3 :

In [1]:
head -n 2 DEgenes.txt picChip.bed hsGRCh38.genes.gff3

==> DEgenes.txt <==
ENSG00000004846
ENSG00000005981

==> picChip.bed <==
22	16192349	16192565	region_1
22	16846630	16870710	region_2

==> hsGRCh38.genes.gff3 <==
1	ensembl_havana	gene	65419	71585	.	+	.	ID=gene:ENSG00000186092;Name=OR4F5;biotype=protein_coding;description=olfactory receptor family 4 subfamily F member 5 [Source:HGNC Symbol%3BAcc:HGNC:14825];gene_id=ENSG00000186092;logic_name=ensembl_havana_gene;version=6
1	ensembl_havana	gene	450703	451697	.	-	.	ID=gene:ENSG00000284733;Name=OR4F29;biotype=protein_coding;description=olfactory receptor family 4 subfamily F member 29 [Source:HGNC Symbol%3BAcc:HGNC:31275];gene_id=ENSG00000284733;logic_name=ensembl_havana_gene;version=1

DEgenes.txt est notre liste de gènes DE :

.

picChip.bed est notre liste de pics issus du chip-seq et positionnés sur le génome :

Pour croiser ces données il faut positionner les gènes DE sur le génome, on fait cela grâce au fichier d'annotation (gènes non DE plus clairs) :

Le positionnement des données sur le génome permet maintenant de les croiser :

Le but sera de compter le nombre gène dans chacune des 4 classes possibles :

Protocole envisagé¶
A. Extraire les intervalles génomiques des gènes différentiellement exprimés (DE)
B. Comparer les intervalles génomiques des gènes DE avec les régions H3K4me3 (ie. l'emplacement des pics de chip-seq)
C. Compter le nombre de chevauchements entre les gènes DE et les pics
D. Trouver les intervalles génomiques des gènes non-différentiellement exprimés
E. Comparer ces intervalles génomiques avec les régions H3K4me3
F. Compter le nombre de chevauchements entre les gènes non DE et les pics
G. Comparer les nombres de chevauchements

A. Extraire les intervalles génomiques des gènes différentiellement exprimés (DE)¶

  1. Chercher (grep) une liste de mots (-f) contenue dans le fichier DEgenes.txt dans le fichier hsGRCh38.genes.gff3 et rediriger la sortie standard dans un fichier nommé DEgenes.gff3
In [4]:
grep -f DEgenes.txt hsGRCh38.genes.gff3 > DEgenes.gff3
head -n2 DEgenes.gff3
echo "-----"
wc -l DEgenes.txt 
wc -l hsGRCh38.genes.gff3
wc -l DEgenes.gff3

1	ensembl_havana	gene	33081104	33120530	.	+	.	ID=gene:ENSG00000142920;Name=AZIN2;biotype=protein_coding;description=antizyme inhibitor 2 [Source:HGNC Symbol%3BAcc:HGNC:29957];gene_id=ENSG00000142920;logic_name=ensembl_havana_gene;version=16
1	ensembl_havana	gene	70852353	71047808	.	-	.	ID=gene:ENSG00000050628;Name=PTGER3;biotype=protein_coding;description=prostaglandin E receptor 3 [Source:HGNC Symbol%3BAcc:HGNC:9595];gene_id=ENSG00000050628;logic_name=ensembl_havana_gene;version=20
-----
162 DEgenes.txt
21492 hsGRCh38.genes.gff3
162 DEgenes.gff3

B. Comparer les intervalles génomiques des gènes DE avec les régions H3K4me3 (ie. l'emplacement des pics de chip-seq)¶

Comme on va utiliser la suite d'outils que propose Bedtools, avant tout, faisons un tour sur la documentation.

Chaque sous-commande fonctionne sur le même principe :

bedtools [sub-command] [options]

Les options varient selon la commande. Il suffit de cliquer sur son nom sur le site pour avoir accès à la fonction de la commande, à ses options et à des exemples.

Les trois qui nous intéresserons dans un premier temps sont intersect, flank et closest

Bedtools accepte quasiment tous les formats de fichiers tabulés (GFF3, BED, etc)

Dans notre cas de figure, nous souhaitons savoir si des gènes DE contiennent des pics et obtenir le nom de ces gènes. Nous utilisons donc la sous-commande "bedtools ....." ?

Le fichier A sera celui de l'annotation des gènes DE (DEgenes.gff3) et le B correspondra aux pics (picChip.bed). Quelle option choisir pour obtenir uniquement les annotations du fichier A qui croisent avec le fichier B ?

Rediriger la sortie standard dans un fichier nommé intersect_DEgenes_picChip.gff3

In [5]:
bedtools intersect -a DEgenes.gff3 -b picChip.bed -wa -u \
         > intersect_DEgenes_picChip.gff3

C. Compter le nombre de chevauchements entre les gènes DE et les pics¶

Nous avons obtenu la liste des gènes dans lesquels se trouvent des pics (intersect_DEgenes_picChip.gff3). Nous avons maintenant besoin de les compter.

In [9]:
wc -l intersect_DEgenes_picChip.gff3

77 intersect_DEgenes_picChip.gff3

D. Trouver les intervalles génomiques des gènes non-différentiellement exprimés¶

grep -v : Chercher l'inverse d'une liste de gènes dans un texte et rediriger la sortie standard dans un fichier nommé notDEgenes.gff3

On récupère donc tous les gènes dont le nom n'est pas dans la liste des gènes DE.

In [10]:
grep -v -f DEgenes.txt hsGRCh38.genes.gff3 > notDEgenes.gff3
wc -l notDEgenes.gff3

21330 notDEgenes.gff3

E. Comparer ces intervalles génomiques avec les régions H3K4me3¶

Faire l'intersection entre le fichier gff3 des gènes non DE et les pics de Chip-seq et écrire la sortie dans un fichier intersect_notDEgenes_picChip.gff3

In [11]:
bedtools intersect -a notDEgenes.gff3 -b picChip.bed -wa -u \
       > intersect_notDEgenes_picChip.gff3
wc -l intersect_notDEgenes_picChip.gff3

9438 intersect_notDEgenes_picChip.gff3

F. Compter le nombre de chevauchements entre les gènes non DE et les pics¶

Afin de pouvoir ensuite faire des ratios "nb de gène avec au moins un chevauchements/nb de gène" nous devons connaitre les valeurs suivantes :

  • Nombre de gènes DE
  • Nombre de gènes DE avec chevauchement
  • Nombre de gènes non DE
  • Nombre de gènes non DE avec chevauchement
In [12]:
wc -l DEgenes.gff3
wc -l intersect_DEgenes_picChip.gff3
wc -l notDEgenes.gff3
wc -l intersect_notDEgenes_picChip.gff3

162 DEgenes.gff3
77 intersect_DEgenes_picChip.gff3
21330 notDEgenes.gff3
9438 intersect_notDEgenes_picChip.gff3

G. Comparer les nombres de chevauchements¶

Résumons nos informations :

Ratios avec pics sans pic total
gene DE 77 162-77 162
gene non DE 9438 21330-9438 21330
total 21492

Par simple calcul de pourcentage nous obtenons le tableau suivant :

Ratios avec pics sans pic total
gene DE 47.53% 52.47% 162
gene non DE 44.24% 55.76% 21330
total 21492

En conclusion : On peut noter une légère hausse du nombre de pics dans les gènes DE vs non DE.

/!\ En revanche on ne peut pas affirmer ou infirmer que cela soit significatif /!\

Bonus n°1¶

Question scientifique¶

Quels sont les variants qui sont dans les promoteurs (2kb en amont du gène) des gènes différentiellement exprimés ?

Données¶
  • DEgenes.gff3 (généré précédemment)
  • variant.vcf
  • chrs.len

Afficher les 2 dernières lignes des fichiers DEgenes.gff3 et variant.vcf :

In [6]:
tail -n2 DEgenes.gff3 variant.vcf

==> DEgenes.gff3 <==
9	havana	gene	133098121	133163914	.	-	.	ID=gene:ENSG00000285245;Name=AL162417.1;biotype=protein_coding;description=novel protein;gene_id=ENSG00000285245;logic_name=havana;version=1
X	ensembl_havana	gene	71107404	71112108	.	-	.	ID=gene:ENSG00000147168;Name=IL2RG;biotype=protein_coding;description=interleukin 2 receptor subunit gamma [Source:HGNC Symbol%3BAcc:HGNC:6010];gene_id=ENSG00000147168;logic_name=ensembl_havana_gene;version=12

==> variant.vcf <==
Y	26614668	rs376925794	TC	T	.	.	RS=376925794;RSPOS=26614671;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9968,0.003244;COMMON=1
Y	26624486	rs771929773	T	TA	.	.	RS=771929773;RSPOS=26624486;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9968,0.003244;COMMON=1
Protocole envisagé¶
A. Extraire la région située à 2kb en amont des gènes différentiellement exprimés (DE)
B. Trouver tous les variants chevauchant les régions trouvées précedemment

.

Nous allons avoir besoin d'une autre sous-commande de bedtools : bedtools flank Regardez bien ses options.

Notons qu'il y a une option obligatoire, -g pour communiquer à l'outil la taille des chromosomes. Pourquoi ?

Lancer la commande et écrire la sortie standard dans un fichier flank2kDEgenes.gff3

In [7]:
bedtools flank -i DEgenes.gff3 -g chrs.len -l 2000 -r 0 -s \
   > flank2kDEgenes.gff3
head -n2 flank2kDEgenes.gff3

1	ensembl_havana	gene	33079104	33081103	.	+	.	ID=gene:ENSG00000142920;Name=AZIN2;biotype=protein_coding;description=antizyme inhibitor 2 [Source:HGNC Symbol%3BAcc:HGNC:29957];gene_id=ENSG00000142920;logic_name=ensembl_havana_gene;version=16
1	ensembl_havana	gene	71047809	71049808	.	-	.	ID=gene:ENSG00000050628;Name=PTGER3;biotype=protein_coding;description=prostaglandin E receptor 3 [Source:HGNC Symbol%3BAcc:HGNC:9595];gene_id=ENSG00000050628;logic_name=ensembl_havana_gene;version=20

Maintenant, c'est du déjà-vu : Trouver tous les variants chevauchant les régions trouvées précedemment. Renvoyer la sortie standard dans un fichier nommé intersect_flank2kDEgenes_variant.txt

N'oubliez pas de regarder les options qui nous intéressent ici et visualisez les 2 premières lignes du fichier quand celui ci est généré

In [10]:
bedtools intersect -a variant.vcf -b flank2kDEgenes.gff3 -wa -u \
    > intersect_flank2kDEgenes_variant.txt
head -n 2 intersect_flank2kDEgenes_variant.txt
echo "-----"
wc -l intersect_flank2kDEgenes_variant.txt

1	33079702	rs10712676	GA	G	.	.	RS=10712676;RSPOS=33079703;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x0501000a0005170126000200;GENEINFO=AZIN2:113451|LOC105378635:105378635;WGT=1;VC=DIV;SLO;INT;R5;ASP;VLD;G5A;G5;GNO;KGPhase3;CAF=0.8706,0.1294;COMMON=1;TOPMED=0.96872610856269113,0.03127389143730886
1	74197393	rs201749320	TTC	T	.	.	RS=201749320;RSPOS=74197394;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x0500000a000504003e000200;GENEINFO=LRRIQ3:127255|FPGT-TNNI3K:100526835|FPGT:8790;WGT=1;VC=DIV;INT;R5;ASP;VLD;KGPhase1;KGPhase3;CAF=0.9958,0.004193;COMMON=1;TOPMED=0.99698967889908256,0.00301032110091743
-----
405 intersect_flank2kDEgenes_variant.txt

Pour aller plus loin avec bedtools intersect, on peut se demander combien y a-t-il de SNP par promoteur ?

Stocker l'information dans un fichier intersect_flanck2kDEgenes_n_variant.txt et visualiser les 2 premières lignes

In [13]:
bedtools intersect -a flank2kDEgenes.gff3 -b variant.vcf -c \
   > intersect_flanck2kDEgenes_n_variant.txt
head -n 2 intersect_flanck2kDEgenes_n_variant.txt

1	ensembl_havana	gene	33079104	33081103	.	+	.	ID=gene:ENSG00000142920;Name=AZIN2;biotype=protein_coding;description=antizyme inhibitor 2 [Source:HGNC Symbol%3BAcc:HGNC:29957];gene_id=ENSG00000142920;logic_name=ensembl_havana_gene;version=16	1
1	ensembl_havana	gene	71047809	71049808	.	-	.	ID=gene:ENSG00000050628;Name=PTGER3;biotype=protein_coding;description=prostaglandin E receptor 3 [Source:HGNC Symbol%3BAcc:HGNC:9595];gene_id=ENSG00000050628;logic_name=ensembl_havana_gene;version=20	0

Bonus n°2¶

Question scientifique¶

Quels sont les gènes différentiellement exprimés les plus proches des pics ChIP ET qui contiennent une mutation ?

Données¶
  • DEgenes.gff3 (généré dans le problème 1)
  • variant.vcf
  • picChip.bed

Afficher les 2 dernières lignes des fichiers DEgenes.gff3, variant.vcf et picChip.bed :

In [14]:
tail -n2 DEgenes.gff3 variant.vcf picChip.bed

==> DEgenes.gff3 <==
9	havana	gene	133098121	133163914	.	-	.	ID=gene:ENSG00000285245;Name=AL162417.1;biotype=protein_coding;description=novel protein;gene_id=ENSG00000285245;logic_name=havana;version=1
X	ensembl_havana	gene	71107404	71112108	.	-	.	ID=gene:ENSG00000147168;Name=IL2RG;biotype=protein_coding;description=interleukin 2 receptor subunit gamma [Source:HGNC Symbol%3BAcc:HGNC:6010];gene_id=ENSG00000147168;logic_name=ensembl_havana_gene;version=12

==> variant.vcf <==
Y	26614668	rs376925794	TC	T	.	.	RS=376925794;RSPOS=26614671;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9968,0.003244;COMMON=1
Y	26624486	rs771929773	T	TA	.	.	RS=771929773;RSPOS=26624486;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9968,0.003244;COMMON=1

==> picChip.bed <==
18	78005736	78005857	region_52421
18	78016442	78016675	region_52422
Protocole envisagé¶
1. Trouver les pics qui contiennent une mutation
2. Trouver le gène le plus proche de chaque région précédemment trouvée

Trouver les pics qui contiennent une mutation¶

Rien de neuf sous le soleil, nous stockons l'output dans un fichier picVariant.bed.

In [15]:
bedtools intersect -a picChip.bed -b variant.vcf -u > picVariant.bed

L'option -u permet de garder l'ensemble de l'information des pics sans ajouter celle des variants

In [16]:
head -n 5 picVariant.bed

22	16846630	16870710	region_2
22	17070620	17114004	region_4
22	17128021	17136091	region_5
22	17147545	17173488	region_6
22	17192673	17206732	region_7

Trouver le gène le plus proche de chaque région précédemment trouvée¶

Nous avons besoin d'une nouvelle sous-commande : bedtools closest

Attention, bien prendre note du message suivant :

L'output est à stocker dans un fichier nommé picVariant_sorted.bed

In [18]:
sort -k1,1 -k2,2n picVariant.bed > picVariant_sorted.bed

Nous tentons un essai en paramètre par défaut. Rediriger la sortie dans un fichier DEgenes_closest_picVariant.bed

In [19]:
bedtools closest -a picVariant_sorted.bed -b DEgenes.gff3 \
   > DEgenes_closest_picVariant.bed

Et si on prenait en compte le brin ?

In [20]:
bedtools closest -a picVariant_sorted.bed -b DEgenes.gff3 -s \
   > DEgenes_closest_picVariant.bed

Compter le nombre de lignes

In [21]:
wc -l DEgenes_closest_picVariant.bed

28080 DEgenes_closest_picVariant.bed

Bonus n°3¶

Devine la commande¶

Ce dernier bonus est l'occasion de voir des commandes qui peuvent être utiles mais que nous n'avons pas pu aborder précédemment

Awk¶

  1. Pour commencer, un peu de Awk. Vous trouverez un fichier nommé .DEgenes_all.gff3 (ls -la pour le voir). Voici comment récupérer uniquement les annotations de gènes :
In [ ]:
awk 'BEGIN {OFS=FS="\t"} $3 == "gene" {print}' ".DEgenes_all.gff3" \
   > DEgenes.gff3
  • FS et OFS : Permettent de définir le séparateur de champ du fichier
  • $3 == "gene" : Permet de vérifier si la colonne 3 est une annotation de gène
  • print : Imprime à l'écran la ligne quand la condition est vraie
  1. La dernière colonne du fichier intersect_flanck2kDEgenes_n_variant.txt contient le nombre de variants par promoteur. Comment s'assurer que nous avons bien pris en compte les 405 variants initialement trouvés en bonus 1 ?
In [25]:
awk -F '\t' 'BEGIN{s=0}{s+=$NF}END{print s}' \
   intersect_flanck2kDEgenes_n_variant.txt

405
  • -F : Permet de définir le séparateur de champ dans l'input
  • s=0 : Initialisation d'une variable s qui vaut 0 avant de lire le fichier
  • s+=$NF : Pour chaque ligne, on ajoute la valeur de la dernière colonne à la somme s
  • print s : Une fois l'ensemble du fichier parcouru, on imprime la dernière valeur de la somme s

Alternative à grep -v¶

Dans l'exercice, point D, on utilise grep -v pour extraire les coordonnées génomiques des gènes non DE.

Il était aussi possible de le faire avec la commande intersect.

Trouverez vous comment ?

Attention aux options !

In [23]:
bedtools intersect -a hsGRCh38.genes.gff3 -b DEgenes.gff3 \
   -v -f 1 -F 1 -s -wa > notDEgenes.gff3
wc -l notDEgenes.gff3

21492 notDEgenes.gff3

L'option la plus importante pour cette question est donc l'option -v.

Mais par défaut, intersect valide une intersection dès qu'il y a une base en commun entre les 2 régions. Dans ce cas d'usage où l'on traite des gènes en entier, il suffit que quelques gènes se chevauchent pour qu'ils soient exclus à tort. La façon correcte est donc de préciser que l'on veut valider les chevauchements que si les gènes se chevauchent à 100% entre eux. Les options pour spécifier un % de chavauchement (exprimés entre [0-1], 1 pour 100%) sont -F pour le fichier A et -f pour le fichier B. Il faut aussi que les 2 gènes soient sur le même brin (option -s).