Characterize apaQTLs

Chromatin regions for QTLs

This will start with the sig QTL bed files:

• /project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuclear.apaQTLs.sort.bed

• /project2/gilad/briana/threeprimeseq/data/ApaQTLs/Total.apaQTLs.sort.bed

Creare a python script using pybedtools:

ChromHMM.QTLs.py

import pybedtools 

sigNuc=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuclear.apaQTLs.sort.bed') 

sigTot=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Total.apaQTLs.sort.bed')

hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")

#map hmm to snps  
Tot_overlapHMM=sigTot.map(hmm, c=4)

Nuc_overlapHMM=sigNuc.map(hmm,c=4)

#save results  

Tot_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Tot_overlapHMM.bed")

Nuc_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuc_overlapHMM.bed")

Run in three-prime-env

chromHmm=read.table("../data/ChromHmmOverlap/chromHMM_regions.txt", col.names = c("number", "name"), stringsAsFactors = F)
TotalOverlapHMM=read.table("../data/ApaQTLs/Tot_overlapHMM.bed", col.names=c("chrom", "start", "end", "pid", "significance", "strand", "number"))
TotalOverlapHMM_names=TotalOverlapHMM %>% left_join(chromHmm, by="number")

NuclearOverlapHMM=read.table("../data/ApaQTLs/Nuc_overlapHMM.bed", col.names=c("chrom", "start", "end", "pid", "significance", "strand", "number"))
NuclearOverlapHMM_names=NuclearOverlapHMM %>% left_join(chromHmm, by="number")


ggplot(TotalOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Total APAQtls" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))

Past versions of unnamed-chunk-3-1.png

ggplot(NuclearOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Nuclear APAQtls" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))

Past versions of unnamed-chunk-3-2.png

Group them to put on plot together

NuclearOverlapHMM_names_byname= NuclearOverlapHMM_names %>% group_by(name) %>% summarise(Nuclear=n())
TotalOverlapHMM_names_byname= TotalOverlapHMM_names %>% group_by(name) %>% summarise(Total=n())


bothFracHMM= TotalOverlapHMM_names_byname %>% full_join(NuclearOverlapHMM_names_byname, by="name")
bothFracHMM$Nuclear= bothFracHMM$Nuclear %>% replace_na(0)

melt and plot

bothFracHMM_melt= melt(bothFracHMM,id.vars="name")
colnames(bothFracHMM_melt)=c("Region", "Fraction", "NQTLs")

apaQTLsChromregion=ggplot(bothFracHMM_melt, aes(x=Region, y=NQTLs, by=Fraction, fill=Fraction)) + geom_bar(stat="identity", position = "dodge")+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +scale_fill_manual(values=c("darkviolet","deepskyblue3"))  + labs(title="apaQTLs by chromatin region")
apaQTLsChromregion

Past versions of unnamed-chunk-5-1.png

ggsave(apaQTLsChromregion, file="../output/plots/apaQTLsbyChromHMM.png", width = 7, height = 5)

bothFracHMM_prop= bothFracHMM %>% mutate(Total_prop=Total/291) %>% mutate(Nuclear_prop=Nuclear/615) %>% select(-Total,-Nuclear) 

bothFracHMM_prop_melt= melt(bothFracHMM_prop,id.vars="name")
colnames(bothFracHMM_prop_melt)=c("Region", "Fraction", "PropQTLs")


propapaQTLsChromregion=ggplot(bothFracHMM_prop_melt, aes(x=Region, y=PropQTLs, by=Fraction, fill=Fraction)) + geom_bar(stat="identity", position = "dodge")+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +scale_fill_manual(values=c("darkviolet","deepskyblue3"))  + labs(y="Proportion of apaQTLs", title="Proportion of apaQTLs by chromatin region")
propapaQTLsChromregion

Past versions of unnamed-chunk-6-1.png

ggsave(propapaQTLsChromregion, file="../output/plots/proportionapaQTLsbyChromHMM.png", width = 7, height = 5)

I will need to get matched snps to look for enrichment in these values.

Overlap apaQTLs between fractions.

Goal: Find the nominal pvalue for the significant snp peak pair in oposite fraction. I can make a dictionary with the total and nuclear QTLs then run through the nominal files to get the pvalues:

Start with apa QTLs:

• /project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuclear.apaQTLs.sort.bed

• /project2/gilad/briana/threeprimeseq/data/ApaQTLs/Total.apaQTLs.sort.bed

NomResFromOppFrac.py

nucQTLs="/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Nuclear.apaQTLs.sort.bed"  
totQTLs="/project2/gilad/briana/threeprimeseq/data/ApaQTLs/Total.apaQTLs.sort.bed"  

nucNom="/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Nuclear.fixed.pheno_5perc.fc.gz.qqnorm_allNomRes.txt"
totNom="/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno_NoMP_sm_quant.Total.fixed.pheno_5perc.fc.gz.qqnorm_allNomRes.txt"

outnuc="/project2/gilad/briana/threeprimeseq/data/QTL_Overlap/TotQTLs_inNucFractionRes.txt"
outtot="/project2/gilad/briana/threeprimeseq/data/QTL_Overlap/NucQTLs_inTotFractionRes.txt"


def oppFract(inRes, inQTL, out):
    fout=open(out, "w")
    qtl_dic={}
    #SNP is key, peak is value
    for ln in open(inQTL,"r"):
        snp=ln.split()[2]
        chrom=ln.split()[0]
        peak=ln.split()[3].split(":")[0]
        qtl=str(chrom) + ":" + str(snp)
        if qtl not in qtl_dic.keys():
            qtl_dic[qtl]=[peak]
        else: 
            qtl_dic[qtl].append(peak)
    #print(qtl_dic)
    for ln in open(inRes, "r"):
        pval=ln.split()[3]
        snp=ln.split()[1]
        peak=ln.split()[0].split(":")[3].split("_")[-1]
        if snp in qtl_dic.keys():
            if peak in qtl_dic[snp]:
                fout.write("%s\t%s\t%s\n"%(snp, peak, pval))   
    fout.close()
    

oppFract(nucNom, totQTLs,outnuc)  
oppFract(totNom, nucQTLs, outtot)  
                

Run in bash: run_NomResFromOppFrac.sh

#!/bin/bash

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


module load Anaconda3
source activate three-prime-env


python NomResFromOppFrac.py

names=c("SNP", "peak", "pval")
NucQTLinTot=read.table("../data/QTL_overlap/NucQTLs_inTotFractionRes.txt", stringsAsFactors = F, col.names = names)
TotQTLinNuc=read.table("../data/QTL_overlap/TotQTLs_inNucFractionRes.txt", stringsAsFactors = F, col.names = names)

Get pi values:

qvalTot= pi0est(NucQTLinTot$pval, pi0.method = "bootstrap")
1-qvalTot$pi0

[1] 0.8424242

qvalNuc= pi0est(TotQTLinNuc$pval, pi0.method = "bootstrap")
1-qvalNuc$pi0

[1] 0.9197861

plots:

par(mfrow=c(1,2))
hist(NucQTLinTot$pval,xlab="Total apaQTL pvalue", main="Nuclear apaQTLs \nin Total Fraction")
text(.6,200, paste("pi_1=", round((1-qvalTot$pi0), digit=3), sep=" "))
hist(TotQTLinNuc$pval,xlab="Nuclear apaQTL pvalue", main="Total apaQTLs \nin Nuclear Fraction")
text(.6,125, paste("pi_1=", round((1-qvalNuc$pi0), digit=3), sep=" "))

Past versions of unnamed-chunk-12-1.png

png("../output/plots/apaFractionOverlapPi1.png", width=1000, height = 500)
par(mfrow=c(1,2))
hist(NucQTLinTot$pval,xlab="Total apaQTL pvalue", main="Nuclear apaQTLs \nin Total Fraction")
text(.6,200, paste("pi_1=", round((1-qvalTot$pi0), digit=3), sep=" "))
hist(TotQTLinNuc$pval,xlab="Nuclear apaQTL pvalue", main="Total apaQTLs \nin Nuclear Fraction")
text(.6,125, paste("pi_1=", round((1-qvalNuc$pi0), digit=3), sep=" "))
dev.off()

quartz_off_screen 
                2 

This provides evidence for high degree of QTL sharing with increased sharing total to nuclear. This demonstrates to me that there are nuclear QTLs that do not persist in the total fraction. I will want to learn more about these.

Overlap with GWAS catelog

I did this analysis with the QTLs in the preprocessed 39 individual analysis. I will follow a similar pipeline here. I will find all of the snps in LD with the QTLs then test for these in the GWAS catelog. The pipeline I used to get the LD for all of the snp is shown here. The plink files are in /project2/gilad/briana/threeprimeseq/data/GWAS_overlap/. There are both map and ped files.

I can now adapt the subset_plink4QTLs.py file to take the current QTLs list. The file just has the QTLs with the chromosome and position. I can make this and put it in:

/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed

The 50mb QTLs are in /project2/gilad/briana/threeprimeseq/data/ApaQTLs.

• NuclearapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt
  
• TotalapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt

The QTL snps are in the 6th column.

cut -f6 -d" " /project2/gilad/briana/threeprimeseq/data/ApaQTLs/NuclearapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt | uniq > /project2/gilad/briana/threeprimeseq/data/ApaQTLs/NuclearQTLs_uniq_50mb.txt
cut -f6 -d" "  /project2/gilad/briana/threeprimeseq/data/ApaQTLs/TotalapaQTLs.GeneLocAnno.noMP.5perc.10FDR.txt | uniq > /project2/gilad/briana/threeprimeseq/data/ApaQTLs/TotalQTLs_uniq_50mb.txt

I can convert these the the way they are in GEU snp files tony made (snp_num_pos)

QTLs2GeuSnps_proc.py

tot_in=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/TotalQTLs_uniq_50mb.txt", "r")  
nuc_in=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/NuclearQTLs_uniq_50mb.txt", "r")

tot_out=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/TotalQTLs_uniq_50mb_GEU.txt", "w") 
nuc_out=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/NuclearQTLs_uniq_50mb_GEU.txt", "w") 


def fix_file(fin, fout):
  for ln in fin:
    chrom, pos = ln.split(":")
    fout.write("snp_%s_%s"%(chrom,pos))
  fout.close()
  

fix_file(tot_in, tot_out)
fix_file(nuc_in, nuc_out)

subset_plink4QTLs_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_processed/%sApaQTL_LD/chr%s.%sQTL.LD.geno.ld"%(fraction,chrom,fraction)
    qtlFile= "/project2/gilad/briana/threeprimeseq/data/ApaQTLs/%sQTLs_uniq_50mb_GEU.txt"%(fraction)
    main(genFile, qtlFile, outFile) 

run_subset_plink4QTLs_proc.sh

#!/bin/bash

#SBATCH --job-name= run_subset_plink4QTLs_proc
#SBATCH --account=pi-yangili1
#SBATCH --time=36:00:00
#SBATCH --output=run_subset_plink4QTLs_proc.out
#SBATCH --error=run_subset_plink4QTLs_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_plink4QTLs_proc.py ${i} "Total"
done

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

This added 2446 total snps and 6258 nuclear snps.

Cat and remove indels:

cat chr* > allChr.TotalQTL.LD.gene.ld
grep -v indel allChr.TotalQTL.LD.gene.ld > allChr.TotalQTL.LD.gene.ld_noIndel

cat chr* > allChr.NuclearQTL.LD.gene.ld
grep -v indel allChr.NuclearQTL.LD.gene.ld > allChr.NuclearQTL.LD.gene.ld_noIndel

Make these bed files:

makeAlloverlapbed_proc.py

#load files:  

QTL_total=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/TotalQTLs_uniq_50mb_GEU.txt", "r")
QTL_nuclear=open("/project2/gilad/briana/threeprimeseq/data/ApaQTLs/NuclearQTLs_uniq_50mb_GEU.txt", "r")
LD_total=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/TotalApaQTL_LD/allChr.TotalQTL.LD.gene.ld_noIndel", "r")
LD_nuclear=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/NuclearApaQTL_LD/allChr.NuclearQTL.LD.gene.ld_noIndel", "r")
outFile= open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllOverlapSnps.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()

Sort this:

sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllOverlapSnps.bed > /project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllOverlapSnps.sort.bed

Overlap with GWAS

I can use the overlapSNPsGWAS.py file I created in the previous rendition of this analysis but run it with these files.

run_overlapSNPsGWAS_proc.sh

#!/bin/bash

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


module load Anaconda3
source activate three-prime-env

python overlapSNPsGWAS.py  "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllOverlapSnps.sort.bed" "/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllSnps_GWASoverlapped.txt"

Total QTLs overlap: rs3117582 6:31620520

Total LD overlap:

• rs2282301 1:155868625
  
• rs3596 12:95696420
  
• rs2277862 20:34152782
  
• rs2517713 6:29918099

Nuclear QTL overlap:

rs7206971 17:45425115

Nucelar LD overlapL

• rs10889353 1:63118196
  
• rs2282301 1:155868625
• rs10859871 12:95711876
  
• rs10133111 14:103377321
• rs17382723 2:242053546
  
• rs2277862 20:34152782
  
• rs2298428 22:21982892
  
• rs13160562 5:96111371
  
• rs29784 5:172595308
  
• rs2517713 6:29918099
  
• rs3077 6:33033022

Are these eQTLs?

Pull these snps from the eQTLs nominal pvalues:

Process the AllSnps_GWASoverlapped.txt file to have the RS id, the snp, and where it was found.

formatGWASOverlap.py

inFile="/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllSnps_GWASoverlapped.txt"

outFile=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllSnps_GWASoverlapped_format.txt","w")

for i in open(inFile, "r"):
   chrom=i.split()[0]
   loc=i.split()[2]
   snp=chrom + ":" + loc
   found=i.split()[3]
   rs=i.split()[7].split(":")[0]
   outFile.write("%s\t%s\t%s\n"%(snp, rs, found))
outFile.close()

This is 14 uniq snps.

/project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.out

I need to go throguh this file and check for the snps. I can do this by making a dictionary of the uniq snps in the overlap then testing the lines ib the nominal file. I want to keep the gene, snp, and pvalue for these associations
getexpPvalOverlap.py

overlapSnps="/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/AllSnps_GWASoverlapped_format.txt"
outF=open("/project2/gilad/briana/threeprimeseq/data/GWAS_overlap_processed/eQTL_pval_GWASOverlap.txt","w")  
nomRes="/project2/gilad/briana/threeprimeseq/data/molecular_QTLs/nom/fastqtl_qqnorm_RNAseq_phase2.fixed.nominal.out"  

snps_dic={}
for ln in open(overlapSnps, "r"):
    snp=ln.split()[0]
    if snp not in snps_dic.keys():
        snps_dic[snp]=""  

for ln in open(nomRes,"r"):
    snp=ln.split()[1]
    if snp in snps_dic.keys():
        gene=ln.split()[0].split(".")[0]
        pval=ln.split()[3]
        outF.write("%s\t%s\t%s\n"%(snp,gene,pval))
outF.close()

I can pull this in and get the lowest pval for each one

eQTLpvalOverlapGWAS=read.table("../data/ApaQTLs/eQTL_pval_GWASOverlap.txt", stringsAsFactors = F, col.names = c("snp", "gene", "pval"))

eQTLpvalOverlapGWAS_min= eQTLpvalOverlapGWAS %>% group_by(snp) %>% summarise(pvalM=min(pval)) %>% mutate(noSig=ifelse(pvalM<.05, "yes", "no"))

FIlter by the non sig ones to look at as exmaples:

eQTLpvalOverlapGWAS_NotSig= eQTLpvalOverlapGWAS_min %>% filter(noSig=="no")

eQTLpvalOverlapGWAS_NotSig

# A tibble: 5 x 3
  snp          pvalM noSig
  <chr>        <dbl> <chr>
1 1:63118196  0.400  no   
2 12:95696420 0.0813 no   
3 12:95711876 0.0813 no   
4 2:242053546 0.0878 no   
5 5:172595308 0.0530 no   

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] bindrcpp_0.2.2  reshape2_1.4.3  qvalue_2.12.0   forcats_0.3.0  
 [5] stringr_1.4.0   dplyr_0.7.6     purrr_0.2.5     readr_1.1.1    
 [9] tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0   tidyverse_1.2.1
[13] workflowr_1.2.0

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