Dose the result have a difference between Rvtest and R SKAT package?
jeong2624 opened this issue · 0 comments
Hi! I'm interested in rare-variant test, such as burden test, SKAT, SKATO and so on.
As I used SKAT test in Rvtest and R SKAT package,
The results of p-value are different by using Rvtest and R SKAT package.
I think the results are same, because of the same test. I don't know why the results are different.
So, what should I do?
example.vcf data are below:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT P1 P2 P3 P4 P5 P6 P7 P8 P9
1 1 . A G 100 . . GT 0/1 0/0 0/0 1/1 0/1 0/0 0/1 0/1 0/0
1 2 . A G 100 . . GT 0/0 0/0 0/0 0/0 0/0 0/0 0/1 0/0 0/1
1 3 . A G 100 . . GT 0/0 0/0 0/0 0/0 0/1 0/0 0/0 0/0 0/0
pheno data are below:
fid iid fatid matid sex y1 y2 y3 y4
P1 P1 0 0 0 1.911 -1.465 -0.817 1
P2 P2 0 0 2 2.146 -2.451 -0.178 2
P3 P3 0 0 2 1.086 -1.194 -0.899 1
P4 P4 0 0 2 0.704 -1.052 -0.237 1
P5 P5 0 0 1 2.512 -3.085 -2.579 1
P6 P6 0 0 2 1.283 -1.294 -0.342 1
P7 P7 0 0 1 2.384 1.070 -1.522 2
P8 P8 0 0 1 3.004 0.680 -1.875 1
P9 P9 0 0 1 0.714 -0.116 -1.604 1
setFile data are below:
1:1-3
The my command in Rvtest
rvtest --noweb --inVcf ./example.vcf.gz --pheno ./pheno --pheno-name y2 --setFile ./setFile --impute hwe --out ./skat/output --kernel skat
The my coomand in R script
library(SKAT)
pheno = read.table("/home/pjw/NGS_Project/Rare_excercise/rvtest_result/excercise/pheno",header = T)
pheno = pheno[,c("y2")]
pheno = as.matrix(pheno)
obj = SKAT_Null_Model(pheno ~ 1, out_type="C", n.Resampling = 10000)
library(vcfR)
geno = read.vcfR("./example.vcf")
gt <- extract.gt(geno)
gt <- extract.gt(geno, element = 'GT', as.numeric = TRUE)
write.table(gt, "./geno_matrix")
geno = read.table("/home/pjw/NGS_Project/Rare_excercise/rvtest_result/excercise/geno_matrix")
geno = as.matrix(t(geno))
SKAT(geno,obj,impute.method = "fixed", weights.beta = c(0.05,1,25))$p.value