APAqtls with Leafcutter

Filter dose files

I need to remove non snps and snps with <.05 from the dosage file.

I will first copy all of the dosage files to my direcory instead of changing tonys.

cp *dose.vcf.gz /project2/gilad/briana/YRI_geno_hg19/

I want to write a python script that will read in the files and perform the filters.

I wrote a python script that take in the dose file and a name of an out file. I will write a bash script to wrap this on all of the chrs.

#!/bin/bash


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

module load python

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 
do 
python filter_vcf.py chr$i.dose.vcf chr$i.dose.filt.vcf
done

Now I can use these for the fastqtl script instead.

I also updated to only use the first 2 pcs as covariates.

Run permuted version

Permutation pass to calculate correctedp-values for molecular phenotypes.

APAqtl_perm_tot.sh

#!/bin/bash


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

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/chr$i.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.2PCs --bed /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr$i.gz --out /project2/gilad/briana/threeprimeseq/data/perm_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr$i.perm.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

APAqtl_perm_nuc.sh

#!/bin/bash


#SBATCH --job-name=APAqtl_nominal_nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=APAqtl_perm_nuc.out
#SBATCH --error=APAqtl_perm_nuc.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END

for i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --vcf /project2/gilad/briana/YRI_geno_hg19/chr$i.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.2PCs --bed /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr$i.gz --out /project2/gilad/briana/threeprimeseq/data/perm_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr$i.perm.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

Try with normal approximation for the chroms that dont work:

APAqtl_perm_norm_tot.sh

#!/bin/bash


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

for i in 13 18 
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --normal --vcf /project2/gilad/briana/YRI_geno_hg19/chr$i.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.2PCs --bed /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr$i.gz --out /project2/gilad/briana/threeprimeseq/data/perm_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr$i.perm.norm.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

APAqtl_perm_norm_nuc.sh

#!/bin/bash


#SBATCH --job-name=APAqtl_nominal_nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=APAqtl_perm_nuc.out
#SBATCH --error=APAqtl_perm_nuc.err
#SBATCH --partition=broadwl
#SBATCH --mem=12G
#SBATCH --mail-type=END

for i in 3 13
do
/home/brimittleman/software/bin/FastQTL/bin/fastQTL.static --permute 1000 --normal --vcf /project2/gilad/briana/YRI_geno_hg19/chr$i.dose.filt.vcf.gz --cov /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.2PCs --bed /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr$i.gz --out /project2/gilad/briana/threeprimeseq/data/perm_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr$i.perm.norm.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

Evaluate the results

The results file has the folowing columns:

• ID of the tested molecular phenotype (in this particular case, the gene ID)
  
• Number of variants tested in cis for this phenotype
  
• MLE of the shape1 parameter of the Beta distribution
  
• MLE of the shape2 parameter of the Beta distribution
  
• Dummy [To be described later]
  
• ID of the best variant found for this molecular phenotypes (i.e. with the smallest p-value)
  
• Distance between the molecular phenotype - variant pair
  
• The nominal p-value of association that quantifies how significant from 0, the regression coefficient is
  
• The slope associated with the nominal p-value of association [only in version > v2-184]
  
• A first permutation p-value directly obtained from the permutations with the direct method. This is basically a corrected version of the nominal p-value that accounts for the fact that multiple variants are tested per molecular phenotype.
  
• A second permutation p-value obtained via beta approximation. We advice to use this one in any downstream analysis.

I can check the experiments as recomended by the FastQTL site.

d = read.table("permutations.all.chunks.txt.gz", hea=F, stringsAsFactors=F)
colnames(d) = c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "ppval", "bpval")
plot(d$ppval, d$bpval, xlab="Direct method", ylab="Beta approximation", main="Check plot")
abline(0, 1, col="red")

I will try this first on the resutls from chr1.

nuc.chr1= read.table("../data/perm_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr1.perm.out",head=F, stringsAsFactors=F, col.names = c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "slope", "ppval", "bpval"))


plot(nuc.chr1$ppval, nuc.chr1$bpval, xlab="Direct method", ylab="Beta approximation", main="Nuclear Check plot")
abline(0, 1, col="red")

tot.chr1=read.table("../data/perm_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr1.perm.out", head=F, stringsAsFactors = F, col.names= c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "slope", "ppval", "bpval"))


plot(tot.chr1$ppval, tot.chr1$bpval, xlab="Direct method", ylab="Beta approximation", main="Total Check plot")
abline(0, 1, col="red")

Correct for multiple testing:

• Bonferonni

nuc.chr1$bonferroni = p.adjust(nuc.chr1$bpval, method="bonferroni")

plot(-log10(nuc.chr1$bonferroni), main="Nuclear chr1 bonferroni corrected pval")

tot.chr1$bonferroni = p.adjust(tot.chr1$bpval, method="bonferroni")

plot(-log10(tot.chr1$bonferroni),  main="Total chr1 bonferroni corrected pval")

< .05 is 1.3 on this plot.

• BH

nuc.chr1$bh=p.adjust(nuc.chr1$bpval, method="fdr")

plot(-log10(nuc.chr1$bh), main="Nuclear chr1 BH corrected pval")

tot.chr1$bh=p.adjust(tot.chr1$bpval, method="fdr")
plot(-log10(tot.chr1$bh), main="Total chr1 BH corrected pval")

10% FDR is 1 on this plot.

Extend to all results:

nuc.res= read.table("../data/perm_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_permQTLresults.out",head=F, stringsAsFactors=F, col.names = c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "slope", "ppval", "bpval"))


plot(nuc.res$ppval, nuc.res$bpval, xlab="Direct method", ylab="Beta approximation", main="Nuclear Check plot")
abline(0, 1, col="red")

tot.res=read.table("../data/perm_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_permQTLresults.out", head=F, stringsAsFactors = F, col.names= c("pid", "nvar", "shape1", "shape2", "dummy", "sid", "dist", "npval", "slope", "ppval", "bpval"))


plot(tot.res$ppval, tot.res$bpval, xlab="Direct method", ylab="Beta approximation", main="Total Check plot")
abline(0, 1, col="red")

• BH

nuc.res$bh=p.adjust(nuc.res$bpval, method="fdr")

plot(-log10(nuc.res$bh), main="Nuclear BH corrected pval")
abline(h=1, col="red")

tot.res$bh=p.adjust(tot.res$bpval, method="fdr")
plot(-log10(tot.res$bh), main="Total BH corrected pval")
abline(h=1, col="red")

Next steps:

• make plots for some of these snps

• /project2/yangili1/yangili/APAqtl/output/ceu.apaqtl.txt.gz.bh.txt (use nominal pvalue)
  
• 1. plot a qqplot with only these SNPs
     
• 2. plot a qqplot with all SNPs that you tested

library(tidyverse)

── Attaching packages ──────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.0.0     ✔ purrr   0.2.5
✔ tibble  1.4.2     ✔ dplyr   0.7.6
✔ tidyr   0.8.1     ✔ stringr 1.3.1
✔ readr   1.1.1     ✔ forcats 0.3.0

── Conflicts ─────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(workflowr)

This is workflowr version 1.1.1
Run ?workflowr for help getting started

library(reshape2)

Attaching package: 'reshape2'

The following object is masked from 'package:tidyr':

    smiths

library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:ggplot2':

    ggsave

ceu_QTL=read.table("../data/nom_QTL/ceu.apaqtl.txt.gz.bh.txt", header = T, stringsAsFactors = F)
nom_nuc=read.table("../data/nom_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_nomQTLresults.out", head=F, stringsAsFactors = F, col.names = c("peakID", "snpID", "dist", "Nuc_pval", "slope"))
nom_tot=read.table("../data/nom_QTL/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_nomQTLresults.out",head=F , stringsAsFactors = F,  col.names = c("peakID", "snpID", "dist", "tot_pval", "slope"))

First I want to filter the CEU data just for snps. Then I want to reformat them to be in the same configuration as the nps in my results.

chr#:pos

ceu_QTL_snp=ceu_QTL %>% filter(grepl("snp", dummy2)) %>% separate(dummy2, c("type", "chr", "loc"), sep="_") %>% unite(snpID, c("chr", "loc"), sep=":")

Join the data frames by the snp ID.

ceuAndTot= ceu_QTL_snp %>% inner_join(nom_tot, by="snpID") %>% select(snpID, bpval, tot_pval)



ceuAndNuc= ceu_QTL_snp %>% inner_join(nom_nuc, by="snpID") %>% select(snpID, bpval, Nuc_pval)

tot_ceuSNPS=runif(nrow(ceuAndTot))

nuc_ceuSNPS=runif(nrow(ceuAndNuc))

par(mfrow=c(1,2))
qqplot(-log10(tot_ceuSNPS), -log10(ceuAndTot$tot_pval), ylab="-log10 Total pvalues", xlab="Uniform expectation", main="Total pvalues for in CEU snps")
abline(0,1)


qqplot(-log10(nuc_ceuSNPS), -log10(ceuAndNuc$Nuc_pva), ylab="-log10 Nuclear pvalues", xlab="Uniform expectation", main="Nuclear pvalues for in CEU snps")
abline(0,1)

Try with all of the snps:

par(mfrow=c(1,2))


qqplot(-log10(runif(nrow(nom_tot))), -log10(nom_tot$tot_pval), ylab="-log10 Total pvalue", xlab="Uniform expectation", main="Total pvalues for all snps")
abline(0,1)

qqplot(-log10(runif(nrow(nom_nuc))), -log10(nom_nuc$Nuc_pval), ylab="-log10 Nuclear pvalue", xlab="Uniform expectation",main= "Nuclear pvalues for all snps")
abline(0,1)

Try this with te permuted pvalues:

par(mfrow=c(1,2))
qqplot(-log10(runif(nrow(tot.res))), -log10(tot.res$bpval),ylab="-log10 Total permuted pvalue", xlab="Uniform expectation", main="Total permuted pvalues for all snps")
abline(0,1)

qqplot(-log10(runif(nrow(nuc.res))), -log10(nuc.res$bpval), ylab="-log10 Nuclear permuted pvalue", xlab="Uniform expectation", main="Nuclear permuted pvalues for all snps")
abline(0,1)

Locus zoom plots to vizualize the top QTLs:

Kenneth gave me this code for making these plots. I can modify this code.

plot_locuszoom <- function(this, gen, xlim, ylim, ...)
{
  
  #this is a r object that will have the results from the fastqtl and the genotypes 
  #this$annotations has gene, snp, dist, pvalue, beta, rsid, chr, pos, bpval, and other extra annotations about the snps
  rbPal <- colorRampPalette(c('lightblue','blue','purple','red'))(101)
  cols <- c()
  
  # gotta figure out how everythign correlates with this snp
  # row <- which(this$annotations$rsid==snp)
  # gen <- as.numeric(this$genotypes[row,10:129])
  nrow <- nrow(this$annotations)
  cors <- sapply(1:nrow, function(j) cor(gen, as.numeric(this$genotypes[j,10:33])))
  
  cols <- c()
  for (j in 1:nrow) cols[j] <- rbPal[round(100*(cors[j])^2)+1]
  
  plot.new()
  plot.window(xlim=xlim, ylim=ylim, xlab='position', ylab='-log10(p-value)', ...)

 
  points(x=this$annotations$pos, y=-log(this$annotations$bpval,10), pch=19, col=cols)
  axis(2)
  box()
  mtext('-log10(p-value)', side=2, line=2, cex=0.7)
}

I will try this with the top total snp first. It is in chrom15, the snip id is 15:76191353. I want to pull genotypes for snp within 50000 bases (window size).

I can write a python script that takes a snp position and filters only the snps 25000 up and 25000 downstream of this snp. I can subset just the individuals in the sample list once i move this into R.

Need to make sure to unzip the specfici vcf file first.

python filter_geno.py  15 76191353 /project2/gilad/briana/threeprimeseq/data/filtered_geno/chrom15pos76191353.vcf

samples=c("NA18486","NA18505", 'NA18508','NA18511','NA18519','NA18520','NA18853','NA18858','NA18861','NA18870','NA18909','NA18916','NA19119','NA19128','NA19130','NA19141','NA19160','NA19209','NA19210','NA19223','NA19225','NA19238','NA19239','NA19257')


chr15.76191353geno=read.table("../data/perm_QTL/chrom15pos76191353.vcf", col.names=c('CHROM', 'POS', 'snpID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'NA18486', 'NA18487', 'NA18488', 'NA18489', 'NA18498', 'NA18499', 'NA18501', 'NA18502', 'NA18504', 'NA18505', 'NA18507', 'NA18508', 'NA18510', 'NA18511', 'NA18516', 'NA18517', 'NA18519', 'NA18520', 'NA18522', 'NA18523', 'NA18852', 'NA18853', 'NA18855', 'NA18856', 'NA18858', 'NA18859', 'NA18861', 'NA18862', 'NA18867', 'NA18868', 'NA18870', 'NA18871', 'NA18873', 'NA18874', 'NA18907', 'NA18909', 'NA18910', 'NA18912', 'NA18913', 'NA18916', 'NA18917', 'NA18923', 'NA18924', 'NA18933', 'NA18934', 'NA19093', 'NA19095', 'NA19096', 'NA19098', 'NA19099', 'NA19101', 'NA19102', 'NA19107', 'NA19108', 'NA19113', 'NA19114', 'NA19116', 'NA19117', 'NA19118', 'NA19119', 'NA19121', 'NA19122', 'NA19127', 'NA19128', 'NA19129', 'NA19130', 'NA19131', 'NA19137', 'NA19138', 'NA19140', 'NA19141', 'NA19143', 'NA19144', 'NA19146', 'NA19147', 'NA19149', 'NA19150', 'NA19152', 'NA19153', 'NA19159', 'NA19160', 'NA19171', 'NA19172', 'NA19175', 'NA19176', 'NA19184', 'NA19185', 'NA19189', 'NA19190', 'NA19197', 'NA19198', 'NA19200', 'NA19201', 'NA19203', 'NA19204', 'NA19206', 'NA19207', 'NA19209', 'NA19210', 'NA19213', 'NA19214', 'NA19222', 'NA19223', 'NA19225', 'NA19226', 'NA19235', 'NA19236', 'NA19238', 'NA19239', 'NA19247', 'NA19248', 'NA19256', 'NA19257'), stringsAsFactors = F) %>% select(one_of(samples))

chr15.76191353geno_anno=read.table("../data/perm_QTL/chrom15pos76191353.vcf", col.names=c('CHROM', 'POS', 'snpID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'NA18486', 'NA18487', 'NA18488', 'NA18489', 'NA18498', 'NA18499', 'NA18501', 'NA18502', 'NA18504', 'NA18505', 'NA18507', 'NA18508', 'NA18510', 'NA18511', 'NA18516', 'NA18517', 'NA18519', 'NA18520', 'NA18522', 'NA18523', 'NA18852', 'NA18853', 'NA18855', 'NA18856', 'NA18858', 'NA18859', 'NA18861', 'NA18862', 'NA18867', 'NA18868', 'NA18870', 'NA18871', 'NA18873', 'NA18874', 'NA18907', 'NA18909', 'NA18910', 'NA18912', 'NA18913', 'NA18916', 'NA18917', 'NA18923', 'NA18924', 'NA18933', 'NA18934', 'NA19093', 'NA19095', 'NA19096', 'NA19098', 'NA19099', 'NA19101', 'NA19102', 'NA19107', 'NA19108', 'NA19113', 'NA19114', 'NA19116', 'NA19117', 'NA19118', 'NA19119', 'NA19121', 'NA19122', 'NA19127', 'NA19128', 'NA19129', 'NA19130', 'NA19131', 'NA19137', 'NA19138', 'NA19140', 'NA19141', 'NA19143', 'NA19144', 'NA19146', 'NA19147', 'NA19149', 'NA19150', 'NA19152', 'NA19153', 'NA19159', 'NA19160', 'NA19171', 'NA19172', 'NA19175', 'NA19176', 'NA19184', 'NA19185', 'NA19189', 'NA19190', 'NA19197', 'NA19198', 'NA19200', 'NA19201', 'NA19203', 'NA19204', 'NA19206', 'NA19207', 'NA19209', 'NA19210', 'NA19213', 'NA19214', 'NA19222', 'NA19223', 'NA19225', 'NA19226', 'NA19235', 'NA19236', 'NA19238', 'NA19239', 'NA19247', 'NA19248', 'NA19256', 'NA19257'), stringsAsFactors = F) %>% select(CHROM, POS, snpID, REF, ALT, QUAL, FILTER, INFO, FORMAT)
chr15.76191353geno_dose=apply(chr15.76191353geno, 2, function(y)sapply(y, function(x)as.integer(strsplit(x,":")[[1]][[2]])))

chr15.76191353geno_dose_full=data.frame(cbind(chr15.76191353geno_anno, chr15.76191353geno_dose))

gen=chr15.76191353geno_dose_full[which(chr15.76191353geno_dose_full$POS==76191353),]
gen

   CHROM      POS       snpID REF ALT QUAL FILTER
84    15 76191353 15:76191353   C   T    .   PASS
                                INFO   FORMAT NA18486 NA18505 NA18508
84 AF=0.08407;MAF=0.08407;R2=0.99998 GT:DS:GP       0       0       1
   NA18511 NA18519 NA18520 NA18853 NA18858 NA18861 NA18870 NA18909 NA18916
84       0       0       0       1       0       0       1       0       0
   NA19119 NA19128 NA19130 NA19141 NA19160 NA19209 NA19210 NA19223 NA19225
84       0       0       0       0       0       0       0       0       0
   NA19238 NA19239 NA19257
84       0       0       0

snps=chr15.76191353geno_dose_full$snpID
in_both_nom= nom_tot %>% filter(snpID %in% snps)


mylist=list(annotations=tot.res,genotypes=chr15.76191353geno_dose_full )

start=76191353 - 25000
end=76191353 + 25000


#plot_locuszoom(mylist, gen, start, end)

I actually need to do this with the nominal snps.

The most sig. in the nominal total is 4:186328829:186328922:NM_018359.3_-_peak260565, 4:186325141

I want to run the python genotype filter.

python filter_geno.py  4 186325141 /project2/gilad/briana/threeprimeseq/data/filtered_geno/chrom4pos186325141.vcf

chrom4pos18632514=read.table("../data/nom_QTL/chrom4pos186325141.vcf", col.names=c('CHROM', 'POS', 'snpID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'NA18486', 'NA18487', 'NA18488', 'NA18489', 'NA18498', 'NA18499', 'NA18501', 'NA18502', 'NA18504', 'NA18505', 'NA18507', 'NA18508', 'NA18510', 'NA18511', 'NA18516', 'NA18517', 'NA18519', 'NA18520', 'NA18522', 'NA18523', 'NA18852', 'NA18853', 'NA18855', 'NA18856', 'NA18858', 'NA18859', 'NA18861', 'NA18862', 'NA18867', 'NA18868', 'NA18870', 'NA18871', 'NA18873', 'NA18874', 'NA18907', 'NA18909', 'NA18910', 'NA18912', 'NA18913', 'NA18916', 'NA18917', 'NA18923', 'NA18924', 'NA18933', 'NA18934', 'NA19093', 'NA19095', 'NA19096', 'NA19098', 'NA19099', 'NA19101', 'NA19102', 'NA19107', 'NA19108', 'NA19113', 'NA19114', 'NA19116', 'NA19117', 'NA19118', 'NA19119', 'NA19121', 'NA19122', 'NA19127', 'NA19128', 'NA19129', 'NA19130', 'NA19131', 'NA19137', 'NA19138', 'NA19140', 'NA19141', 'NA19143', 'NA19144', 'NA19146', 'NA19147', 'NA19149', 'NA19150', 'NA19152', 'NA19153', 'NA19159', 'NA19160', 'NA19171', 'NA19172', 'NA19175', 'NA19176', 'NA19184', 'NA19185', 'NA19189', 'NA19190', 'NA19197', 'NA19198', 'NA19200', 'NA19201', 'NA19203', 'NA19204', 'NA19206', 'NA19207', 'NA19209', 'NA19210', 'NA19213', 'NA19214', 'NA19222', 'NA19223', 'NA19225', 'NA19226', 'NA19235', 'NA19236', 'NA19238', 'NA19239', 'NA19247', 'NA19248', 'NA19256', 'NA19257'), stringsAsFactors = F) %>% select(one_of(samples))

chrom4pos18632514_anno=read.table("../data/nom_QTL/chrom4pos186325141.vcf", col.names=c('CHROM', 'POS', 'snpID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'NA18486', 'NA18487', 'NA18488', 'NA18489', 'NA18498', 'NA18499', 'NA18501', 'NA18502', 'NA18504', 'NA18505', 'NA18507', 'NA18508', 'NA18510', 'NA18511', 'NA18516', 'NA18517', 'NA18519', 'NA18520', 'NA18522', 'NA18523', 'NA18852', 'NA18853', 'NA18855', 'NA18856', 'NA18858', 'NA18859', 'NA18861', 'NA18862', 'NA18867', 'NA18868', 'NA18870', 'NA18871', 'NA18873', 'NA18874', 'NA18907', 'NA18909', 'NA18910', 'NA18912', 'NA18913', 'NA18916', 'NA18917', 'NA18923', 'NA18924', 'NA18933', 'NA18934', 'NA19093', 'NA19095', 'NA19096', 'NA19098', 'NA19099', 'NA19101', 'NA19102', 'NA19107', 'NA19108', 'NA19113', 'NA19114', 'NA19116', 'NA19117', 'NA19118', 'NA19119', 'NA19121', 'NA19122', 'NA19127', 'NA19128', 'NA19129', 'NA19130', 'NA19131', 'NA19137', 'NA19138', 'NA19140', 'NA19141', 'NA19143', 'NA19144', 'NA19146', 'NA19147', 'NA19149', 'NA19150', 'NA19152', 'NA19153', 'NA19159', 'NA19160', 'NA19171', 'NA19172', 'NA19175', 'NA19176', 'NA19184', 'NA19185', 'NA19189', 'NA19190', 'NA19197', 'NA19198', 'NA19200', 'NA19201', 'NA19203', 'NA19204', 'NA19206', 'NA19207', 'NA19209', 'NA19210', 'NA19213', 'NA19214', 'NA19222', 'NA19223', 'NA19225', 'NA19226', 'NA19235', 'NA19236', 'NA19238', 'NA19239', 'NA19247', 'NA19248', 'NA19256', 'NA19257'), stringsAsFactors = F) %>% select(CHROM, POS, snpID, REF, ALT, QUAL, FILTER, INFO, FORMAT)
chrom4pos18632514_dose=apply(chrom4pos18632514, 2, function(y)sapply(y, function(x)as.integer(strsplit(x,":")[[1]][[2]])))

chrom4pos18632514_dose_full=data.frame(cbind(chrom4pos18632514_anno, chrom4pos18632514_dose))


snps=chrom4pos18632514_dose_full$snpID
in_both_nom= nom_tot %>% filter(snpID %in% snps)


gen=chrom4pos18632514_dose_full[which(chrom4pos18632514_dose_full$POS==186325141),]


mylist=list(annotations=in_both_nom,genotypes=chrom4pos18632514_dose_full)


start=18632514- 25000
end=18632514 + 25000

#plot_locuszoom(mylist, gen, start, end)

#problem: the in_both_nom has more values because snps can be associated with more than one peak  w

Try to make a boxplot:

FIrst for the strongest total pval.

geno=chr15.76191353geno_dose_full[which(chr15.76191353geno_dose_full$POS==76191353),10:33]
# find the phentpye values for peak 15:76234771:76234852:NM_138573.3_-_peak118132
#grep -F "15:76234771:76234852:NM_138573.3_-_peak118132" filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.phen_chr15 > ../qtl_example/tot_peak118132
pheno=read.table("../data/perm_QTL/tot_peak118132", stringsAsFactors = F, col.names = c("Chr",  "start",    "end",  "ID",   'NA18486',  'NA18497',  'NA18500',  'NA18505','NA18508' ,'NA18511', 'NA18519',  'NA18520',  'NA18853',  'NA18858',  'NA18861'   ,'NA18870', 'NA18909',  'NA18916',  'NA19092',  'NA19119',  'NA19128'   ,'NA19130', 'NA19141'   ,'NA19160', 'NA19193',  'NA19209'   ,'NA19210', 'NA19223'   ,'NA19225', 'NA19238',  'NA19239'   ,'NA19257')) %>%  select(one_of(samples))


for_plot=data.frame(bind_rows(geno,pheno) %>% t)
colnames(for_plot)=c("Genotype", "PAS")
for_plot$Genotype=as.factor(for_plot$Genotype)


ggplot(for_plot, aes(x=Genotype, y=PAS, fill=Genotype, group=Genotype)) + geom_boxplot() + labs(x="Genotype", title="15:76234771:76234852:NM_138573.3_-_peak118132 QTL") + geom_jitter( aes(x=Genotype, y=PAS))

Generally I will need to grep the correct line from the geno and pheno file then make the plot like this.

Next I will run for the top Nuc QTL.

peak: 12:9092958:9093051:NM_004426.2_+_peak67056 SNP: 12:9049821

#grep -F "12:9092958:9093051:NM_004426.2_+_peak67056" filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.phen_chr12 > ../qtl_example/nuc_peak67056
pheno_names=c("Chr",    "start",    "end",  "ID",   'NA18486',  'NA18497',  'NA18500',  'NA18505','NA18508' ,'NA18511', 'NA18519',  'NA18520',  'NA18853',  'NA18858',  'NA18861'   ,'NA18870', 'NA18909',  'NA18916',  'NA19092',  'NA19119',  'NA19128'   ,'NA19130', 'NA19141'   ,'NA19160', 'NA19193',  'NA19209'   ,'NA19210', 'NA19223'   ,'NA19225', 'NA19238',  'NA19239'   ,'NA19257')
geno_names=c('CHROM', 'POS', 'sid', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO', 'FORMAT', 'NA18486', 'NA18487', 'NA18488', 'NA18489', 'NA18498', 'NA18499', 'NA18501', 'NA18502', 'NA18504', 'NA18505', 'NA18507', 'NA18508', 'NA18510', 'NA18511', 'NA18516', 'NA18517', 'NA18519', 'NA18520', 'NA18522', 'NA18523', 'NA18852', 'NA18853', 'NA18855', 'NA18856', 'NA18858', 'NA18859', 'NA18861', 'NA18862', 'NA18867', 'NA18868', 'NA18870', 'NA18871', 'NA18873', 'NA18874', 'NA18907', 'NA18909', 'NA18910', 'NA18912', 'NA18913', 'NA18916', 'NA18917', 'NA18923', 'NA18924', 'NA18933', 'NA18934', 'NA19093', 'NA19095', 'NA19096', 'NA19098', 'NA19099', 'NA19101', 'NA19102', 'NA19107', 'NA19108', 'NA19113', 'NA19114', 'NA19116', 'NA19117', 'NA19118', 'NA19119', 'NA19121', 'NA19122', 'NA19127', 'NA19128', 'NA19129', 'NA19130', 'NA19131', 'NA19137', 'NA19138', 'NA19140', 'NA19141', 'NA19143', 'NA19144', 'NA19146', 'NA19147', 'NA19149', 'NA19150', 'NA19152', 'NA19153', 'NA19159', 'NA19160', 'NA19171', 'NA19172', 'NA19175', 'NA19176', 'NA19184', 'NA19185', 'NA19189', 'NA19190', 'NA19197', 'NA19198', 'NA19200', 'NA19201', 'NA19203', 'NA19204', 'NA19206', 'NA19207', 'NA19209', 'NA19210', 'NA19213', 'NA19214', 'NA19222', 'NA19223', 'NA19225', 'NA19226', 'NA19235', 'NA19236', 'NA19238', 'NA19239', 'NA19247', 'NA19248', 'NA19256', 'NA19257')

top_nuc_geno=read.table("../data/perm_QTL/genotpye12:904921", stringsAsFactors = F, col.names = geno_names) %>%  select(one_of(samples))

top_nuc_geno_dose=apply(top_nuc_geno, 2, function(y)sapply(y, function(x)as.integer(strsplit(x,":")[[1]][[2]])))


top_nuc_pheo=read.table("../data/perm_QTL/nuc_peak67056", stringsAsFactors = F, col.names = pheno_names) %>% select(one_of(samples))



top_nuc_plot=data.frame(bind_rows(top_nuc_geno_dose, top_nuc_pheo) %>% t)
colnames(top_nuc_plot)=c("Genotype", "PAS")
top_nuc_plot$Genotype=as.factor(top_nuc_plot$Genotype)

ggplot(top_nuc_plot, aes(x=Genotype, y=PAS, fill=Genotype, group=Genotype)) + geom_boxplot() + labs(x="Genotype", title="12:9092958:9093051:NM_004426.2_+_peak67056") + geom_jitter( aes(x=Genotype, y=PAS))

3:119242427:119242509:NM_016589.3_+_peak233134

3:119211867

grep -F "3:119242427:119242509:NM_016589.3_+_peak233134" filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.phen_chr3 > ../qtl_example/nuc_peak233134

YRI_geno_hg19]$ less chr3.dose.filt.vcf.gz | grep   "3:119211867" > ../threeprimeseq/data/qtl_example/genotype3:199211867

top_nuc_geno2=read.table("../data/perm_QTL/genotype3:199211867", stringsAsFactors = F, col.names = geno_names) %>%  select(one_of(samples))

top_nuc_geno2_dose=apply(top_nuc_geno2, 2, function(y)sapply(y, function(x)as.integer(strsplit(x,":")[[1]][[2]])))


top_nuc_pheo2=read.table("../data/perm_QTL/nuc_peak233134", stringsAsFactors = F, col.names = pheno_names) %>% select(one_of(samples))



top_nuc_plot2=data.frame(bind_rows(top_nuc_geno2_dose, top_nuc_pheo2) %>% t)
colnames(top_nuc_plot2)=c("Genotype", "PAS")
top_nuc_plot2$Genotype=as.factor(top_nuc_plot2$Genotype)

ggplot(top_nuc_plot2, aes(x=Genotype, y=PAS, fill=Genotype, group=Genotype)) + geom_boxplot() + labs(x="Genotype", title="12:9092958:9093051:NM_004426.2_+_peak67056") + geom_jitter( aes(x=Genotype, y=PAS))

tot_QTL=tot.res %>% filter(bh < .15 )
nuc_QTL= nuc.res %>% filter(bh< .15)

tot.res = tot.res %>% mutate(QTL=ifelse(bh<.15, "Yes", "No") )
nuc.res = nuc.res %>% mutate(QTL=ifelse(bh<.15, "Yes", "No") )

tot.dist=ggplot(tot.res, aes(x=log10(abs(dist)+1), group=QTL, fill=QTL)) + geom_density(alpha=.4) + labs(title="Distribtuion of density in Total QTLS", x="Log 10 abs. values distance from SNP to peaks")
nuc.dist=ggplot(nuc.res, aes(x=log10(abs(dist)+1), group=QTL, fill=QTL)) + geom_density(alpha=.4) + labs(title="Distribtuion of density in Nuclear QTLS",x="Log 10 abs. values distance from SNP to peaks")
plot_grid(tot.dist, nuc.dist)

nQTL_tot=c()
FDR=seq(.05, .5, .01)
for (i in FDR){
  x=tot.res %>% filter(bh < i ) %>% nrow()
  nQTL_tot=c(nQTL_tot, x)
}

FDR=seq(.05, .5, .01)
nQTL_nuc=c()
for (i in FDR){
  x=nuc.res %>% filter(bh < i ) %>% nrow()
  nQTL_nuc=c(nQTL_nuc, x)
}

nQTL=as.data.frame(cbind(FDR, Total=nQTL_tot, Nuclear=nQTL_nuc))
nQTL_long=melt(nQTL, id.vars = "FDR")

ggplot(nQTL_long, aes(x=FDR, y=value, by=variable, col=variable)) + geom_line(size=1.5) + labs(y="Number of Significant QTLs", title="APAqtls detected by FDR cuttoff", color="Fraction") 

Session information

sessionInfo()

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

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  cowplot_0.9.3   reshape2_1.4.3  workflowr_1.1.1
 [5] forcats_0.3.0   stringr_1.3.1   dplyr_0.7.6     purrr_0.2.5    
 [9] readr_1.1.1     tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0  
[13] tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4  haven_1.1.2       lattice_0.20-35  
 [4] colorspace_1.3-2  htmltools_0.3.6   yaml_2.1.19      
 [7] rlang_0.2.1       R.oo_1.22.0       pillar_1.3.0     
[10] glue_1.3.0        withr_2.1.2       R.utils_2.6.0    
[13] modelr_0.1.2      readxl_1.1.0      bindr_0.1.1      
[16] plyr_1.8.4        munsell_0.5.0     gtable_0.2.0     
[19] cellranger_1.1.0  rvest_0.3.2       R.methodsS3_1.7.1
[22] evaluate_0.11     labeling_0.3      knitr_1.20       
[25] broom_0.5.0       Rcpp_0.12.18      scales_0.5.0     
[28] backports_1.1.2   jsonlite_1.5      hms_0.4.2        
[31] digest_0.6.15     stringi_1.2.4     grid_3.5.1       
[34] rprojroot_1.3-2   cli_1.0.0         tools_3.5.1      
[37] magrittr_1.5      lazyeval_0.2.1    crayon_1.3.4     
[40] whisker_0.3-2     pkgconfig_2.0.1   xml2_1.2.0       
[43] lubridate_1.7.4   assertthat_0.2.0  rmarkdown_1.10   
[46] httr_1.3.1        rstudioapi_0.7    R6_2.2.2         
[49] nlme_3.1-137      git2r_0.23.0      compiler_3.5.1