Mengoreksi nilai p untuk beberapa tes di mana tes berkorelasi (genetika)

24

Saya memiliki nilai p dari banyak tes dan ingin tahu apakah sebenarnya ada sesuatu yang signifikan setelah mengoreksi beberapa pengujian. Komplikasinya: tes saya tidak independen. Metode yang saya pikirkan (varian dari Fisher's Product Method, Zaykin et al., Genet Epidemiol , 2002) membutuhkan korelasi antara nilai-nilai p.

Untuk memperkirakan korelasi ini, saya sedang berpikir tentang kasus bootstrap, menjalankan analisis dan menghubungkan vektor yang dihasilkan dari nilai p. Adakah yang punya ide yang lebih baik? Atau bahkan ide yang lebih baik untuk masalah asli saya (mengoreksi beberapa pengujian dalam tes berkorelasi)?

Latar belakang: Saya secara logistik mengalami kemunduran apakah subyek saya menderita penyakit tertentu atau tidak pada interaksi antara genotipe mereka (AA, Aa atau aa) dan kovariat. Namun, genotipe ini sebenarnya banyak (30-250) dari Single Nucleotide Polymorphisms (SNPs), yang tentunya tidak independen tetapi dalam Linkage Disequilibrium.

S. Kolassa - Reinstate Monica
sumber

Jawaban:

29

Ini sebenarnya topik hangat dalam studi analisis Genomewide (GWAS)! Saya tidak yakin metode yang Anda pikirkan adalah yang paling tepat dalam konteks ini. Kumpulan nilai-p dijelaskan oleh beberapa penulis, tetapi dalam konteks yang berbeda (studi replikasi atau meta-analisis, lihat misalnya (1) untuk ulasan terbaru). Menggabungkan nilai p SNP dengan metode Fisher umumnya digunakan ketika seseorang ingin mendapatkan nilai p unik untuk gen yang diberikan; ini memungkinkan untuk bekerja pada tingkat gen, dan mengurangi jumlah dimensi pengujian berikutnya, tetapi seperti yang Anda katakan ketidakbergantungan antara penanda (yang timbul dari colocation spasial atau disiquilibrium keterkaitan, LD) menimbulkan bias. Alternatif yang lebih kuat bergantung pada prosedur resampling,

Kekhawatiran utama saya dengan bootstraping (dengan penggantian) adalah bahwa Anda memperkenalkan bentuk keterkaitan buatan, atau dengan kata lain Anda membuat kembar virtual, karenanya mengubah keseimbangan Hardy-Weinberg (tetapi juga frekuensi alel minimum dan tingkat panggilan). Ini tidak akan menjadi kasus dengan pendekatan permutasi di mana Anda mengubah label masing-masing dan menjaga data genotyping seperti apa adanya. Biasanya, perangkat lunak plink dapat memberi Anda nilai-p yang baku dan diijinkan, meskipun menggunakan (secara default) strategi pengujian adaptif dengan jendela geser yang memungkinkan untuk berhenti menjalankan semua permutasi (katakanlah 1000 per SNP) jika tampak bahwa SNP di bawah pertimbangannya tidak "menarik"; itu juga memiliki opsi untuk menghitung maxT, lihat bantuan online .

Tetapi mengingat rendahnya jumlah SNP yang Anda pertimbangkan, saya sarankan mengandalkan tes berbasis FDR atau maxT seperti yang diterapkan dalam paket multtest R (lihat mt.maxT), tetapi panduan definitif untuk menyusun kembali strategi untuk aplikasi genom adalah Prosedur Pengujian Berganda dengan Aplikasi untuk Genomics , dari Dudoit & van der Laan (Springer, 2008). Lihat juga buku Andrea Foulkes tentang genetika dengan R , yang diulas dalam JSS. Dia memiliki materi yang bagus tentang berbagai prosedur pengujian.

Catatan selanjutnya

Banyak penulis telah menunjukkan fakta bahwa beberapa metode koreksi pengujian sederhana seperti Bonferroni atau Sidak terlalu ketat untuk menyesuaikan hasil untuk masing-masing SNP. Selain itu, tak satu pun dari metode ini memperhitungkan korelasi yang ada antara SNP karena LD yang menandai variasi genetik di seluruh wilayah gen. Alternatif lain telah diusulkan, seperti turunan dari metode Holm untuk perbandingan berganda (3), Hidden Markov Model (4), FDR bersyarat atau positif (5) atau turunannya (6), untuk menyebutkan beberapa. Apa yang disebut statistik gap atau jendela geser telah terbukti berhasil dalam beberapa kasus, tetapi Anda akan menemukan ulasan yang bagus dalam (7) dan (8).

Saya juga pernah mendengar tentang metode yang memanfaatkan struktur haplotype atau LD secara efektif, misalnya (9), tetapi saya tidak pernah menggunakannya. Namun, mereka tampaknya lebih terkait dengan memperkirakan korelasi antara penanda, bukan nilai p seperti yang Anda maksudkan. Tetapi pada kenyataannya, Anda mungkin berpikir lebih baik dalam hal struktur ketergantungan antara statistik tes berturut-turut, daripada antara p-nilai berkorelasi.

Referensi

  1. Cantor, RM, Lange, K dan Sinsheimer, JS. Memprioritaskan Hasil GWAS: Tinjauan Metode Statistik dan Rekomendasi untuk Aplikasi Mereka . Am J Hum Genet. 2010 86 (1): 6-22.
  2. Corley, RP, Zeiger, JS, Crowley, T et al. Asosiasi gen kandidat dengan ketergantungan obat antisosial pada remaja . Ketergantungan Obat dan Alkohol 2008 96: 90–98.
  3. Dalmasso, C, Génin, E dan Trégouet DA. Sebuah Prosedur Tertimbang-Holm Akuntansi untuk Frekuensi Allele dalam Studi Asosiasi Genomewide . Genetika 2008 180 (1): 697–702.
  4. Wei, Z, Sun, W, Wang, K, dan Hakonarson, H. Beberapa Pengujian dalam Studi Asosiasi Genome-Wide melalui Hidden Markov Models . Bioinformatika 2009 25 (21): 2802-2808.
  5. Broberg, P. Tinjauan perbandingan estimasi proporsi gen yang tidak berubah dan tingkat penemuan palsu . BMC Bioinformatics 2005 6: 199.
  6. Perlu, AC, Ge, D, Weale, ME, dan lain-lain. Investigasi Genom-Wide SNP dan CNV di Skizofrenia . Geno PLoS. 2009 5 (2): e1000373.
  7. Han, B, Kang, HM, dan Eskin, E. Koreksi Pengujian Berganda dan Estimasi Cepat dan Akurat untuk Jutaan Penanda Terkait . PLoS Genetics 2009
  8. Liang, Y dan Kelemen, A. Kemajuan statistik dan tantangan untuk menganalisis data snp dimensi tinggi berkorelasi dalam studi genom untuk penyakit kompleks . Survei Statistik 2008 2: 43–60. - ulasan terbaru terbaik yang pernah ada
  9. Nyholt, DR. Koreksi Sederhana untuk Pengujian Berganda untuk Polimorfisme Nukleotida-Tunggal dalam Ketidakseimbangan Tautan Satu Sama Lain . Am J Hum Genet. 2004 74 (4): 765-769.
  10. Nikodemus, KK, Liu, W, Chase, GA, Tsai, YY, dan Fallin, MD. Perbandingan kesalahan tipe I untuk koreksi beberapa pengujian dalam studi polimorfisme nukleotida tunggal besar menggunakan komponen utama versus algoritma pemblokiran haplotype . BMC Genetics 2005; 6 (Suppl 1): S78.
  11. Tes interval kepercayaan bootstrap bootstrap berbasis Peng, Q, Zhao, J, dan Xue, F. untuk asosiasi penyakit gen yang melibatkan banyak SNP . BMC Genetics 2010, 11: 6
  12. Li, M, Romero, R, Fu, WJ, dan Cui, Y (2010). Memetakan Interaksi Haplotype-haplotype dengan LASSO Adaptif . BMC Genetics 2010, 11:79 - meskipun tidak terkait langsung dengan pertanyaan, itu mencakup analisis berbasis haplotype / efek epistatik
chl
sumber
1
Wow, terima kasih sudah membantu semua masalah ini! Saya mengerti keraguan Anda tentang bootstrap, dan saya hampir yakin. Saya pikir komplikasi utama saya adalah kovariat numerik yang saya miliki yang tentunya akan diperlukan (baik dengan sendirinya atau dalam interaksi dengan genotipe), dan yang tampaknya mengesampingkan mt.maxT dan mem-plink, walaupun saya mungkin perlu memeriksanya lagi. Tapi saya pasti akan menggali referensi yang Anda berikan!
S. Kolassa - Reinstate Monica
Anda selalu dapat bekerja dengan residu GLM Anda untuk mendapatkan kovariat, meskipun Anda kehilangan beberapa Df yang akan sulit untuk diperhitungkan atau diperkenalkan kembali sesudahnya (misalnya untuk menghitung nilai-p).
chl
Hm, residu dari regresi logistik saya? Apakah itu sah?
S. Kolassa - Reinstate Monica
Ya kenapa tidak Tidak jarang menghapus varians yang diperhitungkan oleh kovariat lain dan kemudian melanjutkan analisis tingkat ke-2 dengan data residual Anda. Seringkali lebih cepat (misalnya, plink cukup lambat dengan kovariat kategoris, sementara itu ok dengan yang berkelanjutan; snpMatrixatau hanya glm()berkinerja lebih baik pada titik ini, tetapi Anda tidak dapat menanamkan banyak SNP dalam glm()...); masalahnya adalah mendapatkan nilai-p yang dikoreksi pada akhir analisis ke-2 Anda agak rumit (karena Anda harus memperhitungkan parameter yang sudah diperkirakan).
chl
Untuk ilustrasi tentang bagaimana orang bekerja dengan residu, lihat misalnya hal. 466 dari Heck et al. Investigasi 17 kandidat gen untuk sifat kepribadian mengkonfirmasi efek dari gen HTR2A pada pencarian baru. Gen, otak, dan perilaku (2009) vol. 8 (4) hlm. 464-72
chl
2

Menggunakan metode seperti bonferroni baik-baik saja, masalahnya adalah jika Anda memiliki banyak tes Anda tidak akan menemukan banyak "penemuan".

Anda bisa menggunakan pendekatan FDR untuk tes dependen (lihat di sini untuk detail ) masalahnya adalah bahwa saya tidak yakin apakah Anda dapat mengatakan di muka jika korelasi Anda semuanya positif.

Dalam R Anda dapat melakukan FDR sederhana dengan p.adjust. Untuk hal-hal yang lebih kompleks, saya akan melihatnya multcomp , tetapi saya tidak melihatnya untuk mencari solusi dalam kasus dependensi.

Semoga berhasil.

Tal Galili
sumber
1
Hai Tal, terima kasih! Bonferroni tampaknya tidak cocok untuk saya - jika salah satu SNP saya bersifat kausal dan yang lain terkait dengannya, harus ada sinyal, dan Bonferroni selalu terlihat terlalu konservatif bagi saya (saya biasanya lebih suka koreksi bertahap Holm). FDR yang Anda tautkan dan hal. Sesuaikan saja tidak mempertimbangkan bukti gabungan (dan FDR masih mengharuskan saya untuk memahami korelasi pengujian saya, pertanyaan awal). multcomp dapat membantu, meskipun pada pandangan pertama sepertinya lebih banyak berurusan dengan beberapa tes dalam satu model, sedangkan saya punya banyak model. Saya akan menggali lebih dalam ...
S. Kolassa - Reinstate Monica
Halo Stephan. Saya mengerti, maaf karena tidak membantu lebih banyak. Semoga berhasil! Tal
Tal Galili
Halo Stephan, saya masih berpikir Anda masih dapat menggunakan metode = BY (untuk Prosedur Benjamini Hochberg Yekuteli) di p.adjust di R, seperti yang ditunjukkan oleh Tal. Jelas, menggunakan Bonferroni bisa menjadi konservatif.
suncoolsu
suncoolsu, saya berpikir bahwa metode ini hanya berfungsi ketika korelasi positif (tidak negatif) antara variabel. Tepuk tangan.
Tal Galili
2

Saya pikir Multivariate Normal Model sedang digunakan untuk memodelkan nilai-p yang berkorelasi dan untuk mendapatkan jenis koreksi pengujian yang tepat. Koreksi Pengujian Cepat dan Akurat Beberapa dan Estimasi Daya untuk Jutaan Penanda Terkait. PLoS Genet 2009 berbicara tentang mereka dan juga memberikan referensi lain. Kedengarannya mirip dengan apa yang Anda bicarakan, tapi saya pikir selain mendapatkan koreksi nilai-p global yang lebih akurat, pengetahuan struktur LD juga harus digunakan untuk menghilangkan positif palsu yang muncul dari marker yang berkorelasi dengan marker kausal.

highBandWidth
sumber
2

Saya mencari solusi untuk masalah yang sama persis. Yang terbaik yang saya temukan adalah Null Unrestricted Bootstrap yang diperkenalkan oleh Foulkes Andrea dalam bukunya Applied Statistics Genetics with R (2009) . Bertolak belakang dengan semua artikel dan buku lain yang dia anggap secara khusus regresi. Selain metode lain ia menyarankan Null Unrestricted Bootstrap, yang cocok di mana orang tidak dapat dengan mudah menghitung residu (seperti dalam kasus saya, di mana saya memodelkan banyak regresi independen (pada dasarnya korelasi sederhana), masing-masing dengan variabel respon yang sama dan snip berbeda). Saya menemukan metode ini juga disebut metode maxT .

> attach(fms)
> Actn3Bin <- > data.frame(actn3_r577x!="TT",actn3_rs540874!="AA",actn3_rs1815739!="TT",actn3_1671064!="GG")
> Mod <- summary(lm(NDRM.CH~.,data=Actn3Bin))
> CoefObs <- as.vector(Mod$coefficients[-1,1]) 
> B <-1000
> TestStatBoot <- matrix(nrow=B,ncol=NSnps)
> for (i in 1:B){
+    SampID <- sample(1:Nobs,size=Nobs, replace=T)
+    Ynew <- NDRM.CH[!MissDat][SampID]
+    Xnew <- Actn3BinC[SampID,]
+    CoefBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,1]
+    SEBoot <- summary(lm(Ynew~.,data=Xnew))$coefficients[-1,2]
+    if (length(CoefBoot)==length(CoefObs)){
+       TestStatBoot[i,] <- (CoefBoot-CoefObs)/SEBoot
+    }
+ }

Setelah kita memiliki semua TestStatBootmatriks (dalam baris kita memiliki replikasi bootstrap, dan dalam kolom kita telah bootstrapT^ statistik) kami temukan untuk yang mana Tkritik. kami mengamati dengan tepat α=0,05 persen lebih signifikan T^ statistik (lebih signifikan berarti bahwa dengan nilai absolut lebih besar dari Tkritik.).

Kami melaporkan saya-komponen model yang signifikan, jika Tsaya^>Tkritik.

Langkah terakhir dapat dilakukan dengan kode ini

p.value<-0.05 # The target alpha threshold
digits<-1000000
library(gtools) # for binsearch

pValueFun<-function(cj)
{
   mean(apply(abs(TestStatBoot)>cj/digits,1,sum)>=1,na.rm=T)
}
ans<-binsearch(pValueFun,c(0.5*digits,100*digits),target=p.value)
p.level<-(1-pnorm(q=ans$where[[1]]/digits))*2 #two-sided.
Adam Ryczkowski
sumber