Overlap with Full GWAS Catelog

QTLs

Full GWAS catelog from the table browser. There are 56248699 lines in this file:

/project2/gilad/briana/genome_anotation_data/hg19.GWASCatelog.allsnps

First I want to subset this to a bed file to use. I also want to subset only to SNPs.

Columns: bin chrom chromStart chromEnd name score strand refNCBI refUCSC observed molType class valid avHet avHetSE func locType weight exceptions submitterCount submitters alleleFreqCount alleles alleleNs alleleFreqs bitfields

sed 's/^chr//' /project2/gilad/briana/genome_anotation_data/hg19.GWASCatelog.allsnps  > /project2/gilad/briana/genome_anotation_data/hg19.GWASCatelog.allsnps.bed

overlapSNPsGWAS_fixed.py

import pybedtools as pybedtools

def main(infile, outfile):
    gwas_file=open("/project2/gilad/briana/genome_anotation_data/hg19.GWASCatelog.allsnps.bed","r")
    gwas=pybedtools.BedTool(gwas_file)
    snps_file=open(infile, "r")
    snps=pybedtools.BedTool(snps_file)
    snpOverGWAS=snps.intersect(gwas, wa=True,wb=True)
    snpOverGWAS.saveas(outfile)

if __name__ == "__main__":
    import sys
    import pybedtools
    infile=sys.argv[1]
    outfile=sys.argv[2]
    main(infile, outfile) 

run_overlapSNPsGWASFixed_proc.sh

#!/bin/bash

#SBATCH --job-name=run_overlapSNPsGWASFixed_proc
#SBATCH --account=pi-yangili1
#SBATCH --time=5:00:00
#SBATCH --output=run_overlapSNPsGWASFixed_proc.out
#SBATCH --error=run_overlapSNPsGWASFixed_proc.err
#SBATCH --partition=broadwl
#SBATCH --mem=10G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env


python overlapSNPsGWAS_fixed.py  "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllOverlapSnps.bed" "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/GWASoverlapped_AllOverlapSnps.bed"

This analysis gives 9k overlaps with 7k uniq snps. of (53135726 uniq snps)

This makes more sense.

I want to only look at relevent GWAS for LCLs

QTLOverlap=fread("../data/GWAS_overlap/GWASoverlapped_AllOverlapSnps.bed", header=F, col.names = c("chromSnp", "startSnp", "endSnp", "Set", "chromGwas", "startGWAS", "endGWAS", "rsID", "score", "strand"))

QTL_overlap_Total=QTLOverlap %>% filter(grepl("Total",Set))

QTL_overlap_Nuclear=QTLOverlap %>% filter(grepl("Nuclear",Set))

MAtched snps

I will also need to compare to random snps. (i can use my matched snps/find those in LD)

• /project2/gilad/briana/threeprimeseq/data/MatchedSnp/Nuclear_matched_snps_sort.bed
• /project2/gilad/briana/threeprimeseq/data/MatchedSnp/Total_matched_snps_sort.bed

Switch to format: snp_10_3154947

fixMatchedFormat.py

def fix_format(inbed,outf):
    bed=open(inbed, "r")
    outF=open(outf, "w")
    for ln in bed:
        chrom, start, end = ln.split()
        outF.write("snp_%s_%s\n"%(chrom, end))
    outF.close()


fix_format("/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Nuclear_matched_snps_sort.bed", "/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Nuclear_matched_snps_GEUFormat.txt") 

fix_format("/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Total_matched_snps_sort.bed", "/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Total_matched_snps_GEUFormat.txt") 

subset_plink4Matched_proc.py

def main(genFile, qtlFile, outFile):
  #convert snp file to a list: 
  def file_to_list(file):
    snp_list=[]
    for ln in file:
      snp=ln.strip()
      snp_list.append(snp)
    return(snp_list)

  gen=open(genFile,"r")
  fout=open(outFile, "w")
  qtls=open(qtlFile, "r")
  qtl_list=file_to_list(qtls)
  for ln in gen:
      snp=ln.split()[2]
      if snp in qtl_list:
          fout.write(ln)
  fout.close()
    

if __name__ == "__main__":
    import sys
    chrom=sys.argv[1]
    fraction=sys.argv[2]
    genFile = "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap/geu_plinkYRI_LDchr%s.ld"%(chrom)
    outFile= "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/%sApaMatch_LD/chr%s.%sMatch.LD.geno.ld"%(fraction,chrom,fraction)
    qtlFile= "/project2/gilad/briana/threeprimeseq/data/MatchedSnp/%s_matched_snps_GEUFormat.txt"%(fraction)
    main(genFile, qtlFile, outFile) 

run_subset_plink4Matched_proc.sh

#!/bin/bash

#SBATCH --job-name=run_subset_plink4Matched_proc
#SBATCH --account=pi-yangili1
#SBATCH --time=36:00:00
#SBATCH --output=subset_plink4Matched_proc.out
#SBATCH --error=subset_plink4Matched_proc.err
#SBATCH --partition=broadwl
#SBATCH --mem=30G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env


for i  in {1..22};
do
python subset_plink4Matched_proc.py ${i} "Total"
done

for i  in {1..22};
do
python subset_plink4Matched_proc.py ${i} "Nuclear"
done

Cat and remove indels:
/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/FractionApaMatch_LD/

cat chr* > allChr.TotalMatch.LD.geno.ld
grep -v indel allChr.TotalMatch.LD.geno.ld > allChr.TotalMatch.LD.geno.ld_noIndel

cat chr* > allChr.NuclearMatch.LD.geno.ld
grep -v indel allChr.NuclearMatch.LD.geno.ld > allChr.NuclearMatch.LD.geno.ld_noIndel

make into bed files:

makeAlloverlapbed_Matched.py

#load files:  

QTL_total=open("/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Total_matched_snps_GEUFormat.txt", "r")
QTL_nuclear=open("/project2/gilad/briana/threeprimeseq/data/MatchedSnp/Nuclear_matched_snps_GEUFormat.txt", "r")
LD_total=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/TotalApaMatch_LD/allChr.TotalMatch.LD.geno.ld_noIndel", "r")
LD_nuclear=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/NuclearApaMatch_LD/allChr.NuclearMatch.LD.geno.ld_noIndel", "r")
outFile= open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/AllOverlapMatchSnps.bed", "w")

#function for qtl to bed format
def qtl2bed(fqtl, fraction, fout=outFile):
    for ln in fqtl:
        snp, chrom, pos = ln.split("_")
        start=int(pos)-1
        end= int(pos)
        fout.write("%s\t%d\t%d\tQTL_%s\n"%(chrom, start, end,fraction))

#function for ld to bed format 
def ld2bed(fLD, fraction, fout=outFile):
    for ln in fLD:
        snpID=ln.split()[5]
        snp, chrom, pos= snpID.split("_")
        start=int(pos)-1
        end=int(pos)
        fout.write("%s\t%d\t%d\tLD_%s\n"%(chrom, start, end,fraction))


#I will run each of these for both fractions to get all of the snps in the out file. 


qtl2bed(QTL_nuclear, "Nuclear")
qtl2bed(QTL_total, "Total")
ld2bed(LD_nuclear, "Nuclear")
ld2bed(LD_total, "Total")


outFile.close()

Run the gwas ovelap:

run_overlapSNPsGWASFixed_match.sh

#!/bin/bash

#SBATCH --job-name=run_overlapSNPsGWASFixed_match
#SBATCH --account=pi-yangili1
#SBATCH --time=5:00:00
#SBATCH --output=run_overlapSNPsGWASFixed_match.out
#SBATCH --error=run_overlapSNPsGWASFixed_match.err
#SBATCH --partition=broadwl
#SBATCH --mem=10G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env


python overlapSNPsGWAS_fixed.py  "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/AllOverlapMatchSnps.bed" "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_matched/GWASOverlap_AllOverlapMatchSnps.bed"

MatchOverlap=fread("../data/GWAS_overlap/GWASOverlap_AllOverlapMatchSnps.bed", header=F, col.names = c("chromSnp", "startSnp", "endSnp", "Set", "chromGwas", "startGWAS", "endGWAS", "rsID", "score", "strand"))

MatchOverlap_Total=MatchOverlap %>% filter(grepl("Total",Set))

MatchOverlap_Nuclear=MatchOverlap %>% filter(grepl("Nuclear",Set))

Compare:

MatchTot= MatchOverlap_Total %>% select(rsID) %>% unique() %>% nrow()
MatchNuc= MatchOverlap_Nuclear  %>% select(rsID) %>% unique() %>% nrow()
QTL_tot=QTL_overlap_Total %>% select(rsID) %>% unique() %>% nrow()
QTL_nuc= QTL_overlap_Nuclear %>% select(rsID) %>% unique() %>% nrow()

Number of uniq snps in GWAS catelog

GWAS_snp=53135726

totpval= phyper(QTL_tot, GWAS_snp, GWAS_snp,MatchTot+QTL_tot, lower.tail = F )

nucpval= phyper(QTL_nuc, GWAS_snp, GWAS_snp,MatchNuc+QTL_nuc ,lower.tail = F)

GWASdf=as.data.frame(cbind(Fraction=c("Total", "Total", "Nuclear", "Nuclear"),Type=c("QTL", "Match", "QTL", "Match"),Value=c(QTL_tot,MatchTot, QTL_nuc, MatchNuc)))
GWASdf$Value=as.numeric(as.character(GWASdf$Value))


anno_df=data.frame(Type="QTL",Value=6300,
                      Fraction = factor("Nuclear",levels = c("Nuclear","Total")))
anno_df2=data.frame(Type="QTL",Value=6000,
                      Fraction = factor("Nuclear",levels = c("Nuclear","Total")))

GwasOverlap=ggplot(GWASdf, aes(x=Type, by=Type, fill=Type, y=Value)) + geom_bar(position="dodge", stat="identity") + facet_grid(~Fraction) + labs(y="N overlap snps",title="apaQTLs overlap with GWAS catelog")+ geom_text(data = anno_df,aes(label="P < .0001")) + scale_fill_manual(values= c("Grey", "Blue")) + geom_text(data = anno_df2,aes(label="***"))
GwasOverlap

ggsave(GwasOverlap, file="../output/plots/apaQTLsoverlapGWASCatelog.png")

Saving 7 x 5 in image

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] data.table_1.12.0 forcats_0.4.0     stringr_1.4.0    
 [4] dplyr_0.8.0.1     purrr_0.3.1       readr_1.3.1      
 [7] tidyr_0.8.3       tibble_2.0.1      ggplot2_3.1.0    
[10] tidyverse_1.2.1   workflowr_1.2.0  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0       cellranger_1.1.0 plyr_1.8.4       pillar_1.3.1    
 [5] compiler_3.5.1   git2r_0.24.0     tools_3.5.1      digest_0.6.18   
 [9] lubridate_1.7.4  jsonlite_1.6     evaluate_0.13    nlme_3.1-137    
[13] gtable_0.2.0     lattice_0.20-38  pkgconfig_2.0.2  rlang_0.3.1     
[17] cli_1.0.1        rstudioapi_0.9.0 yaml_2.2.0       haven_2.1.0     
[21] xfun_0.5         withr_2.1.2      xml2_1.2.0       httr_1.4.0      
[25] knitr_1.21       hms_0.4.2        generics_0.0.2   fs_1.2.6        
[29] rprojroot_1.3-2  grid_3.5.1       tidyselect_0.2.5 glue_1.3.0      
[33] R6_2.4.0         readxl_1.3.0     rmarkdown_1.11   reshape2_1.4.3  
[37] modelr_0.1.4     magrittr_1.5     whisker_0.3-2    backports_1.1.3 
[41] scales_1.0.0     htmltools_0.3.6  rvest_0.3.2      assertthat_0.2.0
[45] colorspace_1.4-0 labeling_0.3     stringi_1.3.1    lazyeval_0.2.1  
[49] munsell_0.5.0    broom_0.5.1      crayon_1.3.4