diff --git a/2023/ebaiin1/TP_croisement/TP_croisement_de_donnees.ipynb b/2023/ebaiin1/TP_croisement/TP_croisement_de_donnees.ipynb new file mode 100644 index 0000000..45e788d --- /dev/null +++ b/2023/ebaiin1/TP_croisement/TP_croisement_de_donnees.ipynb @@ -0,0 +1,1186 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "676a7a2f-6578-4b9d-a288-3c338e705f70", + "metadata": { + "tags": [] + }, + "source": [ + "# Croisement des données " + ] + }, + { + "cell_type": "markdown", + "id": "4c1b28ac-acda-466a-9d75-c1157907afc5", + "metadata": {}, + "source": [ + "##### *Claire Toffano-Nioche, I2BC, Paris-Saclay*\n", + "\n", + "##### *Pauline Francois, Anses, Lyon*\n", + "\n", + "##### *Selon les ressources de Matthias Zytnicki, inrae, Toulouse*\n", + "\n", + "##### *et toute la team Roscoff*" + ] + }, + { + "cell_type": "markdown", + "id": "232ac2ed-00aa-4652-9a70-9c944bf00f81", + "metadata": {}, + "source": [ + "## Objectifs du cours\n", + "\n", + "- Permettre de croiser des données génomiques de types différents\n", + "- Renforcer les notions vues précedemment\n", + "- Se familiariser avec bedtools" + ] + }, + { + "cell_type": "markdown", + "id": "c11fd4d0-4c89-4287-8cf4-e1af3cc6bfba", + "metadata": {}, + "source": [ + "Ce que ce cours n'est pas :\n", + "- Une réponse à une vraie question biologique\n", + "- Une méthode statistiquement valide" + ] + }, + { + "cell_type": "markdown", + "id": "3797b525-90a4-4126-8afe-817b8acfe8bc", + "metadata": { + "tags": [] + }, + "source": [ + "## Croisement de données\n", + "\n", + "#### Qu’est-ce c'est ?\n", + "- Une comparaison des positions ou intervalles génomiques: \n", + " - un variant par rapport à des gènes d’intérêt, \n", + " - des pics de ChIP-Seq par rapport à des gènes différentiellement exprimés, etc.\n", + " \n", + "#### À quel moment est-ce valide ?\n", + "- Lorsque vous cherchez des co-occurrences\n", + "- Lorsque vous donnez des distributions (de distance)\n", + "\n", + "#### À quel moment est-ce douteux ?\n", + "- Lorsque les résultats sont présentés comme significatifs.\n", + "\n", + "#### Données disponibles\n", + "- 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)\n", + "- Une liste des sites de fixation de H3K4me3 déterminés par Chip-seq (BED)\n", + "- Une liste de variants (VCF)\n", + "- Un fichier d'annotation humain (GFF3)\n", + "- La liste des chromosomes humains et leur taille (LEN)\n", + "\n", + "## Avant de commencer\n", + "\n", + "1. Vérifier que vous êtes dans votre espace projet :" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "id": "7a0bd48a-14bd-47d5-ade3-03ec8a10a06c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/shared/ifbstor1/projects/primir_orfevre_2/tp_croisement\n" + ] + } + ], + "source": [ + "pwd" + ] + }, + { + "cell_type": "markdown", + "id": "9a629a3f-7de5-46a2-8d73-4d2e5893c662", + "metadata": {}, + "source": [ + "2. Créer un dossier nommé tp_croisement et vous déplacer dedans" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "841914a7-436f-4c22-b098-38782b0ec298", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/shared/ifbstor1/projects/form_2022_32/TP_croisement/tp_croisement\n" + ] + } + ], + "source": [ + "mkdir -p tp_croisement\n", + "cd tp_croisement\n", + "pwd\n" + ] + }, + { + "cell_type": "markdown", + "id": "8bc23cff-130e-4e35-a7ae-7103ac98f04a", + "metadata": {}, + "source": [ + "## Accès aux données du TP" + ] + }, + { + "cell_type": "markdown", + "id": "ca66dbeb-82ea-4fe8-baf2-9a3bbddb62ae", + "metadata": {}, + "source": [ + "1. Regarder la liste des fichiers présents dans /shared/projects/form_2021_26/data/atelier_croisement/" + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "id": "5c7e5721-822f-4e1f-9802-dbc5244ec631", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "chrs.len\t\t differentially_expressed_genes.txt\n", + "common_all_20180418_div.vcf h3k4me3_k562.bed\n" + ] + } + ], + "source": [ + "ls /shared/projects/form_2022_32/data/atelier_croisement/" + ] + }, + { + "cell_type": "markdown", + "id": "c3d3e4f9-e9ad-45e1-b981-92a55e47d1af", + "metadata": {}, + "source": [ + "Ce TP repose sur 4 jeux de données, issues d'analayses précédentes ainsi que sur les annotations du génome humain : \n", + "\n", + "- la longueur des chrs (et des contigs) de la version GRCh27 du génome humain (*.len)\n", + "- une liste de SNP variants (*.vcf)\n", + "- une liste des identifiants de gènes différentiellement exprimés (*.txt)\n", + "- une liste des pics issus d'une analyse de ChIP-seq (*.bed), ici une marque d'histone\n", + "- les annotations du génome humain (*.gff3)\n", + "\n", + "2. Ces données sont \"loin\" (/shared/...). Pour simplifier l'écriture des lignes de commandes, on va utiliser la commande ln -s qui permet de créer des raccourcis. \n", + "Exécuter la liste de commandes suivantes et remplir les commentaires (lignes commençant par \"#\") :" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "710e33e6-a2bd-4249-8748-bc7103cdec84", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "total 20K\n", + "lrwxrwxrwx 1 cpatiou cpatiou 62 Nov 17 09:07 chrs.len -> /shared/projects/form_2022_32/data/atelier_croisement/chrs.len\n", + "lrwxrwxrwx 1 cpatiou cpatiou 88 Nov 17 09:07 DEgenes.txt -> /shared/projects/form_2022_32/data/atelier_croisement/differentially_expressed_genes.txt\n", + "lrwxrwxrwx 1 cpatiou cpatiou 65 Nov 17 09:07 hsGRCh38.gff3 -> /shared/bank/homo_sapiens/GRCh38/gff3/Homo_sapiens.GRCh38.94.gff3\n", + "lrwxrwxrwx 1 cpatiou cpatiou 70 Nov 17 09:07 picChip.bed -> /shared/projects/form_2022_32/data/atelier_croisement/h3k4me3_k562.bed\n", + "lrwxrwxrwx 1 cpatiou cpatiou 81 Nov 17 09:07 variant.vcf -> /shared/projects/form_2022_32/data/atelier_croisement/common_all_20180418_div.vcf\n" + ] + } + ], + "source": [ + "#lien gènes différentiellement exprimés\n", + "ln -s /shared/projects/form_2022_32/data/atelier_croisement/differentially_expressed_genes.txt DEgenes.txt\n", + "#lien fichier variants\n", + "ln -s /shared/projects/form_2022_32/data/atelier_croisement/common_all_20180418_div.vcf variant.vcf\n", + "#lien bed de pics Chip-seq\n", + "ln -s /shared/projects/form_2022_32/data/atelier_croisement/h3k4me3_k562.bed picChip.bed\n", + "#lien contigs \n", + "ln -s /shared/projects/form_2022_32/data/atelier_croisement/chrs.len .\n", + "#lien gff3\n", + "ln -s /shared/bank/homo_sapiens/GRCh38/gff3/Homo_sapiens.GRCh38.94.gff3 hsGRCh38.gff3\n", + "\n", + "ls -lh" + ] + }, + { + "cell_type": "markdown", + "id": "c69e3d9b-8e71-4e83-a0f2-ee6da9090549", + "metadata": {}, + "source": [ + "3. Afin de se familiariser avec le contenu des données, on affiche les 2 premières lignes de l'ensemble des fichiers :" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "id": "bf626e3c-fc81-4e24-b07c-4a92fe62e9f8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==> chrs.len <==\n", + "1\t248956422\n", + "10\t133797422\n", + "\n", + "==> DEgenes.txt <==\n", + "ENSG00000004846\n", + "ENSG00000005981\n", + "\n", + "==> hsGRCh38.gff3 <==\n", + "##gff-version 3\n", + "##sequence-region 1 1 248956422\n", + "\n", + "==> picChip.bed <==\n", + "22\t16192349\t16192565\tregion_1\n", + "22\t16846630\t16870710\tregion_2\n", + "\n", + "==> variant.vcf <==\n", + "##fileformat=VCFv4.0\n", + "##fileDate=20180418\n" + ] + } + ], + "source": [ + "head -n 2 *" + ] + }, + { + "cell_type": "markdown", + "id": "c0d60187-8a3e-4755-bc04-4cba30bafb27", + "metadata": {}, + "source": [ + "4. Pour deux des fichiers on ne voit que des commentaires (souvent des lignes qui commencent par \"#\") ; afficher les deux dernières lignes de ces fichiers :" + ] + }, + { + "cell_type": "code", + "execution_count": 57, + "id": "b330a69b-f820-4176-87d5-88d9028c5617", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==> hsGRCh38.gff3 <==\n", + "Y\thavana\texon\t56855244\t56855488\t.\t+\t.\tParent=transcript:ENST00000431853;Name=ENSE00001794473;constitutive=1;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001794473;rank=1;version=1\n", + "###\n", + "\n", + "==> variant.vcf <==\n", + "Y\t26614668\trs376925794\tTC\tT\t.\t.\tRS=376925794;RSPOS=26614671;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9968,0.003244;COMMON=1\n", + "Y\t26624486\trs771929773\tT\tTA\t.\t.\tRS=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\n" + ] + } + ], + "source": [ + "tail -n 2 hsGRCh38.gff3 variant.vcf" + ] + }, + { + "cell_type": "markdown", + "id": "8874a4ff-f99e-4582-9630-c0651bca04b6", + "metadata": {}, + "source": [ + "5. Afficher l'arborescence des fichiers avec la commande tree.\n", + "Observer les liens symboliques :" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "172f2484-d6a2-4fe6-ba8f-118a7c675db9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ".\n", + "├── chrs.len -> /shared/projects/form_2022_32/data/atelier_croisement/chrs.len\n", + "├── DEgenes.txt -> /shared/projects/form_2022_32/data/atelier_croisement/differentially_expressed_genes.txt\n", + "├── hsGRCh38.gff3 -> /shared/bank/homo_sapiens/GRCh38/gff3/Homo_sapiens.GRCh38.94.gff3\n", + "├── picChip.bed -> /shared/projects/form_2022_32/data/atelier_croisement/h3k4me3_k562.bed\n", + "└── variant.vcf -> /shared/projects/form_2022_32/data/atelier_croisement/common_all_20180418_div.vcf\n", + "\n", + "0 directories, 5 files\n" + ] + } + ], + "source": [ + "tree" + ] + }, + { + "cell_type": "markdown", + "id": "92acd376-3647-4e9e-89b2-00b5b4b860bf", + "metadata": {}, + "source": [ + "6. Activer les outils bedtools (v2.30.0) et bc (1.07.1)" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "id": "49e89bca-a3da-4b88-9699-666dd55ce5a8", + "metadata": {}, + "outputs": [], + "source": [ + "module load bedtools/2.30.0 bc/1.07.1" + ] + }, + { + "cell_type": "markdown", + "id": "107cabf7-3ed2-4fdd-8fa5-c1c24d1d43bb", + "metadata": { + "tags": [] + }, + "source": [ + "## Problème\n", + "\n", + "##### *Question scientifique*\n", + "\n", + "Mes gènes différentiellement exprimés contiennent-ils plus de pics de chip qu'attendu ? (oui/non, pas de quantification)\n", + "\n", + "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 ?\n", + "\n", + "##### *Données*\n", + "- DEgenes.txt\n", + "- picChip.bed\n", + "- hsGRCh38.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 60, + "id": "bb52b1d0-7643-40cd-9db3-263d568ec9f9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==> DEgenes.txt <==\n", + "ENSG00000004846\n", + "ENSG00000005981\n", + "ENSG00000006747\n", + "ENSG00000015568\n", + "ENSG00000047457\n", + "\n", + "==> picChip.bed <==\n", + "22\t16192349\t16192565\tregion_1\n", + "22\t16846630\t16870710\tregion_2\n", + "22\t17067019\t17067283\tregion_3\n", + "22\t17070620\t17114004\tregion_4\n", + "22\t17128021\t17136091\tregion_5\n", + "Y\t.\tbiological_region\t26627457\t26628186\t0.997\t+\t.\texternal_name=rank %3D 1;logic_name=firstef\n", + "Y\thavana\tpseudogene\t56855244\t56855488\t.\t+\t.\tID=gene:ENSG00000235857;Name=CTBP2P1;biotype=processed_pseudogene;description=C-terminal binding protein 2 pseudogene 1 [Source:HGNC Symbol%3BAcc:HGNC:23940];gene_id=ENSG00000235857;logic_name=havana;version=1\n", + "Y\thavana\tpseudogenic_transcript\t56855244\t56855488\t.\t+\t.\tID=transcript:ENST00000431853;Parent=gene:ENSG00000235857;Name=CTBP2P1-201;biotype=processed_pseudogene;tag=basic;transcript_id=ENST00000431853;transcript_support_level=NA;version=1\n", + "Y\thavana\texon\t56855244\t56855488\t.\t+\t.\tParent=transcript:ENST00000431853;Name=ENSE00001794473;constitutive=1;ensembl_end_phase=-1;ensembl_phase=-1;exon_id=ENSE00001794473;rank=1;version=1\n", + "###\n" + ] + } + ], + "source": [ + "head -n 5 DEgenes.txt picChip.bed\n", + "tail -n 5 hsGRCh38.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "95e6c422-9dc5-4e3e-9239-81c0e179f8b5", + "metadata": {}, + "source": [ + "DEgenes.txt est notre liste de gènes DE :\n", + "\n", + "![](images/liste_geneDE.png)\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "picChip.bed est notre liste de pics issus du chip-seq et positionnés sur le génome :\n", + "\n", + "![](images/liste_chip.png)\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "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) :\n", + "\n", + "![](images/genes_position_genome.png)\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "Le positionnement des données sur le génome permet maintenant de les croiser :\n", + "\n", + "![](images/associationGeneDE_chip.png)\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "Le but sera de compter le nombre gène dans chacune des 4 classes possibles :\n", + "\n", + "![](images/geneDE_Chip.png)" + ] + }, + { + "cell_type": "markdown", + "id": "baac3578-e77e-42e3-9775-95771c05bb0e", + "metadata": { + "tags": [] + }, + "source": [ + "##### *Protocole envisagé*\n", + "
\n", + "\n", + " A. Extraire les intervalles génomiques des gènes différentiellement exprimés (DE)\n", + " \n", + " B. Comparer les intervalles génomiques des gènes DE avec les régions H3K4me3 (ie. l'emplacement des pics de chip-seq)\n", + " \n", + " C. Compter le nombre de chevauchements entre les gènes DE et les pics\n", + "\n", + " D. Trouver les gènes non-différentiellement exprimés\n", + "\n", + " E. Comparer ces intervalles génomiques avec les régions H3K4me3\n", + "\n", + " F. Compter le nombre de chevauchements entre les gènes non DE et les pics\n", + "\n", + " G. Comparer les nombres de chevauchements\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "cf068431-71f8-463c-bdb1-9d97525ec1d8", + "metadata": { + "tags": [] + }, + "source": [ + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "#### A. Extraire les intervalles génomiques des gènes différentiellement exprimés (DE)\n", + "\n", + "1. Chercher une liste de mots dans un texte et rediriger la sortie standard dans un fichier nommé DEgenes_all.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "52d15cbf-191b-44dd-a14c-b0a1df960035", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "162 DEgenes.txt\n", + "1106 DEgenes_all.gff3\n", + "2814993 hsGRCh38.gff3\n", + "1\tensembl_havana\tgene\t33081104\t33120530\t.\t+\t.\tID=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\n", + "1\thavana\tlnc_RNA\t33081104\t33093400\t.\t+\t.\tID=transcript:ENST00000478635;Parent=gene:ENSG00000142920;Name=AZIN2-211;biotype=processed_transcript;transcript_id=ENST00000478635;transcript_support_level=3;version=5\n", + "1\tensembl_havana\tmRNA\t33081113\t33120394\t.\t+\t.\tID=transcript:ENST00000294517;Parent=gene:ENSG00000142920;Name=AZIN2-201;biotype=protein_coding;ccdsid=CCDS375.1;tag=basic;transcript_id=ENST00000294517;transcript_support_level=1;version=10\n", + "1\thavana\tlnc_RNA\t33081148\t33120529\t.\t+\t.\tID=transcript:ENST00000481886;Parent=gene:ENSG00000142920;Name=AZIN2-212;biotype=processed_transcript;transcript_id=ENST00000481886;transcript_support_level=2;version=5\n", + "1\thavana\tlnc_RNA\t33081150\t33120526\t.\t+\t.\tID=transcript:ENST00000475935;Parent=gene:ENSG00000142920;Name=AZIN2-208;biotype=processed_transcript;transcript_id=ENST00000475935;transcript_support_level=2;version=5\n", + "1\thavana\tlnc_RNA\t33081155\t33092221\t.\t+\t.\tID=transcript:ENST00000462920;Parent=gene:ENSG00000142920;Name=AZIN2-205;biotype=processed_transcript;transcript_id=ENST00000462920;transcript_support_level=5;version=5\n", + "1\thavana\tlnc_RNA\t33081155\t33120530\t.\t+\t.\tID=transcript:ENST00000471119;Parent=gene:ENSG00000142920;Name=AZIN2-206;biotype=processed_transcript;transcript_id=ENST00000471119;transcript_support_level=2;version=5\n", + "1\thavana\tlnc_RNA\t33081161\t33084127\t.\t+\t.\tID=transcript:ENST00000495135;Parent=gene:ENSG00000142920;Name=AZIN2-217;biotype=processed_transcript;transcript_id=ENST00000495135;transcript_support_level=5;version=5\n", + "1\thavana\tlnc_RNA\t33081163\t33096861\t.\t+\t.\tID=transcript:ENST00000477570;Parent=gene:ENSG00000142920;Name=AZIN2-209;biotype=processed_transcript;transcript_id=ENST00000477570;transcript_support_level=3;version=5\n", + "1\thavana\tlnc_RNA\t33081164\t33094713\t.\t+\t.\tID=transcript:ENST00000483027;Parent=gene:ENSG00000142920;Name=AZIN2-213;biotype=processed_transcript;transcript_id=ENST00000483027;transcript_support_level=5;version=5\n" + ] + } + ], + "source": [ + "grep -f DEgenes.txt hsGRCh38.gff3 > DEgenes_all.gff3\n", + "wc -l DEgenes.txt\n", + "wc -l DEgenes_all.gff3\n", + "wc -l hsGRCh38.gff3\n", + "head DEgenes_all.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "66b60524-98be-400e-8b88-2b428bf302f8", + "metadata": {}, + "source": [ + "2. Se restreindre aux gènes et rediriger la sortie standard dans un fichier nommé DEgenes.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 75, + "id": "e8629629-e562-401d-a95f-3253d06275d4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "gene\n", + "1\tensembl_havana\tgene\t33081104\t33120530\t.\t+\t.\tID=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\n", + "1\tensembl_havana\tgene\t70852353\t71047808\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t74198235\t74544393\t.\t+\t.\tID=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\n", + "1\tensembl_havana\tgene\t111722064\t111755824\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t147756199\t147773362\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t149886975\t149887364\t.\t+\t.\tID=gene:ENSG00000184260;Name=HIST2H2AC;biotype=protein_coding;description=histone cluster 2 H2A family member c [Source:HGNC Symbol%3BAcc:HGNC:4738];gene_id=ENSG00000184260;logic_name=ensembl_havana_gene;version=5\n", + "1\tensembl_havana\tgene\t173903804\t173917378\t.\t-\t.\tID=gene:ENSG00000117601;Name=SERPINC1;biotype=protein_coding;description=serpin family C member 1 [Source:HGNC Symbol%3BAcc:HGNC:775];gene_id=ENSG00000117601;logic_name=ensembl_havana_gene;version=13\n", + "1\tensembl_havana\tgene\t183186238\t183244900\t.\t+\t.\tID=gene:ENSG00000058085;Name=LAMC2;biotype=protein_coding;description=laminin subunit gamma 2 [Source:HGNC Symbol%3BAcc:HGNC:6493];gene_id=ENSG00000058085;logic_name=ensembl_havana_gene;version=14\n", + "1\tensembl_havana\tgene\t203026498\t203078740\t.\t+\t.\tID=gene:ENSG00000143847;Name=PPFIA4;biotype=protein_coding;description=PTPRF interacting protein alpha 4 [Source:HGNC Symbol%3BAcc:HGNC:9248];gene_id=ENSG00000143847;logic_name=ensembl_havana_gene;version=15\n", + "1\tensembl_havana\tgene\t205913048\t205943460\t.\t-\t.\tID=gene:ENSG00000174502;Name=SLC26A9;biotype=protein_coding;description=solute carrier family 26 member 9 [Source:HGNC Symbol%3BAcc:HGNC:14469];gene_id=ENSG00000174502;logic_name=ensembl_havana_gene;version=18\n" + ] + } + ], + "source": [ + "awk -F\"\\t\" 'BEGIN {OFS=FS} $3 == \"gene\" {print}' \"DEgenes_all.gff3\" > DEgenes.gff3 \n", + "cut -f 3 DEgenes.gff3 | sort | uniq\n", + "head DEgenes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "86d3b689-c558-47df-8689-058b16d8cad5", + "metadata": {}, + "source": [ + "Ou, en plus court, générer directement le fichier DEgenes.gff3 (à l'aide d'un pipe) :" + ] + }, + { + "cell_type": "code", + "execution_count": 74, + "id": "0e17bf37-263e-4664-8dd9-7ea923cdbfd7", + "metadata": {}, + "outputs": [], + "source": [ + "grep -f DEgenes.txt hsGRCh38.gff3 > DEgenes_all.gff3 | awk -F\"\\t\" 'BEGIN {OFS=FS} $3 == \"gene\" {print}' \"DEgenes_all.gff3\" > DEgenes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "1696cd01-46de-440f-9492-77c89ce369d5", + "metadata": { + "tags": [] + }, + "source": [ + ".\n", + "\n", + "\n", + "#### B. Comparer les intervalles génomiques des gènes DE avec les régions H3K4me3 (ie. l'emplacement des pics de chip-seq)\n", + "\n", + "Comme on va utiliser la suite d'outils que propose [Bedtools](https://bedtools.readthedocs.io/en/latest/content/bedtools-suite.html), avant tout, faisons un tour sur la documentation.\n", + "\n", + "Chaque sous-commande fonctionne sur le même principe :\n", + "\n", + "bedtools [sub-command] [options]\n", + "\n", + "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.\n", + "\n", + "Les trois qui nous intéresserons dans un premier temps sont [intersect](https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html), [flank](https://bedtools.readthedocs.io/en/latest/content/tools/flank.html) et [closest](https://bedtools.readthedocs.io/en/latest/content/tools/closest.html)\n", + "\n", + "Bedtools accepte quasiment tous les formats de fichiers tabulés (GFF3, BED, etc)\n", + "\n", + "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 .....\" ?\n", + "\n", + "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 ?\n", + "\n", + "Rediriger la sortie standard dans un fichier nommé intersect_DEgenes_picChip.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 88, + "id": "05adfa85-3fc9-41ef-8130-5fa279446f60", + "metadata": {}, + "outputs": [], + "source": [ + "bedtools intersect -a DEgenes.gff3 -b picChip.bed -wa -u > intersect_DEgenes_picChip.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "ffff467e-cdc8-4e25-9eb7-100b1e9feba0", + "metadata": {}, + "source": [ + ".\n", + "\n", + "#### C. Compter le nombre de chevauchements entre les gènes DE et les pics\n", + "\n", + "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." + ] + }, + { + "cell_type": "code", + "execution_count": 89, + "id": "c18ab42b-aed0-424d-81ec-be40be39ea9d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "77 intersect_DEgenes_picChip.gff3\n" + ] + } + ], + "source": [ + "wc -l intersect_DEgenes_picChip.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "b21d36e1-75b0-4026-b73a-116946280ac3", + "metadata": { + "tags": [] + }, + "source": [ + ".\n", + "\n", + "\n", + "#### D. Trouver des gènes non-différentiellement exprimés\n", + "\n", + "1. Sélectionner tous les gènes présents dans l'annotation et les écrire dans un fichier nommé all_genes.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 79, + "id": "fdbdef07-c991-41b8-ba13-407f09b732de", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "gene\n", + "21492 all_genes.gff3\n" + ] + } + ], + "source": [ + "awk -F\"\\t\" 'BEGIN {OFS=FS} $3 == \"gene\" {print}' \"hsGRCh38.gff3\" > all_genes.gff3 \n", + "cut -f 3 all_genes.gff3 | sort | uniq\n", + "wc -l all_genes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "df2d7966-f4f4-43fd-b8ea-eb5b9a1b5996", + "metadata": {}, + "source": [ + "2. Ne garder que les gènes qui ne sont pas dans le fichier des gènes DE (regarder l'option -v) dans un fichier notDEgenes.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 91, + "id": "b80b3af9-7cb5-43bd-b8d5-8bd12e9542f9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "21180 notDEgenes.gff3\n" + ] + } + ], + "source": [ + "bedtools intersect -a all_genes.gff3 -b DEgenes.gff3 -v > notDEgenes.gff3\n", + "wc -l notDEgenes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "e69a0772-f05e-47f1-af42-96a0902652be", + "metadata": {}, + "source": [ + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "L'option la plus importante pour cette question est donc l'option -v. Mais par défaut, intersect valide une interscetion 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 (c'est le cas pour 4 gènes dans ces données). 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). Du coup, une commande plus exacte est :" + ] + }, + { + "cell_type": "code", + "execution_count": 96, + "id": "eea1f27a-4a66-4a40-aeb7-ef48ac7c770e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "21330 notDEgenes.gff3\n" + ] + } + ], + "source": [ + "bedtools intersect -a all_genes.gff3 -b DEgenes.gff3 -v -f 1 -F 1 -s -wa > notDEgenes.gff3\n", + "wc -l notDEgenes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "31dd0b7a-fd24-4430-8f27-5670b78c57bd", + "metadata": {}, + "source": [ + "#### E. Comparer ces intervalles génomiques avec les régions H3K4me3\n", + "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" + ] + }, + { + "cell_type": "code", + "execution_count": 97, + "id": "54320590-5175-455b-bca4-c8bb7dde1823", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "9438 intersect_notDEgenes_picChip.gff3\n" + ] + } + ], + "source": [ + "bedtools intersect -a notDEgenes.gff3 -b picChip.bed -wa -u > intersect_notDEgenes_picChip.gff3\n", + "wc -l intersect_notDEgenes_picChip.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "cd6aeff0-ff25-4ea5-bb8e-ac0154b34e10", + "metadata": {}, + "source": [ + ".\n", + "\n", + "\n", + "#### F. Compter le nombre de chevauchements entre les gènes non DE et les pics\n", + "\n", + "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 :\n", + "- Nombre de gènes DE\n", + "- Nombre de gènes DE avec chevauchement\n", + "- Nombre de gènes non DE\n", + "- Nombre de gènes non DE avec chevauchement\n" + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "id": "1d212113-1ab5-416f-b735-a1be6cd51c96", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "162 DEgenes.gff3\n", + "77 intersect_DEgenes_picChip.gff3\n", + "21330 notDEgenes.gff3\n", + "9438 intersect_notDEgenes_picChip.gff3\n" + ] + } + ], + "source": [ + "wc -l DEgenes.gff3\n", + "wc -l intersect_DEgenes_picChip.gff3\n", + "wc -l notDEgenes.gff3\n", + "wc -l intersect_notDEgenes_picChip.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "2385cd90-6b6a-4613-8587-6f146a9cc21a", + "metadata": {}, + "source": [ + ".\n", + "\n", + "#### G. Comparer les nombres de chevauchements\n", + "\n", + " \n", + " Ratios | avec pics | sans pic | total\n", + "-----------|----------|----------|-------\n", + "gene DE |77 |162-77 |162\n", + "gene non DE|9438 |21330-9438|21330 \n", + "total | | |21492\n", + "\n", + "Sous Bash (commande bc), un calcul peut être réalisé de la façon suivante :\n", + "echo \"CALCUL\" | bc -l\n", + "\n", + "C'est à dire que nous affichons le calcul sur la sortie standard que nous envoyons (pipe) dans la commande bc avec l'option -l. Cette option permet de faire le calcul à l'aide de la bibliothèque mathématique standard." + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "id": "a5603134-6c2f-45bf-96de-92cf193b7d8d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "85\n", + "11892\n" + ] + } + ], + "source": [ + "echo \"162-77\" | bc -l\n", + "echo \"21330 - 9438\" | bc -l" + ] + }, + { + "cell_type": "markdown", + "id": "f9948f3f-ca6b-453c-ba49-3ecd1b145823", + "metadata": {}, + "source": [ + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "## Bonus n°1\n", + "\n", + "##### *Question scientifique*\n", + "\n", + "Quels sont les variants qui sont dans les promoteurs (2kb en amont du gène) des gènes différentiellement exprimés ?\n", + "\n", + "##### *Données*\n", + "- DEgenes.gff3 (généré précédemment)\n", + "- variant.vcf\n", + "- chrs.len" + ] + }, + { + "cell_type": "code", + "execution_count": 102, + "id": "f6382f41-47cc-4f7f-9204-963cb4b60c29", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "==> DEgenes.gff3 <==\n", + "9\tensembl_havana\tgene\t41131306\t41199261\t.\t-\t.\tID=gene:ENSG00000215126;Name=CBWD6;biotype=protein_coding;description=COBW domain containing 6 [Source:HGNC Symbol%3BAcc:HGNC:31978];gene_id=ENSG00000215126;logic_name=ensembl_havana_gene;version=10\n", + "9\tensembl_havana\tgene\t64369394\t64413142\t.\t+\t.\tID=gene:ENSG00000172014;Name=ANKRD20A4;biotype=protein_coding;description=ankyrin repeat domain 20 family member A4 [Source:HGNC Symbol%3BAcc:HGNC:31982];gene_id=ENSG00000172014;logic_name=ensembl_havana_gene;version=12\n", + "9\tensembl_havana\tgene\t83623049\t83644130\t.\t+\t.\tID=gene:ENSG00000148057;Name=IDNK;biotype=protein_coding;description=IDNK%2C gluconokinase [Source:HGNC Symbol%3BAcc:HGNC:31367];gene_id=ENSG00000148057;logic_name=ensembl_havana_gene;version=15\n", + "9\tensembl_havana\tgene\t100442271\t100577636\t.\t+\t.\tID=gene:ENSG00000251349;Name=MSANTD3-TMEFF1;biotype=protein_coding;description=MSANTD3-TMEFF1 readthrough [Source:HGNC Symbol%3BAcc:HGNC:38838];gene_id=ENSG00000251349;logic_name=ensembl_havana_gene;version=3\n", + "9\tensembl_havana\tgene\t114784635\t114806126\t.\t-\t.\tID=gene:ENSG00000181634;Name=TNFSF15;biotype=protein_coding;description=TNF superfamily member 15 [Source:HGNC Symbol%3BAcc:HGNC:11931];gene_id=ENSG00000181634;logic_name=ensembl_havana_gene;version=7\n", + "9\tensembl_havana\tgene\t127738317\t127778741\t.\t-\t.\tID=gene:ENSG00000095370;Name=SH2D3C;biotype=protein_coding;description=SH2 domain containing 3C [Source:HGNC Symbol%3BAcc:HGNC:16884];gene_id=ENSG00000095370;logic_name=ensembl_havana_gene;version=19\n", + "9\tensembl_havana\tgene\t128149071\t128153455\t.\t+\t.\tID=gene:ENSG00000148346;Name=LCN2;biotype=protein_coding;description=lipocalin 2 [Source:HGNC Symbol%3BAcc:HGNC:6526];gene_id=ENSG00000148346;logic_name=ensembl_havana_gene;version=11\n", + "9\tensembl_havana\tgene\t132725578\t132878777\t.\t-\t.\tID=gene:ENSG00000165695;Name=AK8;biotype=protein_coding;description=adenylate kinase 8 [Source:HGNC Symbol%3BAcc:HGNC:26526];gene_id=ENSG00000165695;logic_name=ensembl_havana_gene;version=9\n", + "9\thavana\tgene\t133098121\t133163914\t.\t-\t.\tID=gene:ENSG00000285245;Name=AL162417.1;biotype=protein_coding;description=novel protein;gene_id=ENSG00000285245;logic_name=havana;version=1\n", + "X\tensembl_havana\tgene\t71107404\t71112108\t.\t-\t.\tID=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\n", + "\n", + "==> variant.vcf <==\n", + "Y\t26568803\trs777181523\tAC\tA\t.\t.\tRS=777181523;RSPOS=26568804;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9968,0.003244;COMMON=1\n", + "Y\t26577320\trs776633239\tATTGT\tA\t.\t.\tRS=776633239;RSPOS=26577321;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9984,0.001622;COMMON=1\n", + "Y\t26578605\trs371511160\tATAGT\tA\t.\t.\tRS=371511160;RSPOS=26578606;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000000005150026000200;WGT=1;VC=DIV;ASP;VLD;G5;KGPhase3;CAF=0.9384,0.06164;COMMON=1\n", + "Y\t26579456\trs768223745\tCA\tC\t.\t.\tRS=768223745;RSPOS=26579457;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9976,0.002433;COMMON=1\n", + "Y\t26588409\trs756738440\tTA\tT\t.\t.\tRS=756738440;RSPOS=26588410;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9976,0.002433;COMMON=1\n", + "Y\t26608388\trs761263486\tT\tTA\t.\t.\tRS=761263486;RSPOS=26608388;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005140026000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9862,0.01379;COMMON=1\n", + "Y\t26608388\trs77074025\tT\tTA\t.\t.\tRS=77074025;RSPOS=26608394;dbSNPBuildID=131;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9862,0.01379;COMMON=1\n", + "Y\t26614668\trs760366233\tTC\tT\t.\t.\tRS=760366233;RSPOS=26614669;dbSNPBuildID=144;SSR=0;SAO=0;VP=0x050000000005040024000200;WGT=1;VC=DIV;ASP;VLD;KGPhase3;CAF=0.9968,0.003244;COMMON=1\n", + "Y\t26614668\trs376925794\tTC\tT\t.\t.\tRS=376925794;RSPOS=26614671;dbSNPBuildID=138;SSR=0;SAO=0;VP=0x050000000005000002000200;WGT=1;VC=DIV;ASP;CAF=0.9968,0.003244;COMMON=1\n", + "Y\t26624486\trs771929773\tT\tTA\t.\t.\tRS=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\n" + ] + } + ], + "source": [ + "tail DEgenes.gff3 variant.vcf" + ] + }, + { + "cell_type": "markdown", + "id": "0fb529e4-73e4-4f52-b291-0acc3ea934c2", + "metadata": {}, + "source": [ + "##### *Protocole envisagé*\n", + "
\n", + "\n", + " A. Extraire la région située à 2kb en amont des gènes différentiellement exprimés (DE)\n", + " \n", + " B. Trouver tous les variants chevauchant les régions trouvées précedemment\n", + "\n", + "
\n", + "\n", + ".\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "fe524b39-abbb-4241-aaed-567bfea84066", + "metadata": {}, + "source": [ + "Nous allons avoir besoin d'une autre sous-commande de bedtools : [bedtools flank](https://bedtools.readthedocs.io/en/latest/content/tools/flank.html) Regardez bien ses options.\n", + "\n", + "Notons qu'il y a une option obligatoire, -g pour communiquer à l'outil la taille des chromosomes. Pourquoi ?\n", + "\n", + "Lancer la commande et écrire la sortie standard dans un fichier flank2kDEgenes.gff3" + ] + }, + { + "cell_type": "code", + "execution_count": 108, + "id": "be749754-969d-47ff-b393-a45076292ce5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1\tensembl_havana\tgene\t33079104\t33081103\t.\t+\t.\tID=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\n", + "1\tensembl_havana\tgene\t70850353\t70852352\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t74196235\t74198234\t.\t+\t.\tID=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\n", + "1\tensembl_havana\tgene\t111720064\t111722063\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t147754199\t147756198\t.\t-\t.\tID=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\n", + "1\tensembl_havana\tgene\t149884975\t149886974\t.\t+\t.\tID=gene:ENSG00000184260;Name=HIST2H2AC;biotype=protein_coding;description=histone cluster 2 H2A family member c [Source:HGNC Symbol%3BAcc:HGNC:4738];gene_id=ENSG00000184260;logic_name=ensembl_havana_gene;version=5\n", + "1\tensembl_havana\tgene\t173901804\t173903803\t.\t-\t.\tID=gene:ENSG00000117601;Name=SERPINC1;biotype=protein_coding;description=serpin family C member 1 [Source:HGNC Symbol%3BAcc:HGNC:775];gene_id=ENSG00000117601;logic_name=ensembl_havana_gene;version=13\n", + "1\tensembl_havana\tgene\t183184238\t183186237\t.\t+\t.\tID=gene:ENSG00000058085;Name=LAMC2;biotype=protein_coding;description=laminin subunit gamma 2 [Source:HGNC Symbol%3BAcc:HGNC:6493];gene_id=ENSG00000058085;logic_name=ensembl_havana_gene;version=14\n", + "1\tensembl_havana\tgene\t203024498\t203026497\t.\t+\t.\tID=gene:ENSG00000143847;Name=PPFIA4;biotype=protein_coding;description=PTPRF interacting protein alpha 4 [Source:HGNC Symbol%3BAcc:HGNC:9248];gene_id=ENSG00000143847;logic_name=ensembl_havana_gene;version=15\n", + "1\tensembl_havana\tgene\t205911048\t205913047\t.\t-\t.\tID=gene:ENSG00000174502;Name=SLC26A9;biotype=protein_coding;description=solute carrier family 26 member 9 [Source:HGNC Symbol%3BAcc:HGNC:14469];gene_id=ENSG00000174502;logic_name=ensembl_havana_gene;version=18\n" + ] + } + ], + "source": [ + "bedtools flank -i DEgenes.gff3 -g chrs.len -l 2000 -r 0 > flank2kDEgenes.gff3\n", + "head flank2kDEgenes.gff3" + ] + }, + { + "cell_type": "markdown", + "id": "e946f02b-c74c-4ba2-bd70-58d7a7df7a60", + "metadata": {}, + "source": [ + "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", + "\n", + "N'oubliez pas de regarder les options qui nous intéressent [ici](https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html) et visualisez les 5 premières lignes du fichier quand celui ci est généré" + ] + }, + { + "cell_type": "code", + "execution_count": 110, + "id": "5b9268f8-637c-4a57-91ba-ccf48b2e6419", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1\t33079702\trs10712676\tGA\tG\t.\t.\tRS=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\n", + "1\t70850638\trs545762042\tAT\tA\t.\t.\tRS=545762042;RSPOS=70850639;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000080005170026000200;GENEINFO=LOC102724572:102724572;WGT=1;VC=DIV;INT;ASP;VLD;G5A;G5;KGPhase3;CAF=0.7636,0.2364;COMMON=1;TOPMED=0.69237385321100917,0.30762614678899082\n", + "1\t70850755\trs55931738\tTC\tT\t.\t.\tRS=55931738;RSPOS=70850756;dbSNPBuildID=137;SSR=0;SAO=0;VP=0x05000008000515003e000200;GENEINFO=LOC102724572:102724572;WGT=1;VC=DIV;INT;ASP;VLD;G5;KGPhase1;KGPhase3;CAF=0.752,0.248;COMMON=1;TOPMED=0.80174885321100917,0.19825114678899082\n", + "1\t70850755\trs796218873\tTC\tT\t.\t.\tRS=796218873;RSPOS=70850758;dbSNPBuildID=146;SSR=0;SAO=0;VP=0x050000080005000002000200;GENEINFO=LOC102724572:102724572;WGT=1;VC=DIV;INT;ASP;CAF=0.752,0.248;COMMON=1\n", + "1\t74197393\trs201749320\tTTC\tT\t.\t.\tRS=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\n", + "396 intersect_flank2kDEgenes_variant.txt\n" + ] + } + ], + "source": [ + "bedtools intersect -a variant.vcf -b flank2kDEgenes.gff3 -wa -u > intersect_flank2kDEgenes_variant.txt\n", + "head -n 5 intersect_flank2kDEgenes_variant.txt\n", + "wc -l intersect_flank2kDEgenes_variant.txt" + ] + }, + { + "cell_type": "markdown", + "id": "c68bcce3-1e44-4502-b93c-2b60636f476f", + "metadata": {}, + "source": [ + "Pour aller plus loin avec bedtools intersect, on peut se demander combien y a-t-il de SNP par promoteur ? \n", + "\n", + "Stocker l'information dans un fichier intersect_flanck2kDEgenes_n_variant.txt et visualiser les 5 premières lignes" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4050259b-607e-4d70-9b46-ccd6044399e9", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "6473abf6-e41a-42c6-bdca-769091a32ed9", + "metadata": {}, + "source": [ + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "## Bonus n°2\n", + "\n", + "##### *Question scientifique*\n", + "\n", + "Quels sont les gènes différentiellement exprimés les plus proches des pics ChIP ET qui contiennent une mutation ?\n", + "\n", + "##### *Données*\n", + "- DEgenes.gff3 (généré dans le problème 1)\n", + "- variant.vcf\n", + "- picChip.bed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29e103eb-5651-4e87-b581-ac1ad7967ed2", + "metadata": {}, + "outputs": [], + "source": [ + "tail DEgenes.gff3 variant.vcf picChip.bed" + ] + }, + { + "cell_type": "markdown", + "id": "73162ff3-57f3-4e84-b31b-df7e3b10cea3", + "metadata": {}, + "source": [ + "##### *Protocole envisagé*\n", + "
\n", + "\n", + " 1. Trouver les pics qui contiennent une mutation\n", + " \n", + " 2. Trouver le gène le plus proche de chaque région précédemment trouvée\n", + "\n", + "
\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + ".\n", + "\n", + "#### Trouver les pics qui contiennent une mutation\n", + "\n", + "Rien de neuf sous le soleil, nous stockons l'output dans un fichier picVariant.bed." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34797d92-7e3c-4c8d-b3ee-b7c5443c6e9e", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "c21c87b8-5e9a-49e0-8832-f3bfb767b524", + "metadata": {}, + "source": [ + ".\n", + "\n", + "#### Trouver le gène le plus proche de chaque région précédemment trouvée\n", + "\n", + "Nous avons besoin d'une nouvelle sous-commande : [bedtools closest](https://bedtools.readthedocs.io/en/latest/content/tools/closest.html)\n", + "\n", + "Attention, bien prendre note du message suivant :\n", + "![](images/exo3_graph1.png)\n", + "\n", + "L'output est à stocker dans un fichier nommé picVariant_sorted.bed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e374d7be-9faa-45d9-8d36-edff635648d5", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "402d7339-b2da-45c3-993d-7be80bd5d036", + "metadata": {}, + "source": [ + "Nous tentons un essai en paramètre par défaut. Rediriger la sortie dans un fichier DEgenes_closest_picVariant.bed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc3313c2-0e89-40ef-b14a-3bdefd71897b", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "21505fce-e754-4466-96fa-72a687d20d3c", + "metadata": {}, + "source": [ + "Et si on prenait en compte le brin ?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ab6e978-4f88-4d09-b5c4-6ba55e89e2e0", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "42148205-cee4-455a-b481-666d2028fb09", + "metadata": {}, + "source": [ + "Afficher le nombre de ligne" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9c68d092-019c-4f06-b6c7-afff61222b81", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Bash", + "language": "bash", + "name": "bash" + }, + "language_info": { + "codemirror_mode": "shell", + "file_extension": ".sh", + "mimetype": "text/x-sh", + "name": "bash" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/2023/ebaiin1/TP_croisement/images/associationGeneDE_chip.png b/2023/ebaiin1/TP_croisement/images/associationGeneDE_chip.png new file mode 100644 index 0000000..f6925d3 Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/associationGeneDE_chip.png differ diff --git a/2023/ebaiin1/TP_croisement/images/exo3_graph1.png b/2023/ebaiin1/TP_croisement/images/exo3_graph1.png new file mode 100644 index 0000000..e84515b Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/exo3_graph1.png differ diff --git a/2023/ebaiin1/TP_croisement/images/geneDE_Chip.png b/2023/ebaiin1/TP_croisement/images/geneDE_Chip.png new file mode 100644 index 0000000..ff8cf9f Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/geneDE_Chip.png differ diff --git a/2023/ebaiin1/TP_croisement/images/genes_liste.png b/2023/ebaiin1/TP_croisement/images/genes_liste.png new file mode 100644 index 0000000..f8d94d1 Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/genes_liste.png differ diff --git a/2023/ebaiin1/TP_croisement/images/genes_position_genome.png b/2023/ebaiin1/TP_croisement/images/genes_position_genome.png new file mode 100644 index 0000000..d9eb786 Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/genes_position_genome.png differ diff --git a/2023/ebaiin1/TP_croisement/images/liste_chip.png b/2023/ebaiin1/TP_croisement/images/liste_chip.png new file mode 100644 index 0000000..cc3989a Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/liste_chip.png differ diff --git a/2023/ebaiin1/TP_croisement/images/liste_geneDE.png b/2023/ebaiin1/TP_croisement/images/liste_geneDE.png new file mode 100644 index 0000000..e6bb2cb Binary files /dev/null and b/2023/ebaiin1/TP_croisement/images/liste_geneDE.png differ