Dies ist das R-Paket für GRAPPLE (Genome-wide mR Analysis under Pervasive PLEiotropy), ein umfassendes Framework für die Mendelsche Randomisierung.
Einzelheiten zum GRAPPLE-Framework finden Sie in unserem Manuskript.
Das GRAPPLE-Paket kann in R aus diesem Github-Repository installiert werden:
library(devtools)
install_github("jingshuw/grapple")
Möglicherweise muss man auch PLINK für LD-Clumping installieren und einen geeigneten LD-Clumping-Referenzdatensatz für PLINK vorbereiten, wenn er/sie GRAPPLE für die Datenvorverarbeitung benötigt. Eine Ressource des Referenzdatensatzes der europäischen Bevölkerung zur Verklumpung, der auch im GRAPPLE-Papier verwendet wird, ist das hier heruntergeladene europäische Referenzpanel mit 1000 Genomen: http://fileserve.mrcieu.ac.uk/ld/data_maf0.01_rs_ref.tgz , aus dem MRCIEU-Zitat.
Die Analyse mit GRAPPLE beginnt mit der Vorverarbeitung der rohen GWAS-Zusammenfassungsstatistiken, die die Auswahl genetischer Instrumente, die Harmonisierung der Datensätze, um denselben Effekt und dasselbe Referenzallel zu teilen, die Berechnung einer Korrelationsmatrix des Rauschens in geschätzten Randassoziationen, die von allen SNPs gemeinsam genutzt wird, umfasst. und Auswählen einer Liste von Kandidaten-SNPs für die Modusmarkierungserkennung.
Dazu benötigt GRAPPLE eine Liste der rohen GWAS-Zusammenfassungsstatistikdateien. Jede GWAS-Datei ist eine „.csv“- oder „.txt“-Datei, die einen Datenrahmen mit mindestens 6 Spalten mit den folgenden Spaltennamen enthält: SNP
, effect_allele
, other_allele
, beta
, se
, pval
. Die SNP
Spalte enthält rsID für jeden SNP. Sowohl die Spalten effect_allele
als auch other_allele
müssen Großbuchstaben enthalten. Die beta
-Spalte enthält die geschätzte Effektgröße für kontinuierliche Merkmale und das Log-Odds-Verhältnis für binäre Merkmale, und die se
-Spalte ist die Standardabweichung des entsprechenden beta
.
Alle GWAS-Zusammenfassungsstatistikdateien im oben genannten Format, die im GRAPPLE-Papier verwendet werden, können aus Datensätzen heruntergeladen werden. Die öffentlich zugängliche Originalquelle jeder Datei finden Sie in Abschnitt 3 der Zusatzinformationen zum Manuskript.
Die Deta-Vorverarbeitung von GRAPPLE benötigt eine Datei für die Krankheit (Ergebnisdatei) und zwei Dateien für jeden Risikofaktor, wobei eine für die SNP-Auswahl (Auswahldatei) und die andere für die Schätzung der marginalen Auswirkungen der SNPs auf diesen Risikofaktor dient ( Belichtungsdatei). Wir erlauben überlappende Personen für die GWAS-Daten in den Ergebnis- und Expositionsdateien. Es sollten möglichst wenige überlappende Personen zwischen den Auswahldateien und anderen Dateien vorhanden sein, um SNP-Auswahlverzerrungen zu vermeiden. Manchmal kann es schwierig sein, solche GWAS-Datensätze für die SNP-Auswahl zu finden. Eine Option besteht darin, GWAS-Daten anderer Vorfahren zu verwenden.
Beispielsweise kann man die zusammenfassenden Statistikdateien für BMI und T2D aus Datensätzen herunterladen, außerdem das Referenzpanel herunterladen und die „.tgz“-Datei im aktuellen Ordner entpacken. Dann können wir damit beginnen, SNPs aus diesen Dateien als genetische Instrumente auszuwählen.
library(GRAPPLE)
sel.file <- "BMI-ukb.csv"
exp.file <- "BMI-giant17eu.csv"
out.file <- "T2D-diagram12-M.csv"
plink_refdat <- "./data_maf0.01_rs_ref/data_maf0.01_rs_ref"
Das Extrahieren der unabhängigen SNPs kann je nach Größe der GWAS-Zusammenfassungsstatistikdateien einige Minuten dauern.
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01, plink_exe = 'plink')
Diese Funktion extrahiert alle unabhängigen SNPs, deren Auswahl-p-Werte 0,01 nicht überschreiten. Es gibt drei Elemente zurück. Das eine ist „data“, ein Datenrahmen der zusammenfassenden Statistiken der ausgewählten SNPs, das andere ist „marker.data“, ein Datenrahmen für alle Kandidaten-Marker-SNPs. Der Unterschied zwischen „data“ und „marker.data“ besteht darin, dass wir bei der LD-Verklumpung für die Auswahl genetischer Instrumente einen strengeren r2 (0,001) verwenden, um die Unabhängigkeit zu gewährleisten, als bei der Auswahl von Kandidaten-Marker-SNPs (r2 = 0,05). Das dritte Element ist die geschätzte Korrelationsmatrix für die GWAS-Kohorten.
Wir verwenden PLINK, um LD-Clumping durchzuführen. Der Standardname/Pfad des plink-Befehls ist „plink“ (wird über das Argument plink_exe
übergeben). Für Benutzer mit Linux-/Windows-Systemen würde dieser Befehl nicht funktionieren und man muss den Pfad der exe-Datei angeben, z. B. „./plink“, je nachdem, wo sie plink installieren. Wenn R den plink-Befehl nicht ausführen kann, wird eine Fehlermeldung angezeigt, die besagt, dass die verklumpte Datei nicht gefunden wurde.
Wenn die Anzahl der Risikofaktoren 1 beträgt, kann GRAPPLE die Anzahl der pleiotropen Pfade ermitteln, indem es die Anzahl der Modi in der robustisierten Profilwahrscheinlichkeit ermittelt, wobei ein p-Wert-Auswahlschwellenwert gegeben ist.
Wir können das obige Beispiel ausprobieren.
## Here we take the p-value threshold be 1e-4, but one can try a series of p-values and see if the result is consistent
diagnosis <- findModes(data.list$data, p.thres = 1e-4)
diagnosis$modes
diagnosis$p
Es gibt nur einen Modus in der Wahrscheinlichkeit, daher finden wir keine Hinweise auf mehrere pleiotrope Wege. Es gibt uns eine gewisse Garantie, MR-RAPs zur Abschätzung der kausalen Wirkung des Risikofaktors zu verwenden.
result <- grappleRobustEst(data.list$data, p.thres = 1e-4)
Diese Funktion gibt unsere Schätzungsergebnisse zurück.
Hier ist ein weiteres Beispiel, bei dem wir mehrere Modi erkennen können
library(GRAPPLE)
sel.file <- "CRP-Prins17.csv"
exp.file <- "CRP-Dehghan11.csv"
out.file <- "CAD-Nelson17.csv"
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01)
diagnosis <- findModes(data.list$data, p.thres = 1e-4, marker.data = data.list$marker.data)
diagnosis$p
Wir finden drei Modi in der Profilwahrscheinlichkeit. Wenn wir die Marker-SNPs, Markergene und das kartierte GWAS-Merkmal überprüfen, können wir feststellen, dass LDL-C ein verwirrendes Merkmal sein kann. Anschließend können wir GRAPPLE ausführen, um die LDL-C-Effekte anzupassen
library(GRAPPLE)
sel.file <-c(sel.file, "LDL-gera18.csv")
exp.file <- c(exp.file, "LDL-glgc13.csv")
out.file <- "CAD-Nelson17.csv"
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01)
result <- grappleRobustEst(data.list$data, p.thres = 1e-4)
Nach Bereinigung um LDL-C können wir keine Hinweise darauf finden, dass es einen kausalen Effekt von CRP auf CAD gibt.