Prj01_RL_pheno_qtl2_DO_11092020_69k.Rmd

---
title: "Prj01_RL_pheno_qtl2_DO_11092020_69k"
author: "Hao He"
output:
  workflowr::wflow_html:
    toc: true
    toc_float: true
    code_folding: "hide"
---

```{r knitr-opts, include=FALSE}
knitr::opts_chunk$set(comment = "#",collapse = TRUE)
```

**Last update:** `r Sys.Date()`

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

### loading libraries
```{r library, message=FALSE, warning=FALSE}
library(ggplot2)
library(gridExtra)
library(GGally)
library(parallel)
library(qtl2)
library(parallel)
library(survival)
library(regress)
library(abind)
library(openxlsx)
library(tidyverse)
library(SimDesign)
library(MASS)
library(EnvStats)


rz.transform <- function(y) {
  rankY=rank(y, ties.method="average", na.last="keep")
  rzT=qnorm(rankY/(length(na.exclude(rankY))+1))
  return(rzT)
}
# heritability
# function to simulate data to estimate heritability intervals
h2_sim <- function(trait, K, covar, nsim=100){
  # trait, K and covar are matrices with row.names in qtl2-style
  # covar adjustment applied to data but not in simulated data

  # reduce mice to common set with kinship and trait values
  # mice.names <- intersect(row.names(K), row.names(trait))
  # K <- K[mice.names, mice.names]
  # trait <- as.matrix(trait[mice.names,],ncol=1)
  N <- dim(trait)[1]

  # estimate heritability
  h2 <- as.numeric(est_herit(trait, K, covar))

  # simulate from kinship matrix and re-estimate
  h2_sim <- rep(0,nsim)
  for(i in 1:nsim){
    trait_sim <- t(rmvnorm(n = 1, mean = rep(0,N), sigma = h2*2*K + (1-h2)*diag(1,N)))
    rownames(trait_sim) <- rownames(trait)
    h2_sim[i] <- as.numeric(est_herit(trait_sim, K))
  }

  # return estimated h2 and simulations
  list(N=N, h2=h2, sim=h2_sim)
}
setwd("/projects/csna/csna_workflow/")
```

### Read phenotype data
```{r Read phenotype data, message=FALSE, fig.height=8, fig.width=6}
#Prj01_RL_pheno_qttl2_DO_11092020
pheno <- readr::read_csv("data/pheno/Prj01_RL_pheno_qttl2_DO_11092020.csv", na = c("N/A", "NA"))
pheno <- pheno %>%
  filter(apply(pheno, 1, function(x) sum(is.na(x))) < 5) %>% # remove rows with 5 NAs
  mutate(
    Acq.Anticipatory.Correct.Responses.nooutlier = case_when(
      (Acq.Anticipatory.Correct.Responses > 5 | is.na(Acq.Anticipatory.Correct.Responses)) ~ NA_real_,
      TRUE ~ Acq.Anticipatory.Correct.Responses),
    Rev.Anticipatory.Incorrect.Responses.nooutlier = case_when(
      (Rev.Anticipatory.Incorrect.Responses > 6) ~ NA_real_,
      TRUE ~ Rev.Anticipatory.Incorrect.Responses
    )) %>% #boxcoxtransformation
  mutate(across(c(Rev.Anticipatory.Incorrect.Responses), ~ EnvStats::boxcoxTransform(.x, lambda = EnvStats::boxcox(.x, optimize = TRUE)$lambda),
                .names ="{.col}.boxcox.trans")) %>% #omit values >6 for Acq.Anticipatory.Correct.Responses and Rev.Anticipatory.Incorrect.Responses variables.
  mutate(across(Acq.Total.Trials:Rev.Anticipatory.Incorrect.Responses.nooutlier, ~ rz.transform(.x),
                .names ="{.col}.rz"))

#load gm data
load("data/Jackson_Lab_12_batches/gm_DO3173_qc_newid.RData")
#subset
overlap.id <- intersect(pheno$mouseID, ind_ids(gm_after_qc))

load("data/Jackson_Lab_12_batches/qc_info.RData")
#merge pheno with covar
pheno <- pheno %>%
  left_join(qc_info) %>%
  filter(mouseID %in% overlap.id) %>%
  filter(bad.sample == FALSE) %>%
  mutate(sex = case_when(
    sex == "M" ~ 1,
    sex == "F" ~ 0
  )) %>%
  mutate_at(vars(sex,
                 ngen),
            list(factor))

#pgmap
load("/projects/csna/csna_workflow/data/69k_grid_pgmap.RData")
query_variants <- create_variant_query_func("/projects/csna/csna_workflow/data/cc_variants.sqlite")
query_genes <- create_gene_query_func("/projects/csna/csna_workflow/data/mouse_genes_mgi.sqlite")

```

### Plot qtl mapping on 69k
Plot by using its output Prj01_RL_11092020_cutoff6.qtlout.69k.RData
```{r Plot qtl mapping on 69k, message=FALSE, warning=FALSE}
load("output/Prj01_RL_11092020_cutoff6.qtlout.69k.RData")

#genome-wide plot
for(i in names(m2.qtl.out)){
  par(mar=c(5.1, 4.1, 1.1, 1.1))
  ymx <- maxlod(m2.qtl.out[[i]]) # overall maximum LOD score
  plot(m2.qtl.out[[i]], map=pmap, lodcolumn=1, col="slateblue", ylim=c(0, 12))
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[1]], col="red")
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[2]], col="red")
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[3]], col="red")
  title(main = paste0(i))
}

pdf(file = paste0("output/RL_prj/Prj01_RL_11092020_cutoff6.qtlout.pdf"), width = 16, height =8)
#save genome-wide plo
for(i in names(m2.qtl.out)){
  par(mar=c(5.1, 4.1, 1.1, 1.1))
  ymx <- maxlod(m2.qtl.out[[i]]) # overall maximum LOD score
  plot(m2.qtl.out[[i]], map=pmap, lodcolumn=1, col="slateblue", ylim=c(0, 10))
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[1]], col="red")
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[2]], col="red")
  abline(h=summary(m2.permu[[i]], alpha=c(0.10, 0.05, 0.01))[[3]], col="red")
  title(main = paste0(i))
}
dev.off()
```

### Coeff plot of qtl mapping
```{r Coeff plot of qtl mapping, message=FALSE, warning=FALSE}
#peaks coeff plot
for(i in names(m2.qtl.out)){
  print(i)
  peaks <- find_peaks(m2.qtl.out[[i]], map=pmap, threshold=6, drop=1.5)
  if(nrow(peaks) > 0) {
    print(peaks)
    for(p in 1:dim(peaks)[1]) {
    print(p)
    chr <-peaks[p,3]
    #coeff plot
    par(mar=c(4.1, 4.1, 0.6, 0.6))
    plot_coefCC(m2.sigqtl.coef[[i]][[p]], pmap[chr], scan1_output=m2.qtl.out[[i]], bgcolor="gray95", legend=NULL)
    plot_coefCC(m2.sigqtl.blup[[i]][[p]], pmap[chr], scan1_output=m2.qtl.out[[i]], bgcolor="gray95", legend=NULL)
    #plot(m2.sigqtl.snps[[i]][[p]]$lod, m2.sigqtl.snps[[i]][[p]]$snpinfo, drop_hilit=1.5, genes=m2.sigqtl.genes[[i]][[p]])
    }
  }
}

#save peaks coeff plot
for(i in names(m2.qtl.out)){
  print(i)
  peaks <- find_peaks(m2.qtl.out[[i]], map=pmap, threshold=6, drop=1.5)
  if(nrow(peaks) > 0){
      fname <- paste("output/RL_prj/Prj01_RL_11092020_cutoff6.qtlout_", str_replace_all(i, "[[:punct:]]", "") ,"_coefplot.pdf",sep="")
      pdf(file = fname, width = 16, height =8)
      for(p in 1:dim(peaks)[1]) {
        print(p)
        chr <-peaks[p,3]
        #coeff plot
        par(mar=c(4.1, 4.1, 0.6, 0.6))
        plot_coefCC(m2.sigqtl.coef[[i]][[p]], pmap[chr], scan1_output=m2.qtl.out[[i]], bgcolor="gray95", legend=NULL)
        plot_coefCC(m2.sigqtl.blup[[i]][[p]], pmap[chr], scan1_output=m2.qtl.out[[i]], bgcolor="gray95", legend=NULL)
        #plot(m2.sigqtl.snps[[i]][[p]]$lod, m2.sigqtl.snps[[i]][[p]]$snpinfo, drop_hilit=1.5, genes=m2.sigqtl.genes[[i]][[p]])
        }
      dev.off()
  }
}

#save peaks coeff blup plot
for(i in names(m2.qtl.out)){
  print(i)
  peaks <- find_peaks(m2.qtl.out[[i]], map=pmap, threshold=6, drop=1.5)
  if(nrow(peaks) > 0) {
      fname <- paste("output/RL_prj/Prj01_RL_11092020_cutoff6.qtlout_", str_replace_all(i, "[[:punct:]]", "") ,"_coefplot_blup.pdf",sep="")
      pdf(file = fname, width = 16, height =8)
      for(p in 1:dim(peaks)[1]) {
        print(p)
        chr <-peaks[p,3]
        #coeff plot
        par(mar=c(4.1, 4.1, 0.6, 0.6))
        plot_coefCC(m2.sigqtl.blup[[i]][[p]], pmap[chr], scan1_output=m2.qtl.out[[i]], bgcolor="gray95", legend=NULL)
        #plot(m2.sigqtl.snps[[i]][[p]]$lod, m2.sigqtl.snps[[i]][[p]]$snpinfo, drop_hilit=1.5, genes=m2.sigqtl.genes[[i]][[p]])
        }
      dev.off()
    }
}

```

### heritability by GCTA
```{r heritability by GCTA, eval=FALSE}
#gcta id
gcta.id <- read.table("/projects/csna/csna_workflow/data/GCTA/12_batches.fam", header = F, sep = " ")
#overlap
pheno.overlap <- pheno[pheno$mouseID %in% intersect(gcta.id$V1, pheno$mouseID),]
pheno.overlap$sex <- as.numeric(pheno.overlap$sex)
pheno.overlap$sex[pheno.overlap$sex == 0] <- 2 #female
dim(pheno.overlap)

#subset id
idlist <- cbind(pheno.overlap$mouseID, pheno.overlap$mouseID)
write.table(idlist, file = "/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/update_id.txt", sep = "\t", col.names = FALSE, row.names = FALSE, quote = FALSE)

#update sex
sex <- data.frame(id1 = pheno.overlap$mouseID,
                  id2 = pheno.overlap$mouseID,
                  sex = as.character(pheno.overlap$sex))
write.table(sex, file = "/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/update_sex.txt",sep = "\t", col.names = FALSE, row.names = FALSE, quote = FALSE)

#updated covariate
covar <- data.frame(id1 = pheno.overlap$mouseID,
                    id2 = pheno.overlap$mouseID,
                    sex = as.character(pheno.overlap$sex)
                    )
write.table(covar, file = "/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/update_covar.txt",sep = "\t", col.names = FALSE, row.names = FALSE, quote = FALSE)

#update pheno
update.pheno <- cbind(data.frame(FID = pheno.overlap$mouseID,
                                 IID = pheno.overlap$mouseID),
                      pheno.overlap[,4:17]
)
write.table(update.pheno, file = "/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/update.pheno.txt",sep = "\t", col.names = FALSE, row.names = FALSE, quote = FALSE)


#subset id and update sex
system(paste0("cd /projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/; /projects/csna/csna_workflow/data/GCTA/plink -bfile ", 
       "/projects/csna/csna_workflow/data/GCTA/","12_batches ",
       "--keep update_id.txt ",
       "--update-sex update_sex.txt ",
       "--make-bed --out ",
       "Prj01_RL_pheno_qttl2_DO_11092020"))

# Estimate the GRM
system(paste0("cd /projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/; /projects/csna/csna_workflow/data/GCTA/gcta64  --bfile ",
              "Prj01_RL_pheno_qttl2_DO_11092020",  " --make-grm  --out ",
              "Prj01_RL_pheno_qttl2_DO_11092020"))

for(i in 1:length(colnames(pheno.overlap)[4:17])){
  system(
    paste0("cd /projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/; /projects/csna/csna_workflow/data/GCTA/gcta64 --reml  --grm ",
           "Prj01_RL_pheno_qttl2_DO_11092020",  " --pheno update.pheno.txt --covar update_covar.txt --mpheno ",
           i," --out ",colnames(pheno.overlap)[4:17][i],".out")
  )
}

#plot heritability by GCTA
hsq.gcta <- list()
phe.name <- colnames(pheno.overlap)[4:17]
for (i in phe.name){
  hsq.gcta[[i]] <- read.table(file = paste0("/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/",i,".out.hsq"),sep = "\t", header = FALSE,fill = TRUE, stringsAsFactors = FALSE)
}
h <- data.frame(
  Phenotype = phe.name,
  Heritability = as.numeric(as.vector(unlist(lapply(hsq.gcta, FUN = function(x){x[5,2]})))),
  SE           = as.numeric(as.vector(unlist(lapply(hsq.gcta, FUN = function(x){x[5,3]})))),
  Sample_size  = as.numeric(as.vector(unlist(lapply(hsq.gcta, FUN = function(x){x[11,2]})))),
#  Domain = sub("\\..*", "", phe.name),
  stringsAsFactors = FALSE)
write.csv(h, file = "/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/h2.csv", row.names = F, quote = F)

#histgram
h$Heritability <- round(h$Heritability,2)
pdf(file = paste0("/projects/csna/csna_workflow/data/GCTA/Prj01_RL_pheno_qttl2_DO_11092020/","Prj01_RL_pheno_qttl2_DO_11092020__heritability_by_GCTA.pdf"), height = 10, width = 10)
p<-ggplot(data=h, aes(x=Phenotype, y=Heritability)) + #, fill=Domain, color = Domain)) +
  geom_bar(stat="identity", fill = "blue", color = "blue", show.legend = FALSE) +
  scale_y_continuous(breaks=seq(0.0, 1.0, 0.1)) +
  geom_text(aes(label = Heritability, y = Heritability + 0.005), position = position_dodge(0.9),vjust = 0) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))
p
dev.off()

#plot heritability by qtl2 array --------------
herit <- data.frame(Phenotype = names(unlist(qtl.hsq)),
                    Heritability = round(unlist(qtl.hsq),3))
herit <- herit %>%
  arrange(desc(Heritability))
herit$Phenotype <- factor(herit$Phenotype, levels = herit$Phenotype)
#histgram
p2 <- ggplot(data=herit, aes(x=Phenotype, y=Heritability)) + #, fill=Domain, color = Domain)) +
  geom_bar(stat="identity", fill = "blue", color = "blue", show.legend = FALSE) +
  scale_y_continuous(breaks=seq(0.0, 1.0, 0.1)) +
  geom_text(aes(label = Heritability, y = Heritability + 0.005), position = position_dodge(0.9),vjust = 0) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  ggtitle("Heritability by qtl2 array")
p2

```

### heritability by qtl2 and the bootstrap estimates
```{r}
apr <- readRDS("output/RL_prj/apr.RDS")
#kinship
k_overall = calc_kinship(probs = apr, type = "overall", use_allele_probs = TRUE, cores = 20)

# #heritability -----------------------------------------------------------
# loop through traits
#pheno list
pheno <- as.data.frame(pheno)
pheno.names <- colnames(pheno)[4:18]
pheno.herit <- data.frame(Name = pheno.names, h2=0, N=0, ll=0, ul=0)
covar = model.matrix(~sex, data = pheno)[,-1, drop=F]
rownames(covar) <- pheno$mouseID
colnames(covar)[1] = "sex"
#loop
for(i in 1:length(pheno.names)){
  j = pheno.names[[i]]
  print(j)
  phe = pheno[,j,drop=F]
  rownames(phe) = pheno$mouseID
  tmp <- h2_sim(phe, k_overall, covar, 100)
  pheno.herit$N[i] <- tmp$N
  pheno.herit$h2[i] <- tmp$h2
  pheno.herit$ll[i] <- quantile(tmp$sim, 0.1)
  pheno.herit$ul[i] <- quantile(tmp$sim, 0.9)
}
write.csv(pheno.herit, file = "output/RL_prj/Prj01_RL_11092020_cutoff6.pheno.herit.csv",
          row.names = FALSE, quote = F)
```