zhengxwen/HIBAG

VCF files

Opened this issue · 3 comments

Hello!

Many than for great package!
I wanted to ask how can extract vcf with info score from hla prediction output?

cl <- makeCluster(20)
set.seed(1000)
parseCommandArgs(evaluate=TRUE)
model.obj <- get(load(file1))
model <- hlaModelFromObj(model.obj)
summary(model)
p1=plot(model)
yourgeno <- hlaBED2Geno(bed.fn=bed.file, fam.fn=fam.file, bim.fn=bim.file)
summary(yourgeno)
pred.guess <- hlaPredict(model, yourgeno, match.type="Position")

summary(pred.guess)

hlaAlleleToVCF(hlaAlleleSubset(pred.guess, 1:4),DS=TRUE, verbose=TRUE, outfn=vcf.out)

save(pred.guess, p1, file=file2)

After I execute this script I get a vcf fle with dosages but no infor score.

Many thanks!

HIBAG does not have allele-level confidence scores, instead it provides sample-level confidence scores.
pred.guess$value$prob is sample-level score, see pred.guess$value.

If I understood correctly, then it's better to remove alleles with pred.guess$value$prob < 0.5 ?

Thank you!

pred.guess$value$prob is sample-level (used to remove low-confidence sample prediction).
If you want to remove alleles, consider using the minor allele frequency threshold (e.g., 3 or 5), just like what you do with PLINK if you use PLINK (singleton and doubleton are not considered in single variant tests).