---
title: "pcadjust_exprvals.Rmd"
author: "Sahin Naqvi"
date: "1/15/2019"
output: html_document
---

Read in unadjusted human GTEx counts and TPM, as well as sample metadata
```{r}
data_counts_unadj=as.matrix(read.delim(gzfile("../exprvals/gtex.filt.salmon.tximport.unadj.counts.txt.gz"),header=TRUE,stringsAsFactors = FALSE,row.names = 1))
data_unadj=as.matrix(read.delim(gzfile("../exprvals/gtex.filt.salmon.tximport.unadj.tpm.txt.gz"),header=TRUE,stringsAsFactors = FALSE,row.names = 1))

meta = read.table("../metadata/human.metadata.txt",header=TRUE,stringsAsFactors = FALSE)
meta$SUBJID = sapply(strsplit(meta$SAMPID,"-[0-9]{4}-"),"[[",1)

histeval = readRDS("../metadata/histeval.rds")
```

First pass linear modeling of sex in each tissue to perform PCA 
```{r}
require(limma)
require(edgeR)
fpkm_cut = 1
tissue_resid_pca = list() 

for (tissue in unique(meta$Tissue)){
  print(tissue)
  
  # subset metadata table and TPM values to calculate male and female mean
  meta_subset = subset(meta,Tissue == tissue)
  data_subset = data[,(colnames(data) %in% meta_subset$RUN)]
  
  meta_subset_male = subset(meta_subset,(Tissue == tissue)&(SEX=="male"))
  data_subset_male = data_subset[,(colnames(data_subset) %in% meta_subset_male$RUN)]
  
  meta_subset_female = subset(meta_subset,(Tissue == tissue)&(SEX=="female"))
  data_subset_female = data_subset[,(colnames(data_subset) %in% meta_subset_female$RUN)]
  
  means = log2(apply(2^data_subset-1, 1, mean)+1)
  means_m = log2(apply(2^data_subset_male-1, 1, mean)+1)
  means_f = log2(apply(2^data_subset_female-1, 1, mean)+1)
  vars = apply(2^data_subset-1, 1, var)
  
  # choose genes with expression greater than cutoff in either males or females
  passfilt = which(((means_m > fpkm_cut)|(means_f > fpkm_cut))&(vars > 0))
  
  # create model matrix and perform differential expression
  varnames = c("SEX")
  coeffnames = c("SEXmale")
  tiss_mod = model.matrix(formula(paste(c("~", paste(c("1", varnames), collapse=" + ")), collapse=" ")) ,data=meta_subset)
  
  data_subset = data_counts[,meta_subset$RUN] 
  tmm = DGEList(data_subset[passfilt,])
  tmm = calcNormFactors(tmm,method="TMM")
  v = voom(tmm,design = tiss_mod)
  fit = eBayes(lmFit(v,design = tiss_mod))
  
  # perform PCA on residuals
  tissue_resid_pca[[tissue]] = prcomp(residuals(fit,v$E))

}

```

choose number of PCs for each tissue
```{r}
source("../scripts/choosePCs.R")
tissue_npcs = list()
for (tissue in tissue_list) {tissue_npcs[[tissue]] = choosePCs(tissue_resid_pca[[tissue]])}
```


Second pass linear modeling of sex in each tissue, adjusting for chosen number of PCs (and histological features for adipose and thyroid)
```{r}
data_adj_counts = data_counts

for (tissue in unique(meta$Tissue)){
  print(tissue)
  
  meta_subset = subset(meta,Tissue == tissue)
  data_subset = data[,(colnames(data) %in% meta_subset$RUN)]
  
  meta_subset_male = subset(meta_subset,(Tissue == tissue)&(SEX=="male"))
  data_subset_male = data_subset[,(colnames(data_subset) %in% meta_subset_male$RUN)]
  
  meta_subset_female = subset(meta_subset,(Tissue == tissue)&(SEX=="female"))
  data_subset_female = data_subset[,(colnames(data_subset) %in% meta_subset_female$RUN)]
  
  means = log2(apply(2^data_subset-1, 1, mean)+1)
  means_m = log2(apply(2^data_subset_male-1, 1, mean)+1)
  means_f = log2(apply(2^data_subset_female-1, 1, mean)+1)
  vars = apply(2^data_subset-1, 1, var)
  passfilt = which(((means_m > fpkm_cut)|(means_f > fpkm_cut))&(vars > 0))
  
  varnames = c("SEX")
  coeffnames = c("SEXmale")
  
  # add pcs to meta table
  varnames = c(varnames,paste("PC",1:tissue_npcs[[tissue]],sep=""))
  coeffnames = c(coeffnames,paste("PC",1:tissue_npcs[[tissue]],sep=""))
  temp_pcs = data.frame(tissue_resid_pca[[tissue]]$rotation[,1:tissue_npcs[[tissue]]])
  colnames(temp_pcs) = paste("PC",1:tissue_npcs[[tissue]],sep="")
  meta_subset = merge(meta_subset,temp_pcs,by.x="RUN",by.y=0) 
  
  # if adipose or thyroid, add histology evaluations to meta table
  if (tissue=="Adipose_Visceral"){
     meta_subset = merge(meta_subset,histeval[["Adipose_Visceral"]][,c("SUBJID","PERC_VESS_CONN")],by="SUBJID",by.y=1)
     varnames = c(varnames,"PERC_VESS_CONN")
     coeffnames = c(coeffnames,"PERC_VESS_CONN")
   }
  if (tissue=="Thyroid"){
     meta_subset = merge(meta_subset,histeval[["Thyroid"]][,c("SUBJID","INTERSTIT_FIBROSIS","INFLAMMATION")],by="SUBJID",by.y=1)
     varnames = c(varnames,"INTERSTIT_FIBROSIS")
     coeffnames = c(coeffnames,"INTERSTIT_FIBROSIS")
   }

  
  tiss_mod = model.matrix(formula(paste(c("~", paste(c("1", varnames), collapse=" + ")), collapse=" ")) ,data=meta_subset)
  
  data_subset = data_counts[,meta_subset$RUN] 
  tmm = DGEList(data_subset[passfilt,])
  tmm = calcNormFactors(tmm,method="TMM")
  v = voom(tmm,design = tiss_mod)
  fit = eBayes(lmFit(v,design = tiss_mod))
  
  # adjust counts for PCs and histological criteria
  data_subset_adj = v$E
  colnames(data_subset_adj) = colnames(data_subset)
  data_subset_adj = data_subset_adj - (fit$coefficients[,3:ncol(fit$coefficients)] %*% t(tiss_mod[,3:ncol(tiss_mod)]))
  data_adj_counts[rownames(data_subset_adj),colnames(data_subset_adj)] = t(t(2^data_subset_adj-0.5) * (apply(data_counts[rownames(data_subset_adj),colnames(data_subset_adj)],2,sum)/1000000))
  
}

# some counts after adjustment are negative -- set to 0
data_adj_counts[which(data_adj_counts < 0)] = 0

write.table(data_adj_counts,file="../exprvals/gtex.filt.salmon.tximport.adj.counts.txt.gz",col.names = TRUE,row.names = TRUE,sep="\t",quote=FALSE)
```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
