reanalyse.single.burden

#!/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 condition on.

#plink --bfile $plinkfile --extract <(echo $varidlist | tr ',' '\n') --recode A --out tocondition
grep -w '^'$ensg.$condition $setfile>redux.variantset
Rscript -e 'library(GMMAT)
library(data.table)
pheno=fread("'$pheno'")
m=pheno
matrix_prefix="'$matrix'"
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% 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(SMMAT.out)'