ChromHMM analysis

Compare enrichment between fractions

I want to make a plot with the enrichment by fraction. I am first going to get an enrichemnt score for each bin naively by looking at the QTL/random in each category.

#perChrommHMM_nuc$Random= as.integer(perChrommHMM_nuc$Random)
#perChrommHMM_nuc_enr=perChrommHMM_nuc %>%  mutate(Nuclear=Nuclear_QTLs-Random)

#perChrommHMM_tot_enr=read.table("../data/ChromHmmOverlap/perChrommHMM_Total_enr.txt",stringsAsFactors = F,header = T)

#allenrich=perChrommHMM_tot_enr %>% inner_join(perChrommHMM_nuc_enr, by="name") %>% select(name, Total, Nuclear)

#allenrich_melt=melt(allenrich, id.vars="name")

plot it

#chromenrichBoth=ggplot(allenrich_melt, aes(x=name, by=variable, y=value, fill=variable)) + geom_bar(stat="identity", position = "dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="QTL-Random for each bin by fraction", y="Num QTL SNPs - Num Random SNPs") + scale_fill_manual(values=c("darkviolet", "deepskyblue3"))


#ggsave("../output/plots/ChromHmmEnrich_BothFrac.png", chromenrichBoth)

Remove significant genes

NucQTL_genes=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_transcript_permResBH.txt", stringsAsFactors = F, header=T)  %>% mutate(sig=ifelse(-log10(bh)>=1, 1,0 )) %>%  separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% filter(sig==1) %>% select(gene) %>% distinct(gene)
#715 genes  
#write this out as NucAPAGenes
write.table(NucQTL_genes, "../data/perm_QTL_trans/NucApaGenes.txt", row.names = F, col.names = F, quote=F)

TotQTL_genes=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt", stringsAsFactors = F, header=T)  %>% mutate(sig=ifelse(-log10(bh)>=1, 1,0 )) %>%  separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% filter(sig==1) %>% select(gene) %>% distinct(gene)
#106 genes
#write out as TotAPAGenes

write.table(TotQTL_genes, "../data/perm_QTL_trans/TotApaGenes.txt", row.names = F, col.names = F, quote=F)

I need to find a way to get rid of these from the files I cam pulling from.

/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_NomRes.txt

/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_NomRes.txt

I can create an python script to do this. I will need to seperate the first column and and only write the line out if the gene is in the apaGenes files I just created.

filterSigGenes.py

#python 

#genes with sig ApaQTL
TotGenes=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/sig_genes/TotApaGenes.txt", "r")
NucGenes=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/sig_genes/NucApaGenes.txt", "r")

#nom res (with all snps tested)  
NucRes=open("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_NomRes.txt", "r")
TotRes=open("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_NomRes.txt", "r")

#output files:
Nuc_nonSig=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt", "w")
Tot_nonSig=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt", "w")

#convert genes to list
def file_to_list(file):
    gene_list=[]
    for ln in file:
      gene=ln.strip()
      gene_list.append(gene)
    return(gene_list)
    
Tot_gene_list=file_to_list(TotGenes)
Nuc_gene_list=file_to_list(NucGenes)  

#function that will take in the input, the list, and the output. I want to be able to run this function for total and nuclear  


def filter(fin,fout, sigGenes):
  for ln in fin:
    gene=ln.split()[0].split(":")[3].split("_")[0]
    if gene not in sigGenes:
      fout.write(ln)
  fout.close()


filter(NucRes,Nuc_nonSig,Nuc_gene_list)
filter(TotRes, Tot_nonSig, Tot_gene_list)

Call this in a bash script:
run_filterSigGenes.sh

#!/bin/bash

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


module load Anaconda3
source activate three-prime-env


python filterSigGenes.py

nuc_random880_chromHmm_noSig_set1.sh

#!/bin/bash

#SBATCH --job-name=nuc_random880_chromHmm_noSig_set1
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=nuc_random880_chromHmm_noSig_set1.out
#SBATCH --error=nuc_random880_chromHmm_noSig_set1.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env
#make random 
for i in {1..250};
do
shuf -n 880 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/randomRes_Nuclear_880_noSig_${i}.txt
done

nuc_random880_chromHmm_noSig_set2.sh

#!/bin/bash

#SBATCH --job-name=nuc_random880_chromHmm_noSig_set2
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=nuc_random880_chromHmm_noSig_set2.out
#SBATCH --error=nuc_random880_chromHmm_noSig_set2.err
#SBATCH --partition=bigmem2
#SBATCH --mem=200G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env
#make random 
for i in {251..500};
do
shuf -n 880 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/randomRes_Nuclear_880_noSig_${i}.txt
done

nuc_random880_chromHmm_noSig_set3.sh

#!/bin/bash

#SBATCH --job-name=nuc_random880_chromHmm_noSig_set3
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=nuc_random880_chromHmm_noSig_set3.out
#SBATCH --error=nuc_random880_chromHmm_noSig_set3.err
#SBATCH --partition=bigmem2
#SBATCH --mem=200G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env
#make random 
for i in {501..750};
do
shuf -n 880 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/randomRes_Nuclear_880_noSig_${i}.txt
done

nuc_random880_chromHmm_noSig_set4.sh

#!/bin/bash

#SBATCH --job-name=nuc_random880_chromHmm_noSig_set4
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=nuc_random880_chromHmm_noSig_set4.out
#SBATCH --error=nuc_random880_chromHmm_noSig_set4.err
#SBATCH --partition=bigmem2
#SBATCH --mem=200G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env
#make random 
for i in {751..1000};
do
shuf -n 880 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/randomRes_Nuclear_880_noSig_${i}.txt
done

Same for total:

total_random118_chromHmm_noSig_set1.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


#test with 2 permutations then make it 1000  
#choose random res
for i in {1..250};
do
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/randomRes_Total_118_noSig_${i}.txt
done

total_random118_chromHmm_noSig_set2.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


#test with 2 permutations then make it 1000  
#choose random res
for i in {251..500};
do
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/randomRes_Total_118_noSig_${i}.txt
done

total_random118_chromHmm_noSig_set3.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


#test with 2 permutations then make it 1000  
#choose random res
for i in {501..750};
do
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/randomRes_Total_118_noSig_${i}.txt
done

total_random118_chromHmm_noSig_set4.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


#test with 2 permutations then make it 1000  
#choose random res
for i in {751..1000};
do
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/randomRes_Total_118_noSig_${i}.txt
done

This may not be enough. I may need to change this so it only has uniq snp. I could be sampling the same snps over an over.

Make these files into snp bed files:

randomRes2SNPbed_noSig.py

def main(inFile, outFile):
    fout=open(outFile, "w")
    fin=open(inFile, "r")
    for ln in fin:
          pid, sid, dist, pval, slope = ln.split()
          chrom, pos= sid.split(":")
          name=sid
          start= int(pos)-1
          end=int(pos)
          strand=pid.split(":")[3].split("_")[1]
          pval=float(pval)
          fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, pval, strand))
    fout.close()

if __name__ == "__main__":
    import sys
    fraction=sys.argv[1]
    nsamp=sys.argv[2]
    nsamp=int(nsamp)
    iter=sys.argv[3]
    inFile = "/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_noSig/randomRes_%s_%d_noSig_%s.txt"%(fraction,fraction, nsamp, iter)
    
    outFile= "/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_noSig/snp_bed_noSig/randomRes_%s_%d_noSig_%s.bed"%(fraction,fraction, nsamp, iter)
    main(inFile,outFile)

randomLines2Snp_noSig.sh

#!/bin/bash

#SBATCH --job-name=randomLines2Snp_noSig
#SBATCH --account=pi-gilad
#SBATCH --time=36:00:00
#SBATCH --output=randomLines2Snp_noSig.out
#SBATCH --error=randomLines2Snp_noSig.err
#SBATCH --partition=gilad
#SBATCH --mem=50G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env


#make random 
for i in {1..1000};
do
python randomRes2SNPbed_noSig.py Nuclear 880 ${i} 
done 

#make random 
for i in {1..1000};
do
python randomRes2SNPbed_noSig.py Total 118 ${i}
done 

sortRandomSnps_noSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


for i in $(ls /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/snp_bed_noSig/);
do
sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/snp_bed_noSig/$i > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/snp_bed_sort_noSig/$i.sort.bed
done

for i in $(ls /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/snp_bed_noSig/);
do
sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/snp_bed_noSig/$i > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/snp_bed_sort_noSig/$i.sort.bed
done

overlap_chromHMM_sepfiles_noSig.py

def main(inFile, outFile):
  rand=pybedtools.BedTool(inFile) 
  hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")
  #map hmm to snps
  Rand_overlapHMM=rand.map(hmm, c=4)
  #save results
  Rand_overlapHMM.saveas(outFile)


if __name__ == "__main__":
    import sys
    import pandas as pd
    import pybedtools
    fraction=sys.argv[1]
    nsamp=sys.argv[2]
    niter=sys.argv[3]
    #which itteration we are on 
    inFile ="/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_noSig/snp_bed_sort_noSig/randomRes_%s_%s_noSig_%s.bed.sort.bed"%(fraction,fraction, nsamp, niter)

    outFile= "/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_noSig/chromHMM_overlap/randomres_overlapChromHMM_%s_%s_noSig_%s.txt"%(fraction,fraction,nsamp, niter)
    main(inFile,outFile)

overlap_chromHMM_sepfiles_noSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

for i in {1..1000};
do
python overlap_chromHMM_sepfiles_noSig.py  Nuclear 880 $i
done

for i in {1..1000};
do
python overlap_chromHMM_sepfiles_noSig.py  Total 118 $i
done

groupRandomChromHMM_noSig.sh

#!/bin/bash

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


module load Anaconda3
source activate three-prime-env



for i in {1..1000};
do
Rscript groupRandomByChromHMM.R -f /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/chromHMM_overlap/randomres_overlapChromHMM_Nuclear_880_noSig_${i}.txt -o /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_noSig/chromHMM_overlap_group/randomres_overlapChromHMM_Nuclear_880_noSig_${i}_grouped.txt
done


for i in {1..1000};
do
Rscript groupRandomByChromHMM.R -f /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/chromHMM_overlap/randomres_overlapChromHMM_Total_118_noSig_${i}.txt  -o /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_noSig/chromHMM_overlap_group/randomres_overlapChromHMM_Total_118_noSig_${i}_grouped.txt
done

cut -d$' ' -f 1,2 randomres_overlapChromHMM_Nuclear_880_noSig_549_grouped.txt > Nuc_chromOverlap_noSig.txt

for i in {1..1000};
do
paste -d" " Nuc_chromOverlap_noSig.txt <(cut -d" " -f 3 randomres_overlapChromHMM_Nuclear_880_noSig_${i}_grouped.txt) > tmp
mv tmp Nuc_chromOverlap_noSig.txt
done


cut -d$' ' -f 1,2 randomres_overlapChromHMM_Total_118_noSig_53_grouped.txt > Tot_chromOverlap_noSig.txt

for i in {1..1000};
do
paste -d" " Tot_chromOverlap_noSig.txt <(cut -d" " -f 3 randomres_overlapChromHMM_Total_118_noSig_${i}_grouped.txt) > tmp
mv tmp Tot_chromOverlap_noSig.txt
done

These files are not on my computer so I can work with them.

permutationResTotal_noSig=read.table("../data/ChromHmmOverlap/Tot_chromOverlap_noSig.txt", header=T, stringsAsFactors = F)
permutationResTotal_noSig[is.na(permutationResTotal_noSig)] <- 0
permutationResTotal_noSig= as_data_frame(permutationResTotal_noSig)
permutationResTotal_noSig_noName=permutationResTotal_noSig[,3:ncol(permutationResTotal_noSig)]
totRand_mean_noSig=rowMeans(permutationResTotal_noSig_noName)/1000

permutationResNuclear_noSig=read.table("../data/ChromHmmOverlap/Nuc_chromOverlap_noSig.txt",header = T,stringsAsFactors = F)
permutationResNuclear_noSig[is.na(permutationResNuclear_noSig)] <- 0
permutationResNuclear_noSig_noName=permutationResNuclear_noSig[,3:ncol(permutationResNuclear_noSig)]
nucRand_mean_noSig=rowMeans(permutationResNuclear_noSig_noName)/1000

permutationResTotal_noSig_melt= melt(permutationResTotal_noSig, id.vars=c("number", "name"))

ggplot(permutationResTotal_noSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + labs(x="N random Snps in category", title="Random permutation Total")

For nuclear:

permutationResNuclear_noSig_melt= melt(permutationResNuclear_noSig, id.vars=c("number", "name"))

ggplot(permutationResNuclear_noSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + labs(x="N random Snps in category", title="Random permutation Nuclear")

ggplot(permutationResNuclear_noSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + scale_y_log10()+ labs(x="random Snps in category", title="Random permutation Nuclear")

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 630 rows containing missing values (geom_bar).

ggplot(permutationResTotal_noSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + scale_y_log10()+ labs(x="random Snps in category", title="Random permutation Total")

Warning: Transformation introduced infinite values in continuous y-axis

Warning: Removed 445 rows containing missing values (geom_bar).

Permute only from unique snps tested

This didnt help. I need to choose from the unique snps rather than counting if they were tested multiple times. I will look for the snps tested and get rid of those called as QTLs.

Steps:

• get total and nuclear snps tested

I need to do this from the nominal results in /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/. The snp ID is in the second column

run_testedSnps_noSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


cut -f2 -d" " /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/NucTestedSnps_nonSigGenes.txt | sort | uniq > /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/uniq_tested_Nuclear.txt


cut -f2 -d" " /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/nomRes_nonsig/TotTestedSnps_nonSigGenes.txt | sort | uniq > /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/uniq_tested_Total.txt

• remove QTL snps

The significant QTL snps are:

/project2/gilad/briana/threeprimeseq/data/GWAS_overlap/ApaQTLsigSNPpos_Nuclear.txt /project2/gilad/briana/threeprimeseq/data/GWAS_overlap/ApaQTLsigSNPpos_Total.txt

I can use python to remove the snps from the full lists:

I can write it as a bed file to make the step easier.

testedSnps_noSig.py

total_qtl=open('/project2/gilad/briana/threeprimeseq/data/GWAS_overlap/ApaQTLsigSNPpos_Total.txt', "r")
nuclear_qtl=open('/project2/gilad/briana/threeprimeseq/data/GWAS_overlap/ApaQTLsigSNPpos_Total.txt', "r")

def file_to_list(file):
    snp_list=[]
    for ln in file:
        snp=ln.strip()
        snp_list.append(snp)
    return(snp_list)
    
total_qtl_list= file_to_list(total_qtl)
nuclear_qtl_list= file_to_list(nuclear_qtl)

total_out=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/uniqNonSigSnps/TotalUniqTestedSnp.bed","w")
for ln in open("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/uniq_tested_Total.txt"):
    snp=ln.strip()
    if snp not in total_qtl_list:
        chrom, pos =snp.split(":")
        start=int(pos)-1
        end=int(pos)
        total_out.write("%s\t%d\t%d\t%s\n"%(chrom, start, end, snp))
total_out.close()  
    
nuclear_out=open("/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/uniqNonSigSnps/NuclearUniqTestedSnp.bed","w")

for ln in open("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/uniq_tested_Nuclear.txt"):
    snp=ln.strip()
    if snp not in total_qtl_list:
        chrom, pos =snp.split(":")
        start=int(pos)-1
        end=int(pos)
        nuclear_out.write("%s\t%d\t%d\t%s\n"%(chrom, start, end, snp))
nuclear_out.close()  
    

I will probablly need to run this in a bash script.

run_testedSnps2bed.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

python testedSnps_noSig.py  

• select correct number of snps 1000 times (permutation) This file is not as big. I may be able to run all of the permutations at once:
  random_UniqNoSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


for i in {1..1000};
do
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/uniqNonSigSnps/TotalUniqTestedSnp.bed  > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/snp_bed/randomRes_Total_118_UniqnoSig_${i}.bed
done


for i in {1..1000};
do
shuf -n 880 /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/uniqNonSigSnps/NuclearUniqTestedSnp.bed > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/snp_bed/randomRes_Nuclear_880_UniqnoSig_${i}.bed
done

• sort bed files

sortRandomUniqSnps_noSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


for i in $(ls /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/snp_bed);
do
sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/snp_bed/$i > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/snp_bed_sort_noSig/$i.sort.bed
done

for i in $(ls /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/snp_bed);
do
sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/snp_bed/$i > /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/snp_bed_sort_noSig/$i.sort.bed
done

• overlap each with chrom HMM

overlap_chromHMM_sepfiles_UniqnoSig.py

def main(inFile, outFile):
  rand=pybedtools.BedTool(inFile) 
  hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")
  #map hmm to snps
  Rand_overlapHMM=rand.map(hmm, c=4)
  #save results
  Rand_overlapHMM.saveas(outFile)


if __name__ == "__main__":
    import sys
    import pandas as pd
    import pybedtools
    fraction=sys.argv[1]
    nsamp=sys.argv[2]
    niter=sys.argv[3]
    #which itteration we are on 
    inFile ="/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_Uniq_noSig/snp_bed_sort_noSig/randomRes_%s_%s_UniqnoSig_%s.bed.sort.bed"%(fraction,fraction, nsamp, niter)

    outFile= "/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/%s_Uniq_noSig/chromHMM_overlap/randomres_overlapChromHMM_%s_%s_UniqnoSig_%s.txt"%(fraction,fraction,nsamp, niter)
    main(inFile,outFile)

overlap_chromHMM_sepfiles_UniqnoSig.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

for i in {1..1000};
do
python overlap_chromHMM_sepfiles_UniqnoSig.py  Nuclear 880 $i
done

for i in {1..1000};
do
python overlap_chromHMM_sepfiles_UniqnoSig.py  Total 118 $i
done

• group the chrom HMM calls

I need to modify groupRandomByChromHMM_uniq.R

#!/bin/rscripts

# usage: groupRandomByChromHMM_uniq.R -f infile -o outfile 

#this file will take any of the itterations and output a file with chrom hmm number, name, numberof snps

library(optparse)
library(dplyr)
library(tidyr)
library(ggplot2)
library(readr)

option_list = list(
  make_option(c("-f", "--file"), action="store", default=NA, type='character',
              help="input coverage file"),
  make_option(c("-o", "--output"), action="store", default=NA, type='character',
              help="output file")
)

opt_parser <- OptionParser(option_list=option_list)
opt <- parse_args(opt_parser)


#interrupt execution if no file is  supplied
if (is.null(opt$file)){
  print_help(opt_parser)
  stop("Need input file", call.=FALSE)
}
if (is.null(opt$output)){
  print_help(opt_parser)
  stop("Need output file", call.=FALSE)
}

randomSNPS=read.table(opt$file, col.names=c("chrom", "start", "end", "sid", "number"),stringsAsFactors = F, na.strings = "NA")
hmm_names=read.table("/project2/gilad/briana/genome_anotation_data/chromHMM_regions.txt", col.names = c("number", "name"),stringsAsFactors=F)
randomSNPS$number=as.character(randomSNPS$number)
hmm_names$number=as.character(hmm_names$number)
randomSNPS_names= randomSNPS  %>% left_join(hmm_names, by="number")
#split the name of the file to get the iteration number
fileSplit=strsplit(opt$file, "/")[[1]][10]
iter.txt=strsplit(fileSplit, "_")[[1]][5]
iter=substr(iter.txt, 1, nchar(iter.txt)-4) 

randomSNPS_names_grouped=randomSNPS_names %>%  group_by(number) %>% summarise(!!iter:=n()) %>%  replace_na(list(name="No_annotation")) %>%  dplyr::select(number, !!iter) 
hmm_names$number=as.character(hmm_names$number)
final=hmm_names %>% left_join(randomSNPS_names_grouped,by="number")

write.table(final,opt$output,quote=FALSE, col.names = T, row.names = F)

groupRandomChromHMM_UniqnoSig.sh

#!/bin/bash

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


module load Anaconda3
source activate three-prime-env



for i in {1..1000};
do
Rscript groupRandomByChromHMM_uniq.R -f /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/chromHMM_overlap/randomres_overlapChromHMM_Nuclear_880_UniqnoSig_${i}.txt -o /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/chromHMM_overlap_group/randomres_overlapChromHMM_Nuclear_880_UniqnoSig_${i}_grouped.txt
done


for i in {1..1000};
do
Rscript groupRandomByChromHMM_uniq.R -f /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/chromHMM_overlap/randomres_overlapChromHMM_Total_118_UniqnoSig_${i}.txt  -o /project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/chromHMM_overlap_group/randomres_overlapChromHMM_Total_118_UniqnoSig_${i}_grouped.txt
done

• make the final matrix

#/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Nuclear_Uniq_noSig/chromHMM_overlap_group/
cut -d$' ' -f 1,2 randomres_overlapChromHMM_Nuclear_880_UniqnoSig_549_grouped.txt > Nuc_chromOverlap_UniqnoSig.txt

for i in {1..1000};
do
paste -d" " Nuc_chromOverlap_UniqnoSig.txt <(cut -d" " -f 3 randomres_overlapChromHMM_Nuclear_880_UniqnoSig_${i}_grouped.txt) > tmp
mv tmp Nuc_chromOverlap_UniqnoSig.txt
done

#/project2/gilad/briana/threeprimeseq/data/random_QTLsnps/Total_Uniq_noSig/chromHMM_overlap_group/
cut -d$' ' -f 1,2 randomres_overlapChromHMM_Total_118_UniqnoSig_53_grouped.txt > Tot_chromOverlap_UniqnoSig.txt

for i in {1..1000};
do
paste -d" " Tot_chromOverlap_UniqnoSig.txt <(cut -d" " -f 3 randomres_overlapChromHMM_Total_118_UniqnoSig_${i}_grouped.txt) > tmp
mv tmp Tot_chromOverlap_UniqnoSig.txt
done

permutationResTotal_UniqnoSig=read.table("../data/ChromHmmOverlap/Tot_chromOverlap_UniqnoSig.txt", header=T, stringsAsFactors = F)
permutationResTotal_UniqnoSig[is.na(permutationResTotal_UniqnoSig)] <- 0
permutationResTotal_UniqnoSig= as_data_frame(permutationResTotal_UniqnoSig)
permutationResTotal_UniqnoSig_noName=permutationResTotal_UniqnoSig[,3:ncol(permutationResTotal_UniqnoSig)]
totRand_mean_UniqnoSig=rowMeans(permutationResTotal_UniqnoSig_noName)/1000

permutationResNuclear_UniqnoSig=read.table("../data/ChromHmmOverlap/Nuc_chromOverlap_UniqnoSig.txt",header = T,stringsAsFactors = F)
permutationResNuclear_UniqnoSig[is.na(permutationResNuclear_UniqnoSig)] <- 0
permutationResNuclear_UniqnoSig_noName=permutationResNuclear_UniqnoSig[,3:ncol(permutationResNuclear_UniqnoSig)]
nucRand_mean_UniqnoSig=rowMeans(permutationResNuclear_UniqnoSig_noName)/1000

permutationResTotal_UniqnoSig_melt= melt(permutationResTotal_UniqnoSig, id.vars=c("number", "name"))

ggplot(permutationResTotal_UniqnoSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + labs(x="N random Snps in category", title="Random permutation Total")

permutationResNuclear_UniqnoSig_melt= melt(permutationResNuclear_UniqnoSig, id.vars=c("number", "name"))

ggplot(permutationResNuclear_UniqnoSig_melt, aes(x=value,fill=name)) + geom_histogram(bins=50) + facet_wrap(~number) + labs(x="N random Snps in category", title="Random permutation Total")

This is the model I will move forward with.

Plot the enrichment:

Look at enrichment by using the average

TotalPermMean_uniq=permutationResTotal_UniqnoSig_melt %>% group_by(number) %>% summarise(TotRandPerm=mean(value), TotRandPermSD=sd(value))
TotalPermMean_uniq$number=as.character(TotalPermMean_uniq$number)
NuclearPermMean_uniq=permutationResNuclear_UniqnoSig_melt %>% group_by(number) %>% summarise(NucRandPerm=mean(value),NucRandPermSD=sd(value))
NuclearPermMean_uniq$number=as.character(NuclearPermMean_uniq$number)

Melt SNP values by name and number to get data in same format. I already did this above.

TotalOverlapHMM_names_melt=melt(TotalOverlapHMM_names, id.vars=c("number", "name"))%>% filter(variable=="sid") %>% group_by(number) %>% summarise(TotalQTL=n())
TotalOverlapHMM_names_melt$number=as.character(TotalOverlapHMM_names_melt$number)
NuclearOverlapHMM_names_melt=melt(NuclearOverlapHMM_names, id.vars=c("number", "name")) %>% filter(variable=="sid") %>% group_by(number) %>% summarise(NucQTL=n())
NuclearOverlapHMM_names_melt$number=as.character(NuclearOverlapHMM_names_melt$number)

chromHmm$number=as.character(chromHmm$number)
TotalOverlapHMM_Uniq_enrichment= TotalOverlapHMM_names_melt %>% full_join(TotalPermMean_uniq, by="number") %>%  replace_na(list(TotalQTL=.00001)) %>% full_join(chromHmm, by="number")

TotalOverlapHMM_Uniq_enrichment$TotalQTL=as.double(TotalOverlapHMM_Uniq_enrichment$TotalQTL)
TotalOverlapHMM_Uniq_enrichment = TotalOverlapHMM_Uniq_enrichment %>% mutate(TotEnrich=(TotalQTL-TotRandPerm)/TotRandPermSD)

NuclearOverlapHMM_Uniq_enrichment=NuclearOverlapHMM_names_melt %>% full_join(NuclearPermMean_uniq, by="number")%>% full_join(chromHmm, by="number")

NuclearOverlapHMM_Uniq_enrichment$NucQTL=as.double(NuclearOverlapHMM_Uniq_enrichment$NucQTL)

NuclearOverlapHMM_Uniq_enrichment=NuclearOverlapHMM_Uniq_enrichment %>%mutate(NucEnrich=(NucQTL-NucRandPerm)/NucRandPermSD)

ggplot(NuclearOverlapHMM_Uniq_enrichment, aes(y=NucEnrich, x=number, fill=name)) + geom_bar(stat="identity") + labs(title="Nuclear ApaQTL enrichment Z score", y="Enrichment Z score")

ggplot(TotalOverlapHMM_Uniq_enrichment, aes(y=TotEnrich, x=number, fill=name)) + geom_bar(stat="identity")+ labs(title="Total ApaQTL enrichment Z score", y="Enrichment Z score")

bothEnrich_uniq=NuclearOverlapHMM_Uniq_enrichment %>% full_join(TotalOverlapHMM_Uniq_enrichment, by=c("name", "number")) %>% select(number, name, NucEnrich,TotEnrich)

bothEnrich_uniqmelt=melt(bothEnrich_uniq, id.vars=c("number", "name"))

ggplot(bothEnrich_uniqmelt, aes(x=number, by=variable, fill=variable, y=value,col=name)) + geom_bar(position = "dodge", stat = "identity",alpha=.5) + scale_fill_manual(values=c("darkviolet", "deepskyblue3")) + labs(y="Enrichment from 1000 permutations", title="ChromHMM enrichment for \nTotal and Nuclear ApaQTLs",x="Region")