selectvariants.snplist

#!/bin/bash
ensg=$1
condition=$2
pheno=$(readlink -f $3)
plinkfile=$(readlink -f $4)
matrix=$(readlink -f $5)
gds=$(readlink -f $6)
setfile=$(readlink -f $7)
varidlist=$8
outf=$9
mkdir $outf || exit
cd $outf

echo received variants $varidlist to run on.

#We extract just the variants we want:
grep -w '^'$ensg.$condition $setfile | sed 's/\t/-/;s/\t/:/' | grep -w -f <(echo $varidlist | sed 's/chr//g'| tr ',' '\n') | sed 's/-/\t/;s/:/\t/'>redux.variantset

if [[ ! -s redux.variantset ]]; then
  echo Could not find variants $varidlist for condition $ensg.$condition in $setfile
  sleep 1
  exit
fi

cat redux.variantset

Rscript -e 'library(GMMAT)
library(data.table)
pheno=fread("'$pheno'")
matrix_prefix="'$matrix'"
m=pheno
library(reshape2)
      grm=as.matrix(fread(paste0(matrix_prefix, ".grm.gz"), header=F))
      grm.names=as.character(fread(paste0(matrix_prefix, ".grm.id"), header=F)$V1)
      grm=as.data.table(grm)
      grm[, c("id1", "id2"):=list(grm.names[V1], grm.names[V2])]
      grm.mat=acast(formula=id1~id2, value.var="V4", data=grm)
      grm.mat[is.na(grm.mat)]=t(grm.mat)[is.na(grm.mat)]
      m=m[id %in% unique(grm.names),]
model0=glmmkin(as.formula(paste(colnames(m)[2], "~ 1")), data=m, id="id",  family= gaussian(link = "identity"), kins=grm.mat)
SMMAT.out=SMMAT(null.obj = model0,
            geno.file = "'$gds'",
            group.file = "redux.variantset",
            MAF.range=c(1e-10, 0.05),
            MAF.weights.beta = c(1,1),miss.cutoff=0.01,tests=c("O", "E"),rho=(0:10)/10, use.minor.allele=T,Garbage.Collection = T,meta.file.prefix = "tocondition")

print(t(SMMAT.out))'