#add GTEx GWAS

#add COVID-19 host genetics
#and other Lung function ones
#rerun is allowed


#remove the filter snp step, since I rerun the cacluation of l2 for baseline model, and now snps set matching problem should be solved
#use GWAS atlas data, some redundant summary statistics


#for the baseline they get different number of SNPs as I got.
#need to use the overlapping snps, ignore possible change to .M, .M_5_50
#should guarantee all the new annotations share the same set of SNPs


CHRs = as.character(1:22)


sh.head = paste("#!/bin/bash",
                "#$ -S /bin/bash",
                #"#$ -P compbio_lab", 
                "#$ -binding linear:1", 
                "##$ -V",
                "#$ -cwd",
                "#$ -e ./error.$JOB_NAME.$JOB_ID",
                "#$ -o ./outpt.$JOB_NAME.$JOB_ID",
                "#$ -l h_vmem=20g",
                "#$ -l h_rt=10:00:00", 
                "##$ -pe smp 5",
                "source ~/.bashrc",
                "source ~/.my.bashrc",
                "export PATH=$HOME/lhou.compbio/software/anaconda2/bin:$PATH",
                "source activate ldsc",
                sep="\n")


# dir.in.plink = "~/lhou.compbio/data/LDSC/1000G_EUR_Phase3_plink/"
# dir.in.snp = "~/lhou.compbio/data/LDSC/hapmap3_snps/"
file.in.baseline.pref = "~/lhou.compbio/data/LDSC/1000G_Phase3_baselineLD_v2.1_ldscores/baselineLD."
file.in.weight.pref = "~/lhou.compbio/data/LDSC/weights_hm3_no_hla/weights."
file.in.freq.pref = "~/lhou.compbio/data/LDSC/1000G_Phase3_frq/1000G.EUR.QC."



file.in.gwas.list="/broad/compbio/lhou/data/eGTEX-H3K27ac/eGTEx_GWAS_list.0908.2022.txt"
gwas.info =read.table(file.in.gwas.list, sep="\t", row.names = 2, header = T, stringsAsFactors = F)

source2gwas.nm=split(rownames(gwas.info), gwas.info$source)

files.in.gwasSummary = paste0("~/lhou.compbio/data/GWAS/GTEx.v8.GWAS/GTEx.v8.imputedGWAS.withRSID/imputed_", source2gwas.nm$GTEx, ".bed.gz")
names(files.in.gwasSummary)=source2gwas.nm$GTEx
files.in.gwasSummary["Shrine_30804560_FEV1"] = "~/lhou.compbio/data/GWAS/Lung.function_fromXushen/Shrine_30804560_FEV1_meta-analysis.parse.txt.gz"
files.in.gwasSummary["Shrine_30804560_FVC"] = "~/lhou.compbio/data/GWAS/Lung.function_fromXushen/Shrine_30804560_FVC_meta-analysis.parse.txt.gz"
files.in.gwasSummary["Shrine_30804560_PEF"] = "~/lhou.compbio/data/GWAS/Lung.function_fromXushen/Shrine_30804560_PEF_meta-analysis.parse.txt.gz"
files.in.gwasSummary["BD_PGC"] = "~/lhou.compbio/data/GWAS/pgc_bipolar/daner_PGC_BIP32b_mds7a_0416a.hg38.sumstats.gz"

dir.in.annot = paste0("a_annotations_mAREGrpsByTiss/")

# dir.tmp.annot = "~/hptmp/eGTEx-H3K27ac/ldsc_a_2_haQTL_V1.2/" #annot_ovlpWithBaseline
# if(!dir.exists(dir.tmp.annot))
# {
#   dir.create(dir.tmp.annot)
# }
# dir.tmp.annot.baseline = paste0(dir.tmp.annot, "baseline/")
# if(!dir.exists(dir.tmp.annot.baseline))
# {
#   dir.create(dir.tmp.annot.baseline)
# }


dir.ou = paste0("a_2_heritPartition_byLDSC_mAREGrpsByTiss_V1.1_updatedGWAS/")
dir.create(dir.ou, showWarnings = F, recursive = T)

#annotation of interest by regular expression


annot.nms = levels(factor(sapply(dir(dir.in.annot, paste0("annot.gz"), full.names=F), #"haQTL.*\\.1.annot.gz"  "diffPeak.*\\.1.annot.gz""
FUN=function(f.i)
{
 return(gsub(".\\d+.annot.gz", "", f.i, perl=T))
})))


for(pheno in names(files.in.gwasSummary)) 
{
 
#LDSC --h2
  # f.o.sh = paste0("a_2_LDSCPartit_", annot.nm, ".", TISSUE, ".sh")
  # write(sh.head, f.o.sh)
  # 
  print(pheno)
  f.o.sh = paste0("a_2_LDSCPartit_", pheno, ".sh")
  cmds=""
   
  for(annot.nm in annot.nms)#  #[1:8]  [25:29] [9:16]#[17:24]
  {
    #d.tmp.annot.nm = paste0(dir.tmp.annot, annot.nm, "/")
    f.o.res = paste0(dir.ou, pheno,"_", annot.nm)
    f.o.log = paste0(f.o.res, ".log")
    if(file.exists(f.o.log) &&
       length(system(paste0("grep \"Analysis finished\" ", f.o.log), intern = T))==1 )
    {
      print("done")
      next
    }
      
    
    
    cmd = paste0("ldsc.py",
               " --h2 ", files.in.gwasSummary[pheno],
               #" --ref-ld-chr ", dir.tmp.annot.baseline, "baselineLD.,", d.tmp.annot.nm, annot.nm, ".", 
                " --ref-ld-chr ", file.in.baseline.pref, ",", dir.in.annot, annot.nm, ".", 
               " --w-ld-chr ", file.in.weight.pref,
               " --overlap-annot",
               " --frqfile-chr ", file.in.freq.pref,
               " --out ", f.o.res
                )
               
    cmds = paste(cmds, cmd, sep = "\n")
    
    
  }
  if(cmds!="")
  {
    write(sh.head, f.o.sh)
    write(cmds, f.o.sh, append = T)    
    write("echo done!", f.o.sh, append = T)    
    system(paste("qsub " , f.o.sh, sep=""))
  } 
    
}