Einführung

ChIP-seq (Chromatin-Immunpräzipitation gefolgt von Sequenzierung) ist ein umfassendes Maß dafür, wo und wie oft spezifische Transkriptionsfaktorbindungen und Histonmodifikationen im Genom auftreten. Sie können mit der Umgebungskonstruktion beginnen, anhand der Daten des Papiers analysieren und schließlich die Ergebnisse im Genombrowser anzeigen und Spitzenaufrufe tätigen ~. Beginnen wir sofort mit der Umgebungskonstruktion!

Umgebung

Installation von miniconda3

Installieren Sie zunächst den Paketmanager miniconda. Der Paketmanager erleichtert die Installation und Verwaltung von Tools beim Erstellen der Umgebung.

Für Linux

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

bash Miniconda3-latest-Linux-x86_64.sh

Für Mac

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-MacOSX-x86_64.sh

bash Miniconda3-latest-MacOSX-x86_64.sh

Befolgen Sie anschließend die Anweisungen und drücken Sie die EINGABETASTE oder geben Sie "Ja" ein. Sie können "Ja" für alle Ihre Fragen eingeben. Wenn Sie fertig sind, schließen Sie das Terminal einmal. Nach dem Start fügen wir Miniconda einen Kanal hinzu. Gehen Sie daher wie folgt vor: __ Stellen Sie sicher, dass Sie dies in dieser Reihenfolge tun __.

conda config --add channels defaults
conda config --add channels conda-forge
conda config --add channels bioconda

Wenn Ihnen mitgeteilt wird, dass Sie hier keine Eigentumswohnung haben, lesen Sie Hilfe 1 unten.

Installation der notwendigen Werkzeuge

Jetzt werden wir Miniconda verwenden, um die erforderlichen Tools zu installieren. Dieses Mal werden wir die folgenden Tools installieren.

SRA Toolkit Tools zum Herunterladen von Sequenzdaten von SRA
Trimmomatic Adapter-Trimmwerkzeug
fastQC Ein Tool zur Qualitätskontrolle von Fastq-Dateien
Bowtie2 Zuordnungstool
Picard Tool zum Entfernen von PCR-Duplikaten
samtools An verschiedenen Stellen verwendete Tools
deepTools Tools zum Korrigieren von Daten
HOMER Verschiedene Tools

Installieren Sie mit dem Befehl conda von miniconda wie folgt.

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

Antworte "y" auf alle "Weiter ([y] / n)?". Die Umgebung ist jetzt bereit! Ich werde von nun an auf die Daten eingehen, aber um Verwirrung zu vermeiden, wird Folgendes geschrieben ** unter der Annahme, dass alle Befehle im selben Verzeichnis ausgeführt werden **. Lassen Sie uns die Daten abrufen ~.

Sequenzdaten herunterladen

Auf Daten prüfen

Laden Sie die Sequenzdaten von SRA (Sequence Read Archive) herunter. Diesmal [Kagey MH * et al., *](Https: // ChIP-seq-Daten für Med1 von Maus-ES-Zellen von www.nature.com/articles/nature09380) (GSM560348) und Creyghton * et al., * Wir verwenden auch die ChIP-seq-Daten für H3K27Ac aus Maus-ES-Zellen von /50/21931.long (GSM594579). Wir erhalten auch die Eingabesequenzdaten (GSM560357) aus demselben Papier wie Med1. () Die Zahlen in sind die GEO-Zugangsnummern für jede Daten. Der Gene Expression Omnibus (GEO) von NCBI ist die Nummer, die zum Durchblättern dieser Daten benötigt wird.

Bevor wir die Daten herunterladen, überprüfen wir zunächst die Existenz dieser Daten im Browser. Nehmen wir als Beispiel die Daten von Med1. Zuerst [GEO-Site](https: //www.ncbi) Öffnen Sie .nlm.nih.gov / geo /) in Ihrem Browser und geben Sie "GSM560348" in das Suchfenster ein, das im Bild unten rot eingeschlossen ist.

Bitte schauen Sie sich die verschiedenen Informationen zu diesen Daten genau an, wie im Bild unten gezeigt.

Klicken Sie nun auf die Nummer rechts neben der SRA unten auf dieser Seite.

Anschließend werden Sie zur folgenden Seite weitergeleitet. Die im Bild rot eingeschlossene Nummer wird als SRR-Nummer bezeichnet, die beim Herunterladen mit der Zugriffsnummer dieser Daten erforderlich ist.

Diese Zahlen sind irgendwo in das Papier geschrieben, damit Sie sehen können, wo sich die Daten befinden.

Laden Sie Daten von SRA herunter

Verwenden Sie den sratoolkit-Befehl fastq-dump. Es ist sehr einfach zu verwenden und wenn die Daten Single-Ended sind

fastq-Dump SRR-Nummer der Daten, die Sie herunterladen möchten

Für gepaarte Enden

fastq-Dump SRR-Nummer der Daten, die Sie herunterladen möchten--split-files

Dadurch sollte die FastQ-Datei, bei der es sich um die vom Sequenzer ausgegebenen Rohdaten handelt, in das Verzeichnis heruntergeladen werden, in dem Sie diesen Befehl ausgeführt haben.

Laden wir diese Daten herunter.

fastq-dump SRR058988 #Med1
fastq-dump SRR066767 #H3K27Ac
fastq-dump SRR058997 #Eingang

Sie können sie auch wie folgt schreiben und alle gleichzeitig herunterladen.

fastq-dump SRR058988 SRR066767 SRR058997

Dieser Prozess lädt die SRA-spezifische .sra-Datei im komprimierten Dateiformat herunter und konvertiert sie in eine Fastq-Datei. Es wird einige Zeit dauern. Seien Sie also geduldig. Verwenden Sie diese Zeit für meine Bioinformatik Wie wäre es mit Kommentaren für die breite Öffentlichkeit? ??

https://laborify.net/2019/11/30/michida-bioinformatics/

Wenn der Download abgeschlossen ist, wird der Dateiname durch die SRR-Nummer angegeben. Benennen Sie ihn daher um, um das Verständnis zu erleichtern.

mv SRR058988.fastq Med1.fastq
mv SRR066767.fastq H3K27Ac.fastq
mv SRR058997.fastq Input.fastq

Adapter trimmen

Lassen Sie uns nun die Sequenzergebnisse mit Trimmomatic bereinigen.

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

Obwohl es lang ist, wird es nicht unterbrochen, um Unfälle aufgrund von Unterbrechungen zu vermeiden. Geben Sie unmittelbar nach "trimmomatisch" an, ob die zu analysierenden Daten einseitig (SE) oder paarweise (PE) sind. In den nächsten "Gewinden" Geben Sie die Anzahl der zu verwendenden Threads an. -Phred33 ist ein Zauberspruch. Geben Sie ihn unbedingt ein. Geben Sie anschließend den Namen der zu schneidenden Datei und den Namen der Datei nach dem Zuschneiden ein.

Der Speicherort nach ILLUMINACLIP ist der Speicherort der Adaptersequenzinformationen, die sich im Verzeichnis miniconda3 befinden sollten. Schreiben Sie sie entsprechend Ihrem Miniconda3-Speicherort neu. (Bei der Installation müssen Sie nichts unternehmen. Es sollte sich wie folgt in Ihrem Home-Verzeichnis befinden.) Außerdem steht 2:30:10 für die zulässige Anzahl von Nichtübereinstimmungen, den Palindrom-Clip-Schwellenwert bzw. den einfachen Clip-Schwellenwert. Grundsätzlich denke ich, dass Sie sich damit nicht anlegen müssen. Außerdem bedeutet "LEADING: 3" und "TRAILING: 3", Basen mit einem Qualitätsfaktor von weniger als 3 vom Anfang bzw. Ende des Lesevorgangs zu entfernen. "SLIDING WINDOW: 4: 15" bedeutet alle 4 bp Sehen Sie sich den durchschnittlichen Qualitätsfaktor an und entfernen Sie die Teile, die kleiner als 15 sind. Das letzte "MINLEN: 36" bedeutet, diejenigen mit einer Leitungslänge von weniger als 36 aus der Analyse zu entfernen. Ich habe die "Schnellstart" -Einstellungen auf der Trimmomatic-Seite (http://www.usadellab.org/cms/?page=trimmomatic) verwendet. Wenn Sie fertig sind, wird eine Datei mit dem Namen Med1_trimmed.fastq generiert. Führen Sie die verbleibenden zwei Daten mit denselben Optionen aus.

QC der Fastq-Datei nach dem Zuschneiden

Verwenden Sie fastQC, um die Qualität der zugeschnittenen fastq-Datei zu steuern.

fastqc --threads 4 --nogroup -o . Med1_trimmed.fastq

Schreiben Sie die Anzahl der Threads mit "--threads" unmittelbar nach "fastqc". Wenn Sie die nächste "--nogroup" schreiben, wird auch der 3'end-Lesevorgang analysiert. Geben Sie das Ergebnis an "-o" aus Schreiben Sie das zu erledigende Verzeichnis. Schreiben Sie abschließend den Namen der zu prüfenden Datei.

Eine Datei mit dem Namen "Med1_trimmed_fastqc.html" wird in dem in "-o" angegebenen Verzeichnis erstellt. Öffnen Sie sie also mit einem Browser.

Wenn die Zusammenfassung links fast grün ist, ist die Qualität in Ordnung. Diese Daten sind zu sauber ... Erläutern Sie jeden Index unter https://bi.biopapyrus.jp/rnaseq/qc/fastqc.html Vielen Dank. Vielen Dank für Ihre fortgesetzte Unterstützung dieser Website, Bioinformatik.

Bitte machen Sie die gleiche Option für die verbleibenden zwei Daten.

Kartierung

Zum Schluss werden wir mit Bowtie2 kartieren! Und vorher müssen wir den Index des Genoms erstellen. Das heißt, wir müssen das Referenzgenom vorbereiten, das für die Kartierung benötigt wird.

Laden Sie das gesamte mm10-Array von der UCSC-Seite mit dem Befehl wget herunter. Ich habe beschlossen, einen Ordner mit dem Namen ref_genome zu erstellen und dort abzulegen. Masu.

mkdir ref_genome
cd ref_genome
wget http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz
tar -xzf chromFa.tar.gz

Auch dieses Mal werden wir zufällige und unbekannte Sequenzen und das mitochondriale Genom (chrM.fa) verwerfen.

rm *random.fa
rm chrUn*
rm chrM.fa

Verwenden Sie cat, um die verbleibende Datei in eine Datei mit dem Namen mm10.fa zu verwandeln.

cat *.fa > mm10.fa

Ich werde ein Verzeichnis namens mm10_index in der gleichen Hierarchie wie "ref_genome" erstellen und den Index dort speichern.

cd .. #Nun ref_Wenn Sie im Genom sind
mkdir mm10_index
bowtie2-build -f ./ref_genome/mm10.fa ./mm10_index/mm10_index --threads 4

bowtie2-build ist der Befehl zum Indizieren von Bowtie2. Schreiben Sie den Pfad zum Ort der Genomsequenz in die Option -f und dann den Pfad zum Index. Geben Sie die Anzahl der Threads mit --threads an In diesem Fall sollten Sie 6 Dateien mit dem Namen "mm10_index * .bt2" in dem zuvor erstellten Verzeichnis "mm10_index" haben. Sie müssen dies nur einmal tun.

Es braucht viel Zeit. Lassen Sie uns diese Zeit nutzen, um die R-Sprache zu lernen, die häufig für statistische Analysen in der Bioinformatik verwendet wird. (Danke!)

https://qiita.com/roadricefield/items/001c882f84dd093f4407

Ich werde R in diesem Artikel nicht verwenden, aber ...

Es braucht wirklich Zeit. Es ist eine gute Stufe, den Befehl kurz vor der Heimreise zu drücken und fortzufahren, wenn Sie morgen kommen. Ich werde heute schlafen gehen und mir die Animation ansehen.

.........................

Guten Morgen zusammen! Ich habe 7 Stunden gebraucht, weil ich vergessen habe, "--threads" anzugeben! Lassen Sie uns abbilden.

bowtie2 -p 4 -x ./mm10_index/mm10_index -U Med1_trimmed.fastq -S Med1.sam

-p ist die Anzahl der Threads, -x ist der Index, -U ist die zuzuordnende Fastq-Datei, -S ist der Name der Ausgabedatei. Sie wird einmal als Sam-Datei ausgegeben. Die Zuordnung nimmt übrigens auch viel Zeit in Anspruch. Es ist ungefähr 30 Minuten.

Wenn Sie fertig sind, konvertieren Sie die Sam-Datei mit samtools in eine BAM-Datei.

samtools view -b -o Med1.bam Med1.sam

Damit ist das Mapping abgeschlossen! Bitte machen Sie die gleiche Option für die verbleibenden zwei Daten.

Entfernung von PCR-Duplikaten (optional)

Verwenden Sie schließlich Picard, um die PCR-Duplikate zu entfernen, nicht unbedingt.

samtools sort Med1.bam > Med1_sorted.bam #Sie müssen die BAM-Datei sortieren, um Picard verwenden zu können.

picard MarkDuplicates I=Med1_sorted.bam O=Med1_rm_dups.bam M=Med1_report.txt REMOVE_DUPLICATES=true

Benennen Sie die BAM-Datei, für die Sie PCR-Duplikate entfernen möchten, in "I", den Namen der Datei nach dem Entfernen von Duplikaten in "O", und schreiben Sie den Namen in "M", da ein Bericht erstellt wird, der die Berechnungsergebnisse zusammenfasst. Masu.

Führen Sie die beiden anderen Daten mit denselben Optionen aus.

Das Verzeichnis, an dem ich arbeite, wird immer voller, daher organisiere ich die Daten hier. Med1 ChIP-seq-Daten befinden sich im Verzeichnis Med1_data, H3K27Ac ChIP-seq-Daten befinden sich im Verzeichnis H3K27Ac_data, Verschieben Sie die Daten in ein Verzeichnis namens "Input_data".

mkdir Med1_data
mv Med1* Med1_data

mkdir H3K27Ac_data
mv H3K27Ac* H3K27Ac_data

mkdir Input_data
mv Input* Input_data

Beobachten Sie die Ergebnisse mit einem Genombrowser

Erstellen Sie eine bigWig-Datei

Zuerst konvertieren wir die bam-Datei in ein Format namens bigWig, das mit einem Genombrowser leicht zu erkennen ist. Verwenden Sie dazu bamCoverage von deepTools. Das heißt, die Signalstärke von ChIP-seq für jeden Bin im Genom wird berechnet. Dazu müssen Sie zuerst eine bam.bai-Datei erstellen, bei der es sich um eine bam-Indexdatei handelt. Verwenden Sie daher samtools. Lass es uns machen.

samtools index Med1_data/Med1_rm_dups.bam

Dadurch wird eine Indexdatei mit dem Namen "Med1_rm_dups.bam.bai" erstellt.

Führen Sie nun bamCoverage aus. Stellen Sie sicher, dass Sie die bam-Datei und ihre bam.bai-Datei im selben Verzeichnis ablegen und ausführen.

bamCoverage -b Med1_data/Med1_rm_dups.bam -p 4 --normalizeUsing RPGC --effectiveGenomeSize 2652783500 --binSize 1 -o Med1_data/Med1.bigwig

Schreiben Sie den Namen der BAM-Datei, die in eine BigWig-Datei konvertiert werden soll, in -b. -p ist die Anzahl der Threads. --NormalizeUsing wählt den Typ des Korrekturwerts aus, der in jedem Bin berechnet werden soll. RPKM, Sie können CPM, BPM, RPGC, None auswählen. Wenn Sie None auswählen, entspricht die Anzahl der im Bin enthaltenen Lesevorgänge dem Wert in diesem Bin. --EffectiveGenomeSize ist das Genom. Geben Sie die Länge (bp) des abbildbaren Teils von ein. Für mm10 (auch als GRCm38 bekannt) ist es "2652783500". (Referenz https://deeptools.readthedocs.io/en/latest/content/feature/ effektivGenomeSize.html) Geben Sie die für die Berechnung verwendete Bin-Länge (bp) in --binSize ein. Schreiben Sie den Namen der Ausgabedatei in -o.

Die Berechnung braucht Zeit, installieren Sie also in der Zwischenzeit den Genombrowser.

Installation von IGV (Integrative Genomics Viewer)

Der Genombrowser ist ein Werkzeug, das Sequenzergebnisse visualisiert. Sie sehen häufig das XX-seq-Signal an einer bestimmten Position auf dem visualisierten Genom, oder? Das war's. Lass es uns jetzt installieren!

Laden Sie das Installationsprogramm für Ihr Betriebssystem von der IGV-Download-Seite (https://software.broadinstitute.org/software/igv/download) herunter. IGV ist eine grafische Benutzeroberfläche (Grafik wird ausgegeben und mit Maus und Tastatur bedient Wenn Sie ein Windows-Benutzer sind, der WSL verwendet, wählen Sie hier die Windows-Version aus. Starten Sie das heruntergeladene Installationsprogramm und befolgen Sie die Anweisungen zur Installation. Anschließend wird die folgende IGV-Verknüpfung auf dem Desktop erstellt. Masu.

Doppelklicken Sie darauf, um es zu starten (der Start dauert ca. 30 Sekunden). Nach dem Start wird das folgende Fenster angezeigt.

Nachdem hg19 geladen wurde, laden wir mm10 herunter und laden es. Klicken Sie auf den Abwärtspfeil im roten Feld auf dem Bildschirm oben und Sie sehen "Mehr ...". Klicken Sie darauf. Klicken Sie dann auf "Maus mm10", aktivieren Sie unten links "Download-Sequenz" und klicken Sie auf "OK". Dadurch wird der Download von mm10 gestartet.

Es ist an der Zeit, dass die "bam Coverage" vorbei ist ...?

Wenn Sie fertig sind, ziehen Sie "Med1.bigwig" per Drag & Drop in das IGV-Fenster.

Haben Sie das Med1 ChIP-seq-Profil wie im obigen Bild gezeigt gesehen? In diesem Zustand betrachten Sie alle Chromosomen aus der Vogelperspektive. Da Sie die Details nicht kennen, geben Sie verschiedene Gennamen in das von Rot umgebene Suchfenster ein. Fliegen wir zu diesem Ort des Genkörpers. Hier ist nur ein Beispiel.

Klf4 * ist ein wichtiges Gen für die Aufrechterhaltung der Stammzelligkeit, aber auf diese Weise kann man sehen, dass Med1, ein Protein, das die Transkription aktiviert, an Maus-ES-Zellen gebunden ist. Wie Sie in diesem Bild sehen können, sind Refseq-Gene Ich ziehe die Spur hoch. Wenn die Skalierung auf der vertikalen Achse des Displays seltsam ist, klicken Sie mit der rechten Maustaste auf die Stelle, an der der Name "Med1.bigwig" geschrieben ist, und klicken Sie auf "Autoscale".

Erstellen Sie auf die gleiche Weise eine Bigwig für die verbleibenden zwei Daten und überprüfen Sie sie mit IGV.

Spitzenanruf

Lassen Sie uns nun einen Peak-Aufruf durchführen, um den Peak des Signals anhand statistischer Kriterien zu ermitteln. Dieses Mal verwenden wir "findPeaks" von HOMER. Ein weiterer häufig verwendeter Peak-Aufrufer ist [MACS2](https :: //github.com/taoliu/MACS). Wenn Sie interessiert sind, vergleichen Sie bitte die Ergebnisse.

Bevor wir "findPeaks" ausführen, müssen wir die BAM-Datei in eine Form von TagDirectory konvertieren, die HOMER verarbeiten kann. Verwenden Sie dazu das "makeTagDirectory" von HOMER.

makeTagDirectory Med1_data/Med1_TagDir -single Med1_data/Med1_rm_dups.bam

Schreiben Sie den Namen des TagDirectory, das unmittelbar nach "makeTagDirectory" erstellt werden soll, schreiben Sie dann die Optionen und schließlich den Namen der BAM-Datei. Dieses Mal gab diese Option nur die Option "-single" an, mit der das TagDirectory bereinigt wird. Optionen für finden Sie unter http://homer.ucsd.edu/homer/ngs/tagDir.html. Erstellen Sie auf die gleiche Weise ein Tag-Verzeichnis für die verbleibenden zwei Daten.

Lassen Sie uns nun findPeaks ausführen.

findPeaks Med1_data/Med1_TagDir -style factor -o auto -i Input_data/Input_TagDir

Schreiben Sie den Namen des TagDirectory, das den Peak-Aufruf unmittelbar nach findPeaks ausführt. Für die Option -style haben wir diesmal Med1 als Transkriptionsfaktor betrachtet und die Option factor. -o eingegeben Schreiben Sie den Speicherort, in den das Ergebnis geschrieben werden soll. Wenn Sie es jedoch als "auto" festlegen, wird es im TagDirectory gespeichert, das den Spitzenaufruf ausführt. Schreiben Sie außerdem in der Option "-i" den Namen des TagDirectory der Eingabedaten. Geben Sie Daten ein Wenn dies nicht der Fall ist und Sie die Option "-i" nicht selbst eingeben, wird die Berechnung ohne Eingabe durchgeführt.

Führen Sie als Nächstes einen Spitzenaufruf für die Daten von H3K27Ac ChIP-seq durch.

findPeaks H3K27Ac_data/H3K27Ac_TagDir -style histone -o auto -i Input_data/Input_TagDir

Da es sich bei H3K27Ac um Histon-modifizierte Daten handelt, setzen Sie die Option "-style" auf "Histon".

Werfen wir einen Blick auf die Ergebnisdatei des Peak-Aufrufs. Ich glaube, es gibt eine Datei mit dem Namen "eaks.txt "in TagDirectory (H3K27Ac heißt" region.txt "). Da es sich um Text handelt, ist es leicht zu erkennen, ob Sie ihn in Excel öffnen.

Nachdem verschiedene Informationen oben geschrieben wurden, werden Peakinformationen wie die untere Bildhälfte geschrieben. Wenn Sie sich die Spalten "chr", "start", "end" ansehen, das Genom jedes Peaks Sie können die Position oben sehen. Normalized Tag Count ist die Signalstärke jedes Peaks. Teilen Sie durch 10, um die Drehzahl (Reads per Million) zu erhalten. Weitere Informationen zur Funktion von findPeaks finden Sie unter http: // homer. Siehe ucsd.edu/homer/ngs/peaks.html.

Als nächstes betrachten wir dieses Ergebnis in IGV zusammen mit der zuvor erwähnten Bigwig-Datei. Um IGV das Lesen von Spitzeninformationen zu erleichtern, ist es besser, das Format als Bettdatei zu verwenden. Die Bettdatei ist sehr einfach und eins. Es ist eine Textdatei wie das folgende Bild, die die Positionsinformationen des Peaks aufzeichnet, die aus drei Zeilen ganz links bestehen: "chr", "start", "end". (Zusätzliche Informationen sind nach der 4. Zeile enthalten. Es kann sein.)

Sie können es aus eaks.txt in der Ausgabe von findPeaks mit dem folgenden Befehl erstellen.

sed '/^#/d' Med1_data/Med1_TagDir/peaks.txt | awk -v 'OFS=\t' '{print $2, $3, $4}' > Med1_data/Med1_TagDir/peaks.bed

Versuchen Sie, das soeben erstellte eaks.bed auf das IGV zu ziehen und dort abzulegen.

Sie sollten den Peak mit der Bezeichnung Position so sehen. Leicht zu verstehen ...

Kombinierte Motivanalyse

Da Transkriptionsfaktoren an eine bestimmte Sequenz binden (Bindungsmotivsequenz), binden bei Angabe von Peakregionen welche Art von Transkriptionsfaktoren an sie, indem untersucht wird, welche Art von Sequenz in ihnen angereichert ist. Zum Beispiel ist es durch Durchführen einer Bindungsmotivanalyse in der Peakregion von ChIP-seq von Med1 möglich, vorherzusagen, welche Art von Transkriptionsfaktor an den Ort bindet, an dem Med1 gebunden ist. Ich werde dies mit HOMERs findMotifsGenome.pl versuchen.

Dies erfordert auch eine Vorbereitung und die Installation des Genoms in HOMER.

perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10

Installieren Sie mm10 mit configureHomer.pl über dem Verzeichnis miniconda3. Dies ist der Fall, wenn sich miniconda3 im Ausgangsverzeichnis befindet. Außerdem können Sie für den Teil homomer-4.10-0`conda installieren Bitte beachten Sie, dass sich dies je nach Zeit und Umgebung ändern kann, als HOMER mit homer installiert wurde. (Da die Version von HOMER möglicherweise unterschiedlich ist.) Für WSL-Benutzer, die hier nicht gearbeitet haben, die folgende Hilfe Siehe 2.

Führen Sie nun findMotifsGenome.pl aus.

findMotifsGenome.pl Med1_data/Med1_TagDir/peaks.txt mm10 Med1_data/Med1_motif -size given -p 4

Schreiben Sie die Ergebnisdatei von findPeaks unmittelbar nach findMotifsGenome.pl. (Wenn Sie das Ergebnis anderer Peak-Aufrufer verwenden möchten, lesen Sie Hilfe 3 unten.) Schreiben Sie dann die Version des Genoms. Schreiben Sie außerdem den Namen des Verzeichnisses, um das Ergebnis zu speichern. Die Berechnung wird für den Bereich durchgeführt, der in der Mitte jedes Peaks der in der Option "-size" eingegebenen ganzzahligen Länge zentriert ist. Es ist schwer zu verstehen. Geben Sie daher ein Beispiel an. Und "Größe 100", die Berechnung wird im Bereich von +/- 50 bp von der Mitte jedes Peaks durchgeführt. Dieses Mal wird die Berechnung im tatsächlichen Bereich jeder Peakfläche durchgeführt, die als "Größe angegeben" an das Programm übergeben wird. Geben Sie die Anzahl der Threads ein, die in -p verwendet werden sollen.

Wenn die Berechnung abgeschlossen ist, gibt es in "Med1_data" ein Verzeichnis mit dem Namen "Med1_motif". Schauen wir uns also das Verzeichnis an. Bemerkenswert ist "unknownResults.html". Öffnen Sie dieses Verzeichnis in Ihrem Browser. Die Berechnung endet hier. Wenn Sie eine Fehlermeldung erhalten, lesen Sie bitte Hilfe 2 unten.

Der Bildschirm sieht folgendermaßen aus. Bei dieser Berechnung zeigt der p-Wert, welche Art von Motivarray in der HOMER-Datenbank in der für den automatisch vorbereiteten Zufallsbereich eingegebenen Peakfläche angereichert wurde. Es wird in aufsteigender Reihenfolge angezeigt. Da die Berechnung für eine zufällige Region durchgeführt wird, ändert sich das Ergebnis jedes Mal geringfügig. Bei Betrachtung dieses Ergebnisses dient die Peakregion von Med1 zur Aufrechterhaltung von Stammzelleigenschaften wie KLF, OCT, SOX. Es gibt viele wichtige Transkriptionsfaktormotive, die darauf hindeuten, dass Med1 und diese Transkriptionsfaktoren gleichzeitig lokalisiert sind.

Abgesehen davon fasst homerResults.html die Sequenzen zusammen, die im Eingabepeakbereich häufig vorkommen und dem Motivarray in der HOMER-Datenbank ähnlich sind. Überprüfen Sie grundsätzlich unknownResults.html. Ich sollte es tun.

Weitere Informationen zu findMotifsGenome.pl finden Sie unter http://homer.ucsd.edu/homer/ngs/peakMotifs.html.

Untersuchen Sie die Überlappung zwischen Peakregionen

Abschließend werde ich Ihnen zeigen, wie Sie die Überlappung zwischen Peaks überprüfen. Verwenden Sie dazu die mergePeaks von HOMER. Dieses Mal untersuchen wir die Überlappung zwischen der Peakfläche von Med1 ChIP-seq und der Peakfläche von H3K27Ac ChIP-seq.

mergePeaks -d given Med1_data/Med1_TagDir/peaks.txt H3K27Ac_data/H3K27Ac_TagDir/regions.txt -prefix mergePeaks -venn venn.txt

Wenn "-d" auf "gegeben" gesetzt ist, wird die Überlappung zwischen den Eingangsspitzenbereichen so berechnet, wie sie ist. Die Eingangsspitzenbereiche werden geschrieben und angeordnet. Es können 3 oder mehr vorhanden sein. "-Prefix XXX" Und der Bereich, in dem die überlappenden Bereiche jedes Peaks, beginnend mit "XXX", kombiniert werden und der Bereich, der nur in einem Peakbereich vorhanden ist, wird separat ausgegeben. Wenn Sie ihn als "-venn YYY.txt" festlegen Es wird eine Tabelle erstellt, in der die Anzahl der überlappenden Peakflächen mit der Bezeichnung "YYY.txt" und die Anzahl der Peakflächen, die nur eine davon sind, zusammengefasst werden, um ein Ben-Diagramm zu zeichnen. Einzelheiten zu Optionen usw. finden Sie unter http: //homer.ucsd Siehe .edu / homor / ngs / mergePeaks.html.

Wenn dieser Befehl ausgeführt wird, wird "mergePeaks_H3K27Ac_data_H3K27Ac_TagDir_regions.txt", "mergePeaks_Med1_data_Med1_TagDir_peaks.txt", "mergePeaks_Med1_data_Med1_TagDir_peaks.txt_H3" angezeigt , Die Region, die nur im Peak von Med1 ChIP-seq existiert, die Region, in der die überlappenden Peakregionen von Med1 ChIP-seq und H3K27Ac ChIP-seq kombiniert werden, und die Tabelle zum Zeichnen des Ben-Diagramms.

Zeichnen wir ein Ben-Diagramm mit matplotlib in Python. Da wir ein Paket namens matplotlib_venn verwenden,

conda install matplotlib-venn

Öffnen Sie dann den Editor und schreiben Sie den folgenden Code, was auch immer Sie wollen.

from matplotlib import pyplot as plt
from matplotlib_venn import venn2

#venn.Öffnen Sie txt und Med1 ChIP-Anzahl der nur in seq, 
#H3K27Ac ChIP-Anzahl der nur in seq,Überprüfen Sie die Anzahl der überlappenden Peaks.

venn2(subsets=(770, 25254, 2738), set_labels = ("Med1", "H3K27Ac"))
#subsets=(Nur Med1,Nur H3K27Ac,Überlappende Spitzen)

plt.savefig("./venn.png ")

Wenn Sie es schreiben können, speichern Sie es unter einem Namen wie "venn_plot.py" und führen Sie den folgenden Befehl am gespeicherten Speicherort aus.

python venn_plot.py

Dann wird eine Datei mit dem Namen "venn.png " in diesem Verzeichnis erstellt. Öffnen Sie sie also.

Etwa 80% der Peaks von ChIP-seq von Med1, einem Protein, das die Transkription aktiviert, überlappen sich mit den Peaks von ChIP-seq von H3K27Ac, das auch ein Marker für die Transkriptionsaktivität ist. Trotzdem ist die Anzahl der Peaks von H3K27Ac ChIP-seq groß. Bitte vergleichen Sie die beiden Daten mit IGV.

Am Ende

Vielen Dank, dass Sie bisher gelesen haben. Die diesmal eingeführte Analyse nach dem Spitzenanruf ist nur ein Teil der ChIP-seq-Analyse. Sie können jetzt den Spitzenanruf abbilden und anzeigen. Führen Sie eine zweckgebundene Analyse mit HOMER, R, Python usw. durch. Informationen zur Bioinformatik, einschließlich dieses Artikels, sind jetzt im Internet reichlich vorhanden. Befolgen Sie außerdem das gleiche Verfahren wie diesmal für den Öffnungs- und Schließzustand von Chromatin im gesamten Genom. Es ist auch möglich, ATAC-seq (Assay for Transposase-Accessible Chromatin Sequencing) zu analysieren, das ausführlich untersucht wird. Wir hoffen, dass dieser Artikel Ihrer Forschung hilft. Wenn Sie Fragen haben, helfen Sie uns bitte so weit wie möglich. Es tut mir leid, also würde ich gerne von Ihnen hören!

Hilfe

1. Der Befehl `conda` funktioniert nicht!

Wahrscheinlich liegt die Ursache darin, dass der Pfad nach der Installation von miniconda3 nicht übergeben wird. Führen Sie die folgenden Schritte aus.

cd #In das Ausgangsverzeichnis wechseln

vim .bash_profile #Beschreibe den Pfad.bash_Profil mit vim öffnen

Drücken Sie beim Öffnen zuerst die Esc-Taste und dann die i-Taste. Jetzt können Sie im Einfügemodus bearbeiten. Stellen Sie sicher, dass Sie Folgendes korrekt eingeben.

PATH=$PATH:~/miniconda3/bin

Schreiben Sie den Pfad zu miniconda3 / bin nach PATH = $ PATH:. Das Obige ist, wenn sich miniconda in Ihrem Home-Verzeichnis befindet. Stellen Sie sicher, dass Sie keinen Fehler gemacht haben, und drücken Sie die Esc-Taste erneut. Und: Geben Sie wq ein und drücken Sie die Eingabetaste.

Starten Sie dann das Terminal neu

source .bash_profile

Dies sollte den Pfad passieren und "conda" ausführen.

2. `configureHomer.pl -install mm10` funktioniert nicht!

Wahrscheinlich der einzige Fehler, der bei WSL-Benutzern auftreten kann, aber möglicherweise, weil die zum Ausführen von "configureHomer.pl" erforderlichen Befehle nicht installiert sind. Führen Sie alle folgenden Schritte aus:

which gcc
which g++
which make
which perl
which zip
which gzip
which wget

dieses Haus

/usr/bin/make

Wenn der Pfad zum Befehl nicht wie angezeigt wird

sudo apt install zip #Wenn der Reißverschluss nicht enthalten war

Führen Sie die Installation aus. Versuchen Sie es erneut, nachdem Sie alles installiert haben, was nicht vorhanden war

perl $HOME/miniconda3/share/homer-4.10-0/configureHomer.pl -install mm10

Versuche zu rennen.

3. Verwenden Sie eine Bettdatei, die von einem anderen Programm als HOMER in HOMER erstellt wurde

Wenn Sie eine Bettdatei verwenden möchten, die von einem anderen Programm als HOMER in HOMER erstellt wurde, können Sie diese sicher im Voraus in eine HOMER-Bettdatei konvertieren. Verwenden Sie dazu die "bed2pos.pl" von HOMER.

bed2pos.pl (Bettdatei, die Sie konvertieren möchten) > Converted_file.hb

Die Erweiterung der HOMER-Bettdatei lautet "hb".

References

Kagey MH et al., Mediator and cohesin connect gene expression and chromatin architecture. Nature 467, 430–435 (2010).
Creyghton et al., Histone H3K27ac separates active from poised enhancers and predicts developmental state. PNAS 107, 21931–21936 (2010).
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/

[PYTHON] ChIP-seq-Analyse ab Null