<p>I managed to get it done without the gcta64.
I have written a Jupyter Notebook with the steps: <a href="https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb" rel="nofollow noreferrer">https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb</a></p>

CC BY-SA 4.0
2019-12-01T16:28:56.313Z
2019-12-01T16:28:56.313Z
Bruno Ambrozio
<p>You might try <code>ldscsim</code> which was developed by Nik Baya for this purpose.</p>
<p><a href="https://hail.is/docs/0.2/experimental/ldscsim.html" rel="nofollow noreferrer">https://hail.is/docs/0.2/experimental/ldscsim.html</a></p>
<pre><code>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
)
</code></pre>
<p><code>ldscsim</code> can generate phenotypes under a number of different models.</p>
<p>For the GWAS tutorial, AFAIK, we just simulated betas for each site, then computed the phenotype from that. If you ask on <a href="https://discuss.hail.is" rel="nofollow noreferrer">https://discuss.hail.is</a> and tag Tim Poterba he can probably give you the precise answer.</p>
<pre><code>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 &lt; logistic(hl.agg.sum(beta * mt.GT.n_alt_alleles()))
)
</code></pre>

CC BY-SA 4.0
2023-01-07T01:01:10.520Z
2023-01-07T01:01:10.520Z
Daniel King
<p>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 <a href="http://cnsgenomics.com/software/gcta/#GWASSimulation" rel="nofollow noreferrer">GCTA64</a>, 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.</p>

<p>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.</p>

<p>As mentioned before, I've tried <code>gcta64</code>, but no success. Here's what I've tried:</p>

<p>1 - Downloaded 1000 Genome sample from Plink page:
<a href="https://www.cog-genomics.org/plink/1.9/resources" rel="nofollow noreferrer">Entire dataset as a single .tar.gz (1.12 GB) (A2 allele major, not ref, on chr3 before 15 Oct 2017)</a></p>

<p>2 - Tried to generate the simulate data by:</p>

<pre><code>./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
</code></pre>

<p>Error: <code>Error: --keep test.indi.list not found</code>.</p>

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

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

<snp><gwas><1000genomes>
Simulating phenotype with the 1000 Genomes Project
# Simulating phenotype with the 1000 Genomes Project

(Asked by: username_0 on 2019-11-07)

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 :)

**Tags:** snp, gwas, 1000genomes

---

## 2 answers

### Answer by username_0 on 2019-12-01. Score: 1 | Accepted answer

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>


---

### Answer by username_1 on 2023-01-07. Score: 1

You might try ldscsim which was developed by Nik Baya for this purpose.
https://hail.is/docs/0.2/experimental/ldscsim.html
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()))
)


---

