04.GEMMA_LMM_analysis

Check gene level associations with GEMMA

Gemma uses a LMM with a kinship matrix. That would account for familial structure that PrediXcan doesnt. We’ll run Gemma with the predicted expression and metabolic phenotypes to compare the results.

We first run GEMMA with predicted expression for Nucleus Accumbens Core Tissue generated by PrediXcan. We’ll create GRMS first for all phenotypes

Create gemma format from predicted expression

pred_expr <- fread(data.dir %&% "prediXcan/metabolic_traits/rat_metabolic_Ac_best__predict.txt") %>% select(-c(FID))
pred_expr <- pred_expr %>% pivot_longer(!IID, names_to = "gene", values_to = "count" ) 
pred_expr <- pred_expr %>% pivot_wider(names_from = IID, values_from = count)
pred_expr <- pred_expr %>% mutate(ref = NA, .before = colnames(pred_expr)[1]) %>% mutate(eff = NA, .before = colnames(pred_expr)[1])) 
write_csv(pred_expr, "/gpfs/data/im-lab/nas40t2/natasha/rat_genomics/GEMMA/Ac/LMM/Ac_pred_expr.bimbam", col_names = FALSE)

create phenotype file

pheno <- read_tsv(data.dir %&% "Box_files/processed_obesity_rat_Palmer_phenotypes.tsv")
all_rats <- read_tsv(data.dir %&% "Box_files/all_names.txt", col_names = TRUE)
pheno <- pheno %>% filter(!(rat_rfid %in% all_rats$ID)) %>% select(-c(rat_rfid))

write_tsv(pheno, data.dir %&% "GEMMA/Ac/LMM/processed_obesity_phenotypes_GEMMA.tsv", col_names =FALSE)

create SNP annotation file (we’ll use gene annotation file for Rats)

gtf <- read_tsv(data.dir %&% "Box_files/gtf.txt", col_names = TRUE)
annot <- tibble(gene = pred_expr$gene)
annot <- annot %>% mutate(chr = gtf[na.omit(match(annot$gene, gtf$Gene)), 1]$Chr) %>% mutate(TSS = gtf[na.omit(match(annot$gene, gtf$Gene)), 4]$Start)

#write_tsv(annot, data.dir %&% "GEMMA/Ac/LMM/gene_annotation.tsv", col_names = FALSE)

use gemma/0.98.1 to create GRM kinship matrix

gemma.dir <- data.dir %&% "GEMMA/Ac/LMM/"
setwd(data.dir %&% "GEMMA/Ac/LMM/GRM/")
for(i in 1:ncol(pheno)) {
trait <- colnames(pheno)[i]
runGEMMAgrm <- "gemma -g " %&%  gemma.dir %&% "Ac_pred_expr.bimbam" %&%  " -p " %&% gemma.dir %&% "processed_obesity_phenotypes_GEMMA.tsv" %&% " -n " %&% i %&%  " -notsnp -gk -o " %&% "grm_Ac_" %&% trait
system(runGEMMAgrm)
}

gemma -g [filename] -p [filename] -a [filename] -k [filename] -lmm [num] -o [prefix]

Run Univariate Linear Mixed model

out.dir <- "/gpfs/data/im-lab/nas40t2/natasha/rat_genomics/GEMMA/Ac/LMM/GRM/output/"
setwd(data.dir %&% "GEMMA/Ac/LMM/assoc/")
for(i in 2:ncol(pheno)) {
  trait <- colnames(pheno)[i]
  runGEMMAgrm <- "gemma -g " %&%  gemma.dir %&% "Ac_pred_expr.bimbam" %&%  " -p " %&% gemma.dir %&% "processed_obesity_phenotypes_GEMMA.tsv" %&% " -a " %&% gemma.dir %&% "gene_annotation.tsv"  %&% " -n " %&% i %&%  " -lmm 1 " %&% " -k " %&% out.dir %&%  "grm_Ac_" %&% trait %&% ".cXX.txt -o " %&% "assoc_Ac_" %&% trait
  system(runGEMMAgrm)
}