APAqtls with Leafcutter

module load samtools
#zip file 
gzip filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt 

module load python
#leafcutter script
python /project2/gilad/briana/threeprimeseq/code/prepare_phenotype_table.py filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz 

#source activate three-prime-env
sh filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz_prepare.sh

#run for nuclear as well 
gzip filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt 
#unload anaconda, load python
python /project2/gilad/briana/threeprimeseq/code/prepare_phenotype_table.py filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz 
#load anaconda and env. 
sh filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz_prepare.sh

#filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.PCs
#filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.PCs

#make a sample list  

fout = file("/project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt",'w')

for ln in open("/project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/file_id_mapping_nuc.txt", "r"):
    bam, sample = ln.split()
    line=sample[:-2]
    fout.write("NA"+line + "\n")
fout.close()

#!/bin/bash


#SBATCH --job-name=APAqtl_nominal_nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=APAqtl_nominal_nuc.out
#SBATCH --error=APAqtl_nominal_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 --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/nominal_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear.txt.gz.qqnorm_chr$i.nominal.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

#!/bin/bash


#SBATCH --job-name=APAqtl_nominal_tot
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=APAqtl_nominal_tot.out
#SBATCH --error=APAqtl_nominal_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 --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/nominal_APAqtl/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total.txt.gz.qqnorm_chr$i.nominal.out --chunk 1 1  --window 5e4 --include-samples /project2/gilad/briana/threeprimeseq/data/filt_peak_refGene_cov/SAMPLE.txt
done

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")

Expand here to see past versions of unnamed-chunk-12-1.png:

Version	Author	Date
b6e6ed9	brimittleman	2018-08-21

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")

Expand here to see past versions of unnamed-chunk-12-2.png:

Version	Author	Date
b6e6ed9	brimittleman	2018-08-21

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")

Expand here to see past versions of unnamed-chunk-13-1.png:

Version	Author	Date
bd21c34	brimittleman	2018-08-21

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

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

Expand here to see past versions of unnamed-chunk-13-2.png:

Version	Author	Date
bd21c34	brimittleman	2018-08-21

< .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")

Expand here to see past versions of unnamed-chunk-14-1.png:

Version	Author	Date
bd21c34	brimittleman	2018-08-21

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

Expand here to see past versions of unnamed-chunk-14-2.png:

Version	Author	Date
bd21c34	brimittleman	2018-08-21

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(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', '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'), stringsAsFactors = F) %>% select(one_of(samples))

chr15.76191353geno_anno=read.table("../data/perm_QTL/chrom15pos76191353.vcf", col.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'), stringsAsFactors = F) %>% select(CHROM, POS, sid, 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=as.integer(which(chr15.76191353geno_dose_full$POS==76191353))[1]
gen

[1] 84

#subset snps in window
tot.res_window=tot.res %>% semi_join(chr15.76191353geno_dose_full, by="sid")


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

start=76191353 - 25000
end=76191353 + 25000


#plot_locuszoom(mylist, 84, start, end)

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   workflowr_1.1.1 forcats_0.3.0  
 [5] stringr_1.3.1   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

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     knitr_1.20        broom_0.5.0      
[25] Rcpp_0.12.18      scales_0.5.0      backports_1.1.2  
[28] jsonlite_1.5      hms_0.4.2         digest_0.6.15    
[31] stringi_1.2.4     grid_3.5.1        rprojroot_1.3-2  
[34] cli_1.0.0         tools_3.5.1       magrittr_1.5     
[37] lazyeval_0.2.1    crayon_1.3.4      whisker_0.3-2    
[40] pkgconfig_2.0.1   xml2_1.2.0        lubridate_1.7.4  
[43] assertthat_0.2.0  rmarkdown_1.10    httr_1.3.1       
[46] rstudioapi_0.7    R6_2.2.2          nlme_3.1-137     
[49] git2r_0.23.0      compiler_3.5.1