I'm looking for a way to simulate phenotypes against a real SNP data source, such as the 1000 Genomes. It must be free for commercial purpose (Eg.: MIT license). Any recommendation? I'm trying to use the GCTA64, but I couldn't get it working. The documentation doesn't help much. At the end of the day, I want: 1 - Simulate Case/Control and/or quantitative phenotype 2 - Link it with a real SNP dataset (eg: 1000 Genomes) 3 - Conduct GWAS analysis using Plink and/or Hail.

Looking Hail's and Plink's tutorials for GWAS, I realised both use simulated phenotypes from real SNP data sources (1000 Genomes and HapMap, respectively), but how they created the datasets are beyond of the scope of the tutorials, thus not reported.

As mentioned before, I've tried gcta64, but no success. Here's what I've tried:

1 - Downloaded 1000 Genome sample from Plink page: Entire dataset as a single .tar.gz (1.12 GB) (A2 allele major, not ref, on chr3 before 15 Oct 2017)

2 - Tried to generate the simulate data by:

./gcta64  --bfile 1kGenomesP1/1kg_phase1_all/1kg_phase1_all --simu-qt  --simu-causal-loci causal.snplist  --simu-hsq 0.5 --simu-rep 3  --keep test
.indi.list --out 1kg_phase1_all

Error: Error: --keep test.indi.list not found.

Question: What should be the files:causal.snplist, test.indi.list, teest? I know the property name suggests what should be, but there's no example of file structure nowhere... Any practical example or tutorial?

Btw - I apologise in advance if this is a too trivial question. I'm quite new at it. Appreciate your patience and help :)


2 Answers 2


I managed to get it done without the gcta64. I have written a Jupyter Notebook with the steps: https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb


You might try ldscsim which was developed by Nik Baya for this purpose.


import hail as hl

mt = hl.import_vcf('path/to/your/data.vcf')
mt.write('data.mt') # store in Hail's native format for speed
mt = hl.read_matrix_table('data.mt')
mt = hl.experimental.ldscsim(
    mt, mt.GT, h2=0.5

ldscsim can generate phenotypes under a number of different models.

For the GWAS tutorial, AFAIK, we just simulated betas for each site, then computed the phenotype from that. If you ask on https://discuss.hail.is and tag Tim Poterba he can probably give you the precise answer.

def logistic(x):
    return 1 / (1 + hl.exp(-x))

mt = mt.annotate_rows(beta = hl.rand_norm(0, 1))
mt = mt.annotate_cols(
    is_case = 0.5 < logistic(hl.agg.sum(beta * mt.GT.n_alt_alleles()))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.