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.

Données disponibles¶

  • Une liste de gènes différentiellement exprimés (Cellules avec un KO du gène beta-2-microglobulin, infectées par l'herpes simplex virus type 1 vs non infectées) (TXT)
  • Une liste des sites de fixation de H3K4me3 déterminés par Chip-seq (BED)
  • Une liste de variants (VCF)
  • Un fichier d'annotation humain (GFF3)
  • La liste des chromosomes humains et leur taille (LEN)

Avant de commencer¶

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

In [1]:
# Afficher le chemin
pwd

/shared/ifbstor1/projects/2422_ebaii_n1/cours_commun/TP_croisement

Accès aux données du TP¶

  1. Rendez vous dans le dossier data/ et listez en le contenu
In [2]:
# Se déplacer dans le dossier data
cd data
# Lister son contenu
ls

DEgenes.txt  chrs.len  hsGRCh38.genes.gff3  picChip.bed  variant.vcf

Ce TP repose sur 4 jeux de données, issues d'analyses précédentes ainsi que sur les annotations du génome humain :

  • la longueur des chrs (et des contigs) de la version GRCh27 du génome humain (chrs.len)
  • une liste de SNP variants (variant.vcf)
  • une liste des identifiants de gènes différentiellement exprimés (DEgenes.txt)
  • une liste des pics issus d'une analyse de ChIP-seq (picChip.bed), ici une marque d'histone h3k4me3
  • les annotations correspondant aux gènes uniquement du génome humain GRCh38 (hsGRCh38.genes.gff3)
  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]:
# Afficher les deux premières lignes de tous les fichiers
head -n 2 *

==> DEgenes.txt <==
ENSG00000004846
ENSG00000005981

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

==> 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]:
# Afficher les deux dernières lignes du fichier variant.vcf
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) et bc (1.07.1) dont nous aurons besoin plus loin
In [5]:
# Charger les outils
module load bedtools/2.30.0 bc/1.07.1

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

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

In [6]:
# Trois premières lignes des fichiers listés
head -n 3 DEgenes.txt picChip.bed hsGRCh38.genes.gff3

==> DEgenes.txt <==
ENSG00000004846
ENSG00000005981
ENSG00000006747

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

==> 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
1	ensembl_havana	gene	685679	686673	.	-	.	ID=gene:ENSG00000284662;Name=OR4F16;biotype=protein_coding;description=olfactory receptor family 4 subfamily F member 16 [Source:HGNC Symbol%3BAcc:HGNC:15079];gene_id=ENSG00000284662;logic_name=ensembl_havana_gene;version=1

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

Drawing

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

Drawing

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

Drawing

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

Drawing

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

Drawing

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 [7]:
# Rechercher les gènes contenu dans le fichier txt dans le fichier gff3
grep -f DEgenes.txt hsGRCh38.genes.gff3 > DEgenes.gff3

On peut vérifier ensuite que le nombre de gènes récupérés dans le fichier gff3 est celui attendu (nombre de gènes de notre liste). Sachant qu'une ligne = un gène, quelles commandes pouvez vous exécuter ?

In [8]:
# Compter le nombre de ligne de chacun des fichiers
wc -l DEgenes.txt
wc -l DEgenes.gff3
wc -l hsGRCh38.genes.gff3

162 DEgenes.txt
162 DEgenes.gff3
21492 hsGRCh38.genes.gff3

Il est aussi possible de visualiser les 5 premières lignes de notre document nouvellement créé pour voir ce qu'il contient

In [9]:
# Afficher les 5 premières lignes du fichier
head -n 5 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
1	ensembl_havana	gene	74198235	74544393	.	+	.	ID=gene:ENSG00000259030;Name=FPGT-TNNI3K;biotype=protein_coding;description=FPGT-TNNI3K readthrough [Source:HGNC Symbol%3BAcc:HGNC:42952];gene_id=ENSG00000259030;logic_name=ensembl_havana_gene;version=7
1	ensembl_havana	gene	111722064	111755824	.	-	.	ID=gene:ENSG00000284755;Name=AL049557.1;biotype=protein_coding;description=inka box actin regulator 2 [Source:NCBI gene%3BAcc:55924];gene_id=ENSG00000284755;logic_name=ensembl_havana_gene;version=1
1	ensembl_havana	gene	147756199	147773362	.	-	.	ID=gene:ENSG00000265107;Name=GJA5;biotype=protein_coding;description=gap junction protein alpha 5 [Source:HGNC Symbol%3BAcc:HGNC:4279];gene_id=ENSG00000265107;logic_name=ensembl_havana_gene;version=2

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(cliquez sur les sous-commandes pour voir leur aide)

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 [ ]:
# Intersection du fichier A et du fichier B
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 [11]:
# Comptage des lignes
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 du fichier d'annotation hsGRCh38.genes.gff3 dont le nom n'est pas dans la liste des gènes DE DEgenes.txt puis on compte combien de gènes sont concernés.

In [12]:
# Recherche inversée du contenu du fichier dans le gff3
grep -v -f DEgenes.txt hsGRCh38.genes.gff3 > notDEgenes.gff3
# Comptage des lignes
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 [13]:
# Nouvel intersect mais avec le fichier de gène non DE
bedtools intersect -a notDEgenes.gff3 -b picChip.bed \
-wa -u > intersect_notDEgenes_picChip.gff3
# Comptage des lignes
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 [14]:
# Comptage des lignes de tous les fichiers
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 deux dernières lignes des fichiers DEgenes.gff3, variant.vcf et chrs.len :

In [15]:
# Visualisation de la fin des fichiers
tail -n2 DEgenes.gff3 variant.vcf chrs.len

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

==> chrs.len <==
KI270392.1	971
KI270394.1	970
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 [16]:
# Recherche de régions flanquantes
bedtools flank -i DEgenes.gff3 -g chrs.len \
-l 2000 -r 0 -s > flank2kDEgenes.gff3
# Visualisation des trois premières lignes
head -n3 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
1	ensembl_havana	gene	74196235	74198234	.	+	.	ID=gene:ENSG00000259030;Name=FPGT-TNNI3K;biotype=protein_coding;description=FPGT-TNNI3K readthrough [Source:HGNC Symbol%3BAcc:HGNC:42952];gene_id=ENSG00000259030;logic_name=ensembl_havana_gene;version=7

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 5 premières lignes du fichier quand celui ci est généré

In [17]:
# Nouvelle intersection
bedtools intersect -a variant.vcf -b flank2kDEgenes.gff3 \
-wa -u > intersect_flank2kDEgenes_variant.txt
# Visualisation des 5 premières lignes
head -n 5 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
1	147773875	rs72546659	C	CT	.	.	RS=72546659;RSPOS=147773875;dbSNPBuildID=130;SSR=0;SAO=0;VP=0x0501100a000504003e000200;GENEINFO=GJA5:2702|LOC102723321:102723321;WGT=1;VC=DIV;TPA;SLO;INT;R5;ASP;VLD;KGPhase1;KGPhase3;CAF=0.9968,0.003195;COMMON=1;TOPMED=0.99143890163098878,0.00856109836901121
1	147775166	rs10641806	A	AAT	.	.	RS=10641806;RSPOS=147775166;dbSNPBuildID=135;SSR=0;SAO=0;VP=0x0500000a000515003e000200;GENEINFO=GJA5:2702|LOC102723321:102723321;WGT=1;VC=DIV;INT;R5;ASP;VLD;G5;KGPhase1;KGPhase3;CAF=0.05811,0.9419;COMMON=1
1	173918736	rs577103838	G	GGAGGCA	.	.	RS=577103838;RSPOS=173918736;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000020005040026000200;GENEINFO=SERPINC1:462;WGT=1;VC=DIV;R5;ASP;VLD;KGPhase3;CAF=0.9988,0.001198;COMMON=1;TOPMED=0.99892488532110091,0.00107511467889908

Combien avons-nous de promoteurs contenant des variants ?

In [18]:
# Comptage des lignes
wc -l intersect_flank2kDEgenes_variant.txt

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 ?

On utilise l'option -c qui permet de compter le nombre de croisements et de le stocker dans un fichier (ici intersect_flanck2kDEgenes_n_variant.txt). On visualise les 5 premières lignes.

In [19]:
# Nouvelle intersection avec comptage
bedtools intersect -a flank2kDEgenes.gff3 -b variant.vcf \
-c > intersect_flanck2kDEgenes_n_variant.txt
# Visualisation des 5 premières lignes
head -n 5 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
1	ensembl_havana	gene	74196235	74198234	.	+	.	ID=gene:ENSG00000259030;Name=FPGT-TNNI3K;biotype=protein_coding;description=FPGT-TNNI3K readthrough [Source:HGNC Symbol%3BAcc:HGNC:42952];gene_id=ENSG00000259030;logic_name=ensembl_havana_gene;version=7	1
1	ensembl_havana	gene	111755825	111757824	.	-	.	ID=gene:ENSG00000284755;Name=AL049557.1;biotype=protein_coding;description=inka box actin regulator 2 [Source:NCBI gene%3BAcc:55924];gene_id=ENSG00000284755;logic_name=ensembl_havana_gene;version=1	0
1	ensembl_havana	gene	147773363	147775362	.	-	.	ID=gene:ENSG00000265107;Name=GJA5;biotype=protein_coding;description=gap junction protein alpha 5 [Source:HGNC Symbol%3BAcc:HGNC:4279];gene_id=ENSG00000265107;logic_name=ensembl_havana_gene;version=2	2

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 (créé dans le problème 1)
  • variant.vcf
  • picChip.bed

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

In [20]:
# Visualisation des 2 dernières lignes
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 [21]:
# Nouvelle intersection
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

# Visualisation des 5 premières lignes
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 :

Drawing

Trions le fichier picVariant.bed.

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

In [22]:
# Tri selon les colonnes 1 et 2 (2 en ordre numérique)
sort -k1,1 -k2,2n picVariant.bed > picVariant_sorted.bed

Nous tentons un essai avec le paramètrage par défaut. Rediriger la sortie dans un fichier DEgenes_closest_picVariant.bed

In [23]:
# Recherche des gènes les plus proches des pics
bedtools closest -a picVariant_sorted.bed \
-b DEgenes.gff3 > DEgenes_closest_picVariant.bed

Combien avait-on de pic de Chip avec une mutation et combien obtient-on de gènes à proximité ?

In [24]:
# Comptage des lignes
wc -l picVariant_sorted.bed DEgenes_closest_picVariant.bed

  28080 picVariant_sorted.bed
  28088 DEgenes_closest_picVariant.bed
  56168 total

On obtient 8 gènes de plus que notre nombre initial de pic. En regardant de plus près les 2 premières lignes, on voit que le premier pic est aux alentours des bases 713,000 alors que le gène est vers 3,000,000...

In [25]:
# Visualisation des deux premières lignes
head -n2 DEgenes_closest_picVariant.bed

1	713056	713670	region_8999	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	760194	766425	region_9005	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

La pertinence de l'impact d'une protéine fixée à l'ADN a plus de 3 millions de paire de base d'un gène se pose. Il serait plus pertinent de rechercher les pics situés à proximité immédiate des gènes DE (entre 2000 bases avant ou après un gène par exemple).

Quelle option peut-on utiliser pour obtenir l'information de la distance entre le pic et le gène ? Stockons l'information dans un fichier nommé DEgenes_closest_picVariant_distance.bed

In [26]:
# Recherche des gènes les plus proches 
bedtools closest -a picVariant_sorted.bed \
-b DEgenes.gff3 -D b > DEgenes_closest_picVariant_distance.bed

L'utilisation de l'option -D b permet de donner la distance par rapport aux gènes qui sont les seuls à avoir une information sur le brin (les pics étant non orientés). Une valeur négative signifie que le pic est en amont du gène et en aval pour les valeurs positives.

A ce stade, nous avons l'information de la distance mais aucun filtre.

In [27]:
# Visualisation des deux premières lignes
head -n2 DEgenes_closest_picVariant_distance.bed

1	713056	713670	region_8999	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	-32367434
1	760194	766425	region_9005	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	-32314679

Il est alors possible de mettre une commande supplémentaire pour filter. On utilise le langage Awk

N.B. : Awk est un langage à part entière, puissant mais qui peut être complexe de prime abord. Nous n'avons pas le temps de vous en expliquer toute la syntaxe mais de très bon site l'on fait, notamment ici et là.

In [28]:
# Recherche des gènes les plus proches et dans la zone définie entre [-2000; +2000]
bedtools closest -a picVariant_sorted.bed \
-b DEgenes.gff3 -D b | \
awk 'BEGIN{FS="\t"}{if($14>=-2000 && $14<=2000){print $0}}' \
> DEgenes_closest_picVariant_borne2000.bed

Rapide aperçu de la commande :

  • awk '' : On définie que tout ce qui est entre guillements simples correspond à du langage awk et non plus du bash.
  • BEGIN{FS="\t"} : On définie le séparateur de nos colonnes, ici la tabulation "\t"
  • {if($14>=-2000 && $14<=2000){print $0}} : Pour chaque ligne (donc chaque gène), on vérifie (if) si la colonne 14 ($14), qui contient la distance, a une valeur comprise entre -2000 ($14>=-2000) et (&&) +2000 ($14<=2000) inclus. Si c'est le cas, on affiche (print) la ligne entière ($0).

Affichons maintenant les deux premières lignes :

In [29]:
# Afficher les deux premières lignes
head -n2 DEgenes_closest_picVariant_borne2000.bed

1	33112919	33119940	region_10192	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	0
1	70875092	70879897	region_11035	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	0

On rapproche à nouveau ces données de la réalité biologique. Il est surement pertinent de ne regarder que les pics présents en amont des gènes DE. Pour cela, on change simplement les filtres sur la distance dans awk.

A vous de jouer !

Rappel : Les bornes de la zone en amont se définissent ainsi [-2000;0[

On écrit dans un fichier DEgenes_closest_picVariant_amont_gene_DE.bed

In [30]:
# Recherche des gènes les plus proches et dans la zone définie entre [-2000; 0[
bedtools closest -a picVariant_sorted.bed \
-b DEgenes.gff3 -D b | \
awk 'BEGIN{FS="\t"}{if($14>=-2000 && $14<0){print $0}}' \
> DEgenes_closest_picVariant_amont_gene_DE.bed

Notre liste est devenu bien plus raisonnable, affichons entièrement ce fichier

In [31]:
# Affichage des fichiers
cat DEgenes_closest_picVariant_amont_gene_DE.bed

10	71202570	71212017	region_4870	10	ensembl_havana	gene	71212570	71302864	.	+	.	ID=gene:ENSG00000107731;Name=UNC5B;biotype=protein_coding;description=unc-5 netrin receptor B [Source:HGNC Symbol%3BAcc:HGNC:12568];gene_id=ENSG00000107731;logic_name=ensembl_havana_gene;version=12	-553
10	71210834	71211813	region_4871	10	ensembl_havana	gene	71212570	71302864	.	+	.	ID=gene:ENSG00000107731;Name=UNC5B;biotype=protein_coding;description=unc-5 netrin receptor B [Source:HGNC Symbol%3BAcc:HGNC:12568];gene_id=ENSG00000107731;logic_name=ensembl_havana_gene;version=12	-757
14	77500630	77515764	region_50749	14	ensembl_havana	gene	77474394	77498850	.	-	.	ID=gene:ENSG00000100593;Name=ISM2;biotype=protein_coding;description=isthmin 2 [Source:HGNC Symbol%3BAcc:HGNC:23176];gene_id=ENSG00000100593;logic_name=ensembl_havana_gene;version=17	-1781

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

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 !

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. Comme ci-dessous :

Drawing

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 chevauchement (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).

Drawing

In [32]:
# Intersection inversée
bedtools intersect -a hsGRCh38.genes.gff3 -b DEgenes.gff3 \
-v -f 1 -F 1 -s -wa > notDEgenes.gff3
# Comptage des lignes
wc -l notDEgenes.gff3

21330 notDEgenes.gff3