GEMMA requires data to be presented in a certain way to run successfully. In this document we describe how you can prepare genotype and phenotype data so it passes through GEMMA correctly.
Note that GEMMA uses spaces, tabs and comma’s as valid field separators. Here we use tabs.
In this example we have a genotype spreadsheet containing something like this
“` @type:riset @mat:B @pat:D @het:H @unk:U Chr Locus cM Mb BXD1 BXD2 BXD5 BXD6 BXD8 BXD9 BXD11 … 1 rs31443144 1.50 3.010274 B B D D D … 1 rs6269442 1.50 3.492195 B B D D D … (…) “`
GEMMA, instead, requires a BIMBAM ‘mean genotyped’ formatted file which looks like
“` rs31443144, X, Y, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, … rs6269442, X, Y, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, … (…) “`
where the first three positions refer to SNP id, minor allele, major allele (both unused) and the mean genotypes where B has the value 1, D has the value 0, H has the value 0.5 and U is marked as NA. Note there is no header line.
At this point you may decide to use a spreadsheet or editor to modify these by hand. The problem is that such a procedure is error prone and hard to reproduce. Always choose a programming method if you can! For further information check out these [suggestions](http://kbroman.org/steps2rr/pages/scripts.html).
For most people using R or Python would be the best way to convert the data using a small script. We will add those instructions later, but here I am going to use a few Unix tools and the [bio-table](https://github.com/pjotrp/bioruby-table) tool to manage this conversion from the command line. First, make sure Ruby is available
“`sh ruby -v ruby 2.4.0p0 (2016-12-24 revision 57164) [x86_64-linux] “`
and install bioruby-table
“`sh gem install bio-table “`
On success you should be able to run
“` bio-table “`
In the first step we need to strip the header lines. This can be done with –skip, e.g.
“`sh bio-table –skip 20 BXD.geno > BXD_noheader.geno “`
We’ll keep the header line for some column processing. Next we drop the first column and keep the rest with –columns. For this exercise we are only interested in a subset of BXDs, so
“`sh bio-table –columns Locus,BXD100,BXD101,BXD102,BXD24,BXD32,BXD34,BXD40,BXD43,BXD44,BXD45,BXD48,BXD48a,BXD49,BXD50,BXD51,BXD55,BXD60,BXD61,BXD62,BXD63,BXD65,BXD65b,BXD66,BXD68,BXD69,BXD70,BXD71,BXD73,BXD73a,BXD73b,BXD75,BXD77,BXD78,BXD79,BXD83,BXD86,BXD87,BXD90,BXD98,C57BL/6J,DBA/2J BXD_noheader.geno > BXD_cols.geno “`
selects those columns by name. If a name is a mis-match bio-table will balk. In this case the parents C57BL/6J,DBA/2J are missing and had to be taken out.
It is also possible to use regex’s to select the columns (with –column-filter). Note that bio-table takes tab delimited files by default, to use a different delimiter use something like –in-format split –split-on ‘,’ or a regex. Next we inject columns 1 and 2 by
“`sh bio-table –rewrite ‘field[0] = “X\tY\t” + field[0]’ BXD_cols.geno > BXD_cols2.geno “`
Next rewrite the field contents with something like
“`sh bio-table BXD_cols2.geno –rewrite ‘fields = fields.map { |f| h={“X”=>”X”, “Y”=>”Y”, “D”=> 0, “B”=>1, “H”=>0.5} ; h[f] }’ > BXD_mine.geno “`
And the resulting BIMBAM format looks like
“` rs31443144 X Y B B B B B B B B D B B B D D B B B D D D D D B D B B B D D D D D D D B D B B B rs6269442 X Y 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1 1 1 0 0 0 0 0 1 0 1 1 1 0 0 0 0 0 0 0 1 0 1 1 1 “`
[bio-table](https://github.com/pjotrp/bioruby-table) has an extensive README and many other options, check it out.
To run gemma an annotation file is optional. This file again has no header and looks like
“` rs31443144 3010274 1 rs6269442 3492195 1 rs32285189 3511204 1 “`
where columns refer to the SNP name (same as in the genotype file), the position and the chromosome number. If you don’t supply this information and run gemma with the -debug switch it will complain
“`
“`
This is harmless until you want to run LOCO.
The phenotypes are simply a matrix. Again, without headers. Missing values are marked as NA. Save the file so it looks like [this](https://github.com/genetics-statistics/GEMMA/blob/master/example/mouse_hs1940.pheno.txt). When running GEMMA the column number gets passed in.
The first step is to compute a kinship matrix K. With our files we try
“` gemma -g BXD_mine.geno -p BXD_pheno.csv -gk GEMMA 0.97 (2017/12/19) by Xiang Zhou and team (C) 2012-2017 Reading Files … Segmentation fault “`
Unfortunately there is an immediate problem. Try running with -debug
“` gemma -g BXD_mine.geno -p BXD_pheno.csv -gk -debug GEMMA 0.97 (2017/12/19) by Xiang Zhou and team (C) 2012-2017 Reading Files …
strtok failed in ReadFile_geno in src/io.cpp at line 706 “`
and we get a slightly more informative error. There is a parsing problem. GEMMA is not very good at data format checking so we need to do some checks. First see if the number of individuals matches in the phenotype file and genotype file. In this case the phenotype file was too large. In still included a header line and the parent phenotypes.
After removing those it was fine.
Gemma contains helpful information, particularly
“`sh gemma -h 2 “`
where we can select the column in the phenotype file with -n. This is important because (at this point) K accounts for missing values:
“`sh gemma -g BXD_mine.geno -p BXD_pheno.csv -gk -n 2 “`
by default it writes K to output/result.cXX.txt.
Gemma also has a leave one chromosoe out (LOCO) option, but it is best to use [gemma-wrapper](https://github.com/genetics-statistics/gemma-wrapper) for that because gemma-wrapper iterates through all chromosomes. More on that below.
Now we have K, let’s run an LMM:
Running an LMM on the phenotype in column 2 (they are numbered 1,2,…)
“` gemma -g BXD_mine.geno \ -p BXD_pheno.csv \ -a example/BXD_snps.txt \ -k output/result.cXX.txt \ -lmm 1 -maf 0.1 -n 2 \ -o n2_assoc.txt “`
gets a result file in output/result.assoc.txt
For LOCO we use gemma-wrapper because it facilitates running through the chromosomes. Essentially the command is the same (after – gets passed to gemma):
“` gemma-wrapper –json \ –loco 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,X – \ -g BXD_mine.geno \ -p BXD_pheno.csv \ -a example/BXD_snps.txt \ -gk \ -n 2 \ -debug > K.json “`
Note that the annotation file now is required to find the chromosomes. Gemma-wrapper computes the K’s for every chromosome and the filenames are passed into K.json. A file which can be fed to the LMM:
“` gemma-wrapper –loco –input K.json – \ -g BXD_mine.geno \ -p BXD_pheno.csv \ -a example/BXD_snps.txt \ -lmm 1 -maf 0.1 -n 2 \ -debug > GWA.json “`
Inside the resulting GWA.json file there is a list of results. Check out the paths and compile the result into one file with, for example,
“` bio-table /tmp/aa2a1c1a67fe2289d6a23afcc025818402f97521.*.assoc.txt.assoc.txt > LOCO_n2_assoc.txt “`
To get an idea of what is a significant hit we can also use run gemma 1000x after shuffling the phenotypes. gemma-wrapper also has an option for that. First create K
“` gemma-wrapper –json – \ -g BXD_mine.geno \ -p BXD_pheno.csv \ -a example/BXD_snps.txt \ -gk \ -n 2 \ -debug > K.json “`
Now run once with
“` gemma-wrapper –json –input K.json – \ -g BXD_mine.geno \ -p BXD_pheno.csv \ -a example/BXD_snps.txt \ -lmm 2 -maf 0.1 -n 2 \ -debug > GWA.json “`
GWA.json should point to the result set.
To run the permutations add one option and move the -p option to –phenotypes before the double dash
“` gemma-wrapper –permutate 1000 –phenotypes BXD_pheno.csv –input K.json – \ -g BXD_mine.geno \ -a example/BXD_snps.txt \ -lmm 1 -maf 0.1 -n 2 \ -debug > GWA.json “`
gemma-wrapper prints out the 95 and 67 percentiles where the first may be considered the ‘significant’ threshold and the latter is the ‘suggestive’ threshold. For example
“` [“95 percentile (significant) “, 3.382905e-05, 4.5] [“67 percentile (suggestive) “, 0.0003651852, 3.4] “`
where 4.5 is the LOD, i.e., -log10(3.382905e-05)
Once there is a list of associations out of GEMMA you may use the gemma-annotate tool as described [here](https://github.com/pjotrp/gemma-annotate).