[PYTHON] Analyse ChIP-seq à partir de zéro
introduction
ChIP-seq (immunoprécipitation de la chromatine suivie d'un séquençage) est une mesure complète de l'endroit et de la fréquence à laquelle des liaisons de facteurs de transcription spécifiques et des modifications d'histones se produisent dans le génome. Vous pourrez partir de la construction de l'environnement, analyser en utilisant les données du papier, et enfin voir les résultats sur le navigateur du génome et faire des appels de pointe ~ Commençons par la construction de l'environnement immédiatement!
Environnement
Installation de miniconda3
Tout d'abord, installez le gestionnaire de paquets, miniconda. Le gestionnaire de paquets facilite l'installation et la gestion des outils lors de la construction de l'environnement.
Pour Linux
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
Pour Mac
wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh
bash Miniconda3-latest-MacOSX-x86_64.sh
Une fois que cela est fait, suivez les instructions et appuyez sur ENTRÉE ou tapez «oui». Vous pouvez taper «oui» pour toutes les questions. Lorsque vous avez terminé, fermez le terminal une fois. Après l'avoir démarré, nous ajouterons un canal à miniconda, veuillez donc procéder comme suit: __ Assurez-vous de le faire dans cet ordre __.
conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda
Si on vous dit que vous n'avez pas de conda ici, consultez l'Aide 1 ci-dessous.
Installation des outils nécessaires
Maintenant, utilisez miniconda pour installer les outils nécessaires. Cette fois, installez les outils suivants.
- SRA Toolkit Outils utilisés pour télécharger les données de séquence depuis SRA
- Trimmomatic Outil de coupe de l'adaptateur
- fastQC Un outil de contrôle qualité des fichiers fastq
- Bowtie2 Outil de cartographie
- Picard Outil pour supprimer les doublons PCR
- samtools Outils utilisés à divers endroits
- deepTools Outils utilisés pour corriger les données
- HOMER Divers outils
Installez en utilisant la commande conda de miniconda comme suit.
conda install sra-tools
conda install trimmomatic
conda install fastqc
conda install bowtie2
conda install picard
conda install samtools
conda install deeptools
conda install homer
Répondez y à tous Continuer ([y] / n)?.
L'environnement est maintenant prêt! Je vais aborder les données à partir de maintenant, mais pour éviter toute confusion, le ** suivant est écrit en supposant que toutes les commandes sont exécutées dans le même répertoire **. Allons chercher les données ~.
Télécharger les données de séquence
Vérifier les données
Téléchargez les données de séquence depuis SRA (Sequence Read Archive). Cette fois [Kagey MH * et al., *](Https: // Données ChIP-seq pour Med1 de cellules ES de souris sur www.nature.com/articles/nature09380) (GSM560348) et Creyghton * et al., * Nous utilisons également les données ChIP-seq pour H3K27Ac des cellules ES de souris de /50/21931.long) (GSM594579). Nous obtenons également les données de séquence d'entrée (GSM560357) du même papier que Med1. () Les numéros indiqués sont les numéros d'accès GEO pour chaque donnée. Le Gene Expression Omnibus (GEO) de NCBI est le numéro nécessaire pour parcourir ces données.
Avant de télécharger les données, vérifions d'abord l'existence de ces données sur le navigateur. Prenons les données de Med1 comme exemple. Premièrement, [site GEO](https: //www.ncbi Ouvrez .nlm.nih.gov / geo /) dans votre navigateur. Tapez «GSM560348» dans la fenêtre de recherche entourée en rouge dans l'image ci-dessous.
Veuillez examiner de près les différentes informations sur ces données, comme indiqué dans l'image ci-dessous.
Maintenant, cliquez sur le numéro à droite de l'endroit où SRA est écrit en bas de cette page.
Ensuite, vous serez redirigé vers la page suivante Le numéro encadré en rouge dans l'image est appelé le numéro SRR, qui est requis lors du téléchargement avec le numéro d'accès de ces données.
Ces chiffres sont écrits quelque part dans le papier afin que vous puissiez voir où se trouvent les données.
Télécharger les données de SRA
Utilisez la commande sratoolkit fastq-dump. Elle est très facile à utiliser et si les données sont asymétriques
fastq-vider le numéro SRR des données que vous souhaitez télécharger
Pour les extrémités jumelées
fastq-vider le numéro SRR des données que vous souhaitez télécharger--split-files
Cela devrait télécharger le fichier fastq, qui est la sortie des données brutes du séquenceur, dans le répertoire où vous avez exécuté cette commande.
Téléchargeons ces données.
fastq-dump SRR058988 #Med1
fastq-dump SRR066767 #H3K27Ac
fastq-dump SRR058997 #contribution
Vous pouvez également les écrire comme suit et les télécharger tous en même temps.
fastq-dump SRR058988 SRR066767 SRR058997
Ce processus télécharge le fichier .sra au format de fichier compressé spécifique à SRA et le convertit en fichier fastq. Cela prendra un certain temps, alors soyez patient. Utilisez ce temps pour ma bioinformatique Que diriez-vous des commentaires pour le grand public? ??
https://laborify.net/2019/11/30/michida-bioinformatics/
Après le téléchargement, le nom du fichier est donné par le numéro SRR, alors renommez-le pour le rendre plus facile à comprendre.
mv SRR058988.fastq Med1.fastq
mv SRR066767.fastq H3K27Ac.fastq
mv SRR058997.fastq Input.fastq
Coupe de l'adaptateur
Nettoyons maintenant les résultats de la séquence à l'aide de Trimmomatic.
trimmomatic SE -threads 4 -phred33 Med1.fastq Med1_trimmed.fastq ILLUMINACLIP:$HOME/miniconda3/share/trimmomatic/adapters/TruSeq3-SE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
Bien qu'il soit long, il n'est pas interrompu pour éviter les accidents dus à des ruptures. Immédiatement après trimmomatic, spécifiez si les données à analyser sont asymétriques (SE) ou par paires (PE). Dans les fils suivants` Spécifiez le nombre de fils à utiliser. «-Phred33» est un sort. Veillez à le saisir. Ensuite, saisissez le nom du fichier à découper et le nom du fichier après le découpage.
L'emplacement après ILLUMINACLIP est l'emplacement des informations de séquence de l'adaptateur, qui doivent se trouver sous le répertoire miniconda3. Réécrivez-le pour qu'il corresponde à l'emplacement de votre miniconda3 (vous n'avez rien à faire lors de l'installation. Il devrait être dans votre répertoire personnel comme ceci.) De plus, 2:30:10 représente respectivement le nombre de discordances autorisées, le seuil de clip palindrome et le seuil de clip simple. Fondamentalement, je pense que vous n'avez pas à vous soucier de cela. De plus, «LEADING: 3» et «TRAILING: 3» signifient supprimer les bases avec un score de qualité inférieur à 3 respectivement au début et à la fin de la lecture. «SLIDING WINDOW: 4: 15» signifie tous les 4 pb Regardez le score de qualité moyen et supprimez les parties qui sont inférieures à 15; et le dernier «MINLEN: 36» signifie supprimer de l'analyse celles dont les longueurs de dérivation sont inférieures à 36. J'ai utilisé les paramètres "Démarrage rapide" de la page Trimmomatic (http://www.usadellab.org/cms/?page=trimmomatic). Une fois terminé, un fichier appelé Med1_trimmed.fastq sera généré. Exécutez les deux autres données avec les mêmes options.
QC du fichier fastq après le rognage
Utilisez fastQC pour contrôler la qualité du fichier fastq découpé.
fastqc --threads 4 --nogroup -o . Med1_trimmed.fastq
Écrivez le nombre de threads avec --threads immédiatement après fastqc. Si vous écrivez le prochain --nogroup, la lecture à l'extrémité 3 sera également analysée. Renvoyez le résultat à -o Écrivez le répertoire à faire, puis écrivez le nom du fichier à faire.
Un fichier appelé Med1_trimmed_fastqc.html sera créé dans le répertoire spécifié dans -o, alors ouvrons-le avec un navigateur.
Si le résumé de gauche est presque vert, la qualité est bonne. Ces données sont trop propres ... Expliquez chaque index à https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html Merci beaucoup pour votre soutien continu à ce site, la bioinformatique.
Veuillez faire la même option pour les deux données restantes.
cartographie
Enfin, nous mapperons en utilisant Bowtie2! Et avant cela, nous devons construire l'index du génome, c'est-à-dire préparer le génome de référence nécessaire à la cartographie.
Téléchargez le tableau mm10 entier depuis la page UCSC avec la commande wget. J'ai décidé de créer un dossier appelé ref_genome et de le déposer là-bas. Masu.
mkdir ref_genome
cd ref_genome
wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz
tar -xzf chromFa.tar.gz
De plus, cette fois, nous éliminerons les séquences aléatoires et inconnues et le génome mitochondrial (chrM.fa).
rm *random.fa
rm chrUn*
rm chrM.fa
Utilisez cat pour transformer le fichier restant en un seul fichier appelé mm10.fa.
cat *.fa > mm10.fa
Je vais créer un répertoire appelé mm10_index dans la même hiérarchie que ref_genome et y enregistrer l'index.
cd .. #Maintenant réf_Si vous êtes dans le génome
mkdir mm10_index
bowtie2-build -f ./ref_genome/mm10.fa ./mm10_index/mm10_index --threads 4
bowtie2-build est la commande pour indexer Bowtie2. Écrivez le chemin vers l'emplacement de la séquence du génome dans l'option -f, puis le chemin vers l'index. Spécifiez le nombre de threads avec --threads Vous pouvez. Dans ce cas, vous devriez avoir 6 fichiers nommés mm10_index * .bt2 dans le répertoire mm10_index que vous avez créé précédemment. Vous ne devez le faire qu'une seule fois.
Cela prend beaucoup de temps. Utilisez ce temps pour étudier le langage R, qui est souvent utilisé pour l'analyse statistique en bioinformatique. (Merci!)
https://qiita.com/roadricefield/items/001c882f84dd093f4407
Je n'utiliserai pas R dans cet article, mais ...
- Cela prend vraiment du temps. C'est un bon niveau pour frapper la commande juste avant de rentrer à la maison et continuer si tu viens demain. Je vais dormir en regardant l'animation jusqu'ici aujourd'hui.
.........................
Bon matin tout le monde! Cela m'a pris 7 heures car j'ai oublié de spécifier --threads! Cartographions.
bowtie2 -p 4 -x ./mm10_index/mm10_index -U Med1_trimmed.fastq -S Med1.sam
-p est le nombre de threads, -x est l'index, -U est le fichier fastq à mapper, -S est le nom du fichier de sortie. C'est environ 30 minutes.
Lorsque vous avez terminé, utilisez samtools pour convertir le fichier sam en fichier bam.
samtools view -b -o Med1.bam Med1.sam
Ceci termine la cartographie! Veuillez faire la même option pour les deux données restantes.
Suppression des doublons PCR (facultatif)
Enfin, utilisez Picard pour supprimer les doublons de PCR, pas nécessairement.
samtools sort Med1.bam > Med1_sorted.bam #Vous devez trier le fichier bam pour utiliser Picard.
picard MarkDuplicates I=Med1_sorted.bam O=Med1_rm_dups.bam M=Med1_report.txt REMOVE_DUPLICATES=true
Le nom du fichier bam pour lequel vous voulez supprimer les doublons PCR dans ʻI, le nom du fichier après avoir supprimé les doublons dans ʻO, et écrivez le nom dans M car cela créera un rapport qui résume les résultats du calcul. Masu.
Faites les deux autres données avec les mêmes options.
Le répertoire sur lequel je travaille est de plus en plus encombré, je vais donc organiser les données ici. Les données Med1 ChIP-seq se trouvent dans le répertoire Med1_data et les données H3K27Ac ChIP-seq se trouvent dans le répertoire H3K27Ac_data. Déplacez les données vers un répertoire appelé ʻInput_data`.
mkdir Med1_data
mv Med1* Med1_data
mkdir H3K27Ac_data
mv H3K27Ac* H3K27Ac_data
mkdir Input_data
mv Input* Input_data
Observez les résultats avec un navigateur génomique
Créer un fichier bigWig
Commençons par convertir le fichier bam dans un format appelé bigWig qui est facile à voir dans le navigateur de génome. Pour cela, utilisez bamCoverage de deepTools. Dans le format bigWig, la valeur de correction du nombre de lectures dans la largeur de bac spécifiée est calculée pour tout le génome. Autrement dit, la puissance du signal de ChIP-seq pour chaque bin du génome est calculée. Pour ce faire, vous devez d'abord créer un fichier bam.bai, qui est un fichier d'index bam, utilisez donc samtools. Faisons-le.
samtools index Med1_data/Med1_rm_dups.bam
Cela créera un fichier d'index appelé Med1_rm_dups.bam.bai.
Maintenant, lancez bamCoverage. Assurez-vous de mettre le fichier bam et son fichier bam.bai dans le même répertoire et exécutez-le.
bamCoverage -b Med1_data/Med1_rm_dups.bam -p 4 --normalizeUsing RPGC --effectiveGenomeSize 2652783500 --binSize 1 -o Med1_data/Med1.bigwig
Écrivez le nom du fichier bam à convertir en fichier bigWig dans -b. -p est le nombre de threads. --NormalizeUsing sélectionne le type de valeur de correction à calculer dans chaque bac. RPKM, Vous pouvez sélectionner CPM, BPM, RPGC, Aucun. Si vous sélectionnez Aucun, le nombre de lectures incluses dans le chutier sera la valeur de ce chutier. --EffectiveGenomeSize est le génome. Entrez la longueur (bp) de la partie cartographiable de. Pour mm10 (également appelé GRCm38), il s'agit de «2652783500». (Référence https://deeptools.readthedocs.io/en/latest/content/feature/ effectiveGenomeSize.html) Entrez la longueur du bac (bp) utilisée pour le calcul dans --binSize. Écrivez le nom du fichier de sortie dans -o.
Le calcul prend du temps, alors installez Genome Browser en attendant.
Installation d'IGV (Integrative Genomics Viewer)
Le navigateur de génomes est un outil qui visualise les résultats des séquences. Vous voyez souvent le signal XX-seq à une position spécifique du génome visualisé, n'est-ce pas? Voilà. Installons-le maintenant!
Téléchargez le programme d'installation de votre système d'exploitation sur la page de téléchargement IGV (https://software.broadinstitute.org/software/igv/download). IGV est une interface graphique (le graphique sort et fonctionne avec la souris et le clavier Si vous êtes un utilisateur Windows utilisant WSL, veuillez sélectionner la version de Windows ici. Lancez le programme d'installation téléchargé et installez selon les instructions. Ensuite, le raccourci IGV suivant sera créé sur le bureau. Masu.
Double-cliquez dessus pour le démarrer (il faut environ 30 secondes pour démarrer). Après le démarrage, la fenêtre suivante apparaîtra.
Maintenant que hg19 est chargé, téléchargeons et chargeons mm10. Cliquez sur la flèche vers le bas dans la case rouge sur l'écran ci-dessus et vous verrez "Plus ...". Cliquez dessus. Cliquez ensuite sur "Souris mm10", cochez "Séquence de téléchargement" en bas à gauche et cliquez sur "OK". Le téléchargement de mm10 commencera.
Il est temps que la «couverture bam» soit terminée ...?
Lorsque vous avez terminé, faites glisser et déposez Med1.bigwig dans la fenêtre IGV.
Avez-vous vu le profil Med1 ChIP-seq comme indiqué dans l'image ci-dessus? Dans cet état, vous regardez tous les chromosomes à vol d'oiseau, et comme vous ne connaissez pas les détails, entrez différents noms de gènes dans la fenêtre de recherche entourée de rouge. Volons vers cet emplacement du corps du gène. Voici juste un exemple.
- Klf4 * est un gène important pour maintenir la cellité souche, mais de cette façon, on peut voir que Med1, une protéine qui active la transcription, est liée aux cellules ES de souris. Comme vous pouvez le voir sur cette image, les gènes Refseq Je suis en train de remonter la piste. De plus, si l'échelle sur l'axe vertical de l'écran est étrange, cliquez avec le bouton droit de la souris à l'endroit où le nom "Med1.bigwig" est écrit et cliquez sur "Autoscale".
Créez un bigwig pour les deux données restantes de la même manière et vérifiez-le avec IGV.
Appel de pointe
Faisons maintenant un appel de pic pour détecter le pic du signal basé sur des critères statistiques. Cette fois, nous utiliserons findPeaks de HOMER. Un autre appelant de pic couramment utilisé est [MACS2](https :: //github.com/taoliu/MACS). Si vous êtes intéressé, veuillez comparer les résultats.
Maintenant, avant de faire findPeaks, nous devons convertir le fichier bam en une forme de TagDirectory que HOMER peut gérer. Pour ce faire, utilisez makeTagDirectory de HOMER.
makeTagDirectory Med1_data/Med1_TagDir -single Med1_data/Med1_rm_dups.bam
Écrivez le nom du TagDirectory à créer immédiatement après makeTagDirectory, puis écrivez les options et enfin écrivez le nom du fichier bam. Cette fois, l'option spécifiait uniquement l'option -single qui nettoie le TagDirectory. Voir http://homer.ucsd.edu/homer/ngs/tagDir.html pour les options de. Créez un répertoire de balises pour les deux données restantes de la même manière.
Maintenant, exécutons findPeaks.
findPeaks Med1_data/Med1_TagDir -style factor -o auto -i Input_data/Input_TagDir
Écrivez le nom du TagDirectory qui effectue l'appel de pointe immédiatement après findPeaks. Pour l'option -style, cette fois, nous avons considéré Med1 comme un facteur de transcription et entré factor. -o option Écrivez l'emplacement dans lequel écrire le résultat, mais si vous le définissez comme ʻauto, il sera enregistré dans le TagDirectory qui effectue l'appel de pointe. De plus, dans l'option -i, écrivez le nom du TagDirectory des données d'entrée. Données d'entrée S'il n'y a pas, si vous n'entrez pas l'option -i` elle-même, le calcul sera effectué sans entrée.
Ensuite, effectuez un appel de pointe pour les données de H3K27Ac ChIP-seq.
findPeaks H3K27Ac_data/H3K27Ac_TagDir -style histone -o auto -i Input_data/Input_TagDir
Puisque H3K27Ac est une donnée modifiée par histone, définissez l'option -style sur histone.
Jetons un coup d'œil au fichier de résultat de l'appel de crête. Je pense qu'il y a un fichier nommé «peaks.txt» dans TagDirectory (H3K27Ac est nommé «regions.txt»). Puisqu'il s'agit de texte, il est facile de voir si vous l'ouvrez dans Excel.
Une fois que diverses informations sont écrites en haut, les informations sur les pics comme la moitié inférieure de l'image sont écrites. Si vous regardez les colonnes «chr», «start», «end», le génome de chaque pic Vous pouvez voir la position ci-dessus. Normalized Tag Count est la force du signal de chaque pic. Divisez par 10 pour obtenir le rpm (Reads per Million). Pour plus d'informations sur la fonction de findPeaks, voir http: // homer. Voir ucsd.edu/homer/ngs/peaks.html.
Ensuite, regardons ce résultat dans IGV avec le fichier bigwig mentionné précédemment. Pour faciliter la lecture des informations de crête par IGV, il est préférable d'utiliser le format comme un fichier de lit. Le fichier de lit est très simple et un C'est un fichier texte comme l'image ci-dessous qui enregistre les informations de position du pic composé de 3 lignes à partir de l'extrême gauche, chr, start, ʻend`. (Des informations supplémentaires sont incluses après la 4ème ligne. Il peut y avoir.)
Vous pouvez le créer à partir de peaks.txt dans la sortie de findPeaks avec la commande suivante.
sed '/^#/d' Med1_data/Med1_TagDir/peaks.txt | awk -v 'OFS=\t' '{print $2, $3, $4}' > Med1_data/Med1_TagDir/peaks.bed
Essayez de faire glisser et déposer le peaks.bed que vous venez de créer sur l'IGV.
Vous devriez voir le pic appelé position comme ceci. Facile à comprendre ...
Analyse de motif combinée
Étant donné que les facteurs de transcription se lient à une certaine séquence (séquence de motifs de liaison), lorsque des régions de pic sont données, quel type de facteurs de transcription s'y lie en examinant le type de séquence qui y est enrichie. Par exemple, en effectuant une analyse de motif de liaison dans la région de pic de ChIP-seq de Med1, il est possible de prédire quel type de facteur de transcription se liera à l'endroit où Med1 est lié. Je vais essayer cela en utilisant findMotifsGenome.pl de HOMER.
Cela nécessite également une préparation et l'installation du génome dans HOMER.
perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10
Installez mm10 avec configureHomer.pl situé au-dessus du répertoire miniconda3. Ce qui précède est le cas lorsque miniconda3 est dans le répertoire personnel. De plus, pour la partie homomer-4.10-0, vous pouvezconda install Veuillez noter que cela peut changer en fonction de l'heure et de l'environnement où HOMER a été installé avec homer (car la version de HOMER peut être différente.) Pour les personnes WSL qui n'ont pas travaillé ici, l'aide ci-dessous Voir 2.
Exécutez maintenant findMotifsGenome.pl.
findMotifsGenome.pl Med1_data/Med1_TagDir/peaks.txt mm10 Med1_data/Med1_motif -size given -p 4
Écrivez le fichier de résultat de findPeaks immédiatement après findMotifsGenome.pl. (Si vous voulez utiliser le résultat d'autres appelants de pointe, reportez-vous à l'Aide 3 ci-dessous.) Ensuite, écrivez la version du génome. De plus, écrivez le nom du répertoire pour enregistrer le résultat. Le calcul est effectué pour la zone centrée sur le centre de chaque pic de la longueur entière saisie dans l'option -size. Et "taille 100", le calcul est effectué dans la plage de +/- 50 pb à partir du centre de chaque pic. Cette fois, le calcul est effectué dans la plage réelle de chaque zone de pic transmise au programme en tant que "taille donnée". Entrez le nombre de threads à utiliser dans -p.
Lorsque le calcul est terminé, il y aura un répertoire appelé «Med1_motif» dans «Med1_data», alors jetons un œil à l'intérieur. À noter: «knownResults.html». Ouvrez-le dans votre navigateur. Le calcul se termine ici. Si vous obtenez une erreur, reportez-vous à l'aide 2 ci-dessous.
L'écran ressemble à ceci: Dans le cas de ce calcul, la valeur p montre quel type de tableau de motifs de la base de données HOMER a été enrichi dans l'entrée de la zone de pic pour la zone aléatoire préparée automatiquement. Il est affiché dans l'ordre croissant. Comme le calcul est effectué pour une région aléatoire, le résultat change légèrement à chaque fois. En regardant ce résultat, la région de pic de Med1 sert à conserver les propriétés des cellules souches telles que KLF, OCT, SOX. Il existe de nombreux motifs de facteurs de transcription importants, suggérant que Med1 et ces facteurs de transcription sont co-localisés.
En dehors de cela, homerResults.html résume les séquences qui sont abondantes dans la zone de pic d'entrée et sont similaires au tableau de motifs de la base de données HOMER. En gros, vérifiez knownResults.html. Je devrais le faire.
Pour plus d'informations sur findMotifsGenome.pl, veuillez visiter http://homer.ucsd.edu/homer/ngs/peakMotifs.html.
Examiner le chevauchement entre les régions de pic
Enfin, je vais vous montrer comment vérifier le chevauchement entre les pics. Pour cela, utilisez le mergePeaks de HOMER. Cette fois, nous étudierons le chevauchement de l'aire des pics de Med1 ChIP-seq et de l'aire des pics de H3K27Ac ChIP-seq.
mergePeaks -d given Med1_data/Med1_TagDir/peaks.txt H3K27Ac_data/H3K27Ac_TagDir/regions.txt -prefix mergePeaks -venn venn.txt
Si «-d» est défini sur «donné», le chevauchement entre les zones de pic d'entrée sera calculé tel quel. Les zones de pic d'entrée seront écrites et arrangées. Il peut y en avoir 3 ou plus. «-Prefix XXX» Et la zone où les zones de chevauchement de chaque pic commençant par «XXX» sont combinées et la zone qui n'existe que dans une seule zone de pic est sortie séparément. Si vous la définissez comme «-venn YYY.txt» Il crée un tableau qui résume le nombre d'aires de pics qui se chevauchent appelé YYY.txt et le nombre d'aires de pics qui ne sont que l'une d'elles pour dessiner un diagramme de Ben. Pour plus de détails sur les options, etc. http: //homer.ucsd Voir .edu / homor / ngs / mergePeaks.html.
Lorsque cette commande est exécutée, mergePeaks_H3K27Ac_data_H3K27Ac_TagDir_regions.txt, mergePeaks_Med1_data_Med1_TagDir_peaks.txt, `mergePeaks_Med1_data_Med1_TagDir_peaks. , La région existant uniquement dans le pic de Med1 ChIP-seq, la région où les régions de pic se chevauchant de Med1 ChIP-seq et H3K27Ac ChIP-seq sont combinées, et le tableau pour dessiner le diagramme de Ben.
Dessinons un diagramme Ben avec matplotlib en Python. Puisque nous utilisons un package appelé matplotlib_venn,
conda install matplotlib-venn
Ouvrez ensuite l'éditeur et écrivez le code suivant, comme vous le souhaitez.
from matplotlib import pyplot as plt
from matplotlib_venn import venn2
#venn.Ouvrez txt et Med1 ChIP-Nombre de pics présents uniquement dans la séquence,
#H3K27Ac ChIP-Nombre de pics présents uniquement dans la séquence,Vérifiez le nombre de pics qui se chevauchent.
venn2(subsets=(770, 25254, 2738), set_labels = ("Med1", "H3K27Ac"))
#subsets=(Med1 uniquement,H3K27Ac uniquement,Chevauchement de pics)
plt.savefig("./venn.png ")
Si vous pouvez l'écrire, enregistrez-le sous un nom tel que venn_plot.py et exécutez la commande suivante à l'emplacement enregistré.
python venn_plot.py
Ensuite, un fichier appelé venn.png sera créé dans ce répertoire, alors ouvrez-le.
Environ 80% des pics de ChIP-seq de Med1, qui est une protéine qui active la transcription, se chevauchent avec les pics de ChIP-seq de H3K27Ac, qui est également un marqueur d'activité transcriptionnelle. Même ainsi, le nombre de pics de H3K27Ac ChIP-seq est important. Veuillez comparer les deux données avec IGV.
À la fin
Merci d'avoir lu jusqu'ici. L'analyse après l'appel de pointe introduite cette fois n'est qu'une partie de l'analyse ChIP-seq. Vous pouvez maintenant cartographier et atteindre un pic d'appel. Effectuez des analyses spécifiques en utilisant HOMER, R, Python, etc. Les informations sur la bioinformatique, y compris cet article, sont désormais abondantes sur le net. Suivez également la même procédure que cette fois pour l'état d'ouverture et de fermeture de la chromatine dans tout le génome. Il est également possible d'analyser ATAC-seq (Assay for Transposase-Accessible Chromatin Sequencing), qui est étudié de manière exhaustive. Nous espérons que cet article vous aidera dans vos recherches. Si vous avez des questions, veuillez nous aider autant que possible. Je suis désolé, alors j'aimerais avoir de vos nouvelles!
Aidez-moi
1. La commande conda ne fonctionne pas!
La cause est probablement que le chemin ne passe pas après l'installation de miniconda3. Veuillez procéder comme suit.
cd #Déplacer vers le répertoire de base
vim .bash_profile #Décrivez le chemin.bash_Profil ouvert avec vim
Lorsqu'il s'ouvre, appuyez d'abord sur la touche Échap. Ensuite, appuyez sur la touche I. Vous pouvez maintenant modifier en mode Insertion. Assurez-vous de saisir correctement les informations suivantes.
PATH=$PATH:~/miniconda3/bin
Écrivez le chemin vers miniconda3 / bin après PATH = $ PATH:. Ce qui précède est lorsque miniconda est dans votre répertoire personnel. Assurez-vous que vous n'avez pas fait d'erreur et appuyez à nouveau sur la touche Échap. Et: Tapez wq et appuyez sur Entrée.
Puis redémarrez le terminal
source .bash_profile
Cela devrait passer le chemin et exécuter conda.
2. configureHomer.pl -install mm10 ne fonctionne pas!
Probablement la seule erreur qui peut se produire pour les utilisateurs WSL, mais cela peut être dû au fait que les commandes nécessaires pour exécuter configureHomer.pl ne sont pas installées. Effectuez toutes les opérations suivantes:
which gcc
which g++
which make
which perl
which zip
which gzip
which wget
cette maison
/usr/bin/make
Si le chemin d'accès à la commande ne s'affiche pas comme
sudo apt install zip #Lorsque le zip n'était pas inclus
Veuillez lancer l'installation. Après avoir installé tout ce qui n'était pas là, réessayez
perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10
Essayez de courir.
3. Utilisez un fichier de lit créé par un programme autre que HOMER dans HOMER
Si vous souhaitez utiliser un fichier de lit créé par un programme autre que HOMER dans HOMER, vous pouvez le convertir à l'avance en fichier de lit HOMER. Pour ce faire, utilisez bed2pos.pl de HOMER.
bed2pos.pl (Fichier de lit que vous souhaitez convertir) > Converted_file.hb
L'extension du dossier de lit HOMER est "hb".
References
-
http://rnakato.hatenablog.jp/entry/2017/07/06/110926
-
https://bi.biopapyrus.jp/rnaseq/qc/trimmomatic.html
-
http://www.usadellab.org/cms/?page=trimmomatic
-
https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html
-
https://bi.biopapyrus.jp/rnaseq/mapping/bowtie2/
Recommended Posts