pi1_calc.R

library(data.table)
library(dplyr)
library(ggplot2)
library(qvalue)
"%&%" = function(a,b) paste(a,b,sep="")

#List METS and MESA pops
discovery_cohorts<-c("METS_FDR0.05_PC0_PF0.txt.gz","METS_FDR0.05_PC0_PF10.txt.gz","METS_FDR0.05_PC0_PF20.txt.gz","METS_FDR0.05_PC10_PF0.txt.gz","METS_FDR0.05_PC10_PF10.txt.gz","METS_FDR0.05_PC10_PF20.txt.gz","METS_FDR0.05_PC3_PF0.txt.gz","METS_FDR0.05_PC3_PF10.txt.gz","METS_FDR0.05_PC3_PF20.txt.gz")
replication_cohorts<-as.matrix(read.table(file="~/mets_analysis/meqtl/replication_pops/pop_list"))
discovery_dir<-"~/mets_analysis/meqtl/combined_pop/"
replication_dir<-"~/mets_analysis/meqtl/replication_pops/"

#initialize an empty pi1 Matrix
pi1_matrix<-matrix(NA, nrow = length(discovery_cohorts), ncol = length(replication_cohorts))
rownames(pi1_matrix)<-discovery_cohorts
colnames(pi1_matrix)<-replication_cohorts

for (i in c(1:length(discovery_cohorts))){
  #read in discovery files
  discovery_file<-read.table(file = discovery_dir %&% discovery_cohorts[i], stringsAsFactors = F)
  colnames(discovery_file)<-c("snps","gene","statistic","pvalue","FDR","beta")
  row<-discovery_cohorts[i]
  
  for (j in c(1:length(replication_cohorts))){
    #read in replication files
    replication_file<-read.table(file = replication_dir %&% replication_cohorts[j], header = T)
    colnames(replication_file)<-c("snps","gene","statistic","pvalue","FDR","beta")
    #select the  snp-gene pairs that replicate
    replicated<-inner_join(discovery_file, replication_file, by=c("snps","gene"))
    str(replicated)
    #grab statistics for pi1 calculation
    over<-dim(replicated)
    replicated_pvals<-replicated$pvalue.y
    str(replicated_pvals)
    #calculate pi1
    qobj_replicated<-qvalue(p = replicated_pvals)
    pi1<- 1 - qobj_replicated$pi0
    #truncate pi1
    pi12<- signif(pi1,4)
    col<-replication_cohorts[j]
    #add to matrix
    pi1_matrix[row,col] <- pi12
    cat("finished combination " %&% i %&% " " %&% j)
  }
}

write.table(pi1_matrix, file="/home/ryan/mets_analysis/meqtl/replication_pop",quote=F, row.names = F, sep = '\t')