New QTL heatmap

Compare QTLs to those found with previous batch data

I have about double the QTLs hear compared to before resequencing batch 4. I will look at the new QTL to see if there is evidence for them being false positives. I am going to see if there is structure in the genotypes for these QTLs.

The old QTLs are from the threeprimeseq repository.

oldtot=read.table("../../threeprimeseq/data/perm_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Total.fixed.pheno_5perc_permResBH.txt", header=T,stringsAsFactors = F) %>% separate(pid, into=c("Chr", "Start", "End", "PeakID"), sep=":") %>% separate(PeakID, into=c("Gene", "Strand","Peak"), sep="_")

Warning: Expected 3 pieces. Additional pieces discarded in 3 rows [886,
887, 888].

OldTotQTLs= oldtot %>% filter(-log10(bh)>=1)
nrow(OldTotQTLs)

[1] 291

newTotQTLs=read.table("../data/apaQTLs/Total_apaQTLs_5fdr.txt", stringsAsFactors = F, header = T)
nrow(newTotQTLs)

[1] 502

UniqueNewTot=newTotQTLs %>% anti_join(OldTotQTLs, by=c("sid","Gene"))

write.table(UniqueNewTot, file="../data/apaQTLs/Total_apaQTLs_5fdr_NewUnique.txt", quote = F, col.names = F, row.names = F)

oldnuc=read.table("../../threeprimeseq/data/perm_APAqtl_GeneLocAnno_noMP_5percUs/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Nuclear.fixed.pheno_5perc_permResBH.txt", header=T,stringsAsFactors = F) %>% separate(pid, into=c("Chr", "Start", "End", "PeakID"), sep=":") %>% separate(PeakID, into=c("Gene", "Strand","Peak"), sep="_")

Warning: Expected 3 pieces. Additional pieces discarded in 3 rows [1056,
1057, 1058].

OldNucQTLs= oldnuc %>% filter(-log10(bh)>=1)
nrow(OldNucQTLs)

[1] 615

newNucQTLs=read.table("../data/apaQTLs/Nuclear_apaQTLs_5fdr.txt", stringsAsFactors = F, header = T)
nrow(newNucQTLs)

[1] 1070

UniqueNewNuc=newNucQTLs %>% anti_join(OldNucQTLs, by=c("sid","Gene"))

write.table(UniqueNewNuc, file="../data/apaQTLs/Nuclear_apaQTLs_5fdr_NewUnique.txt", quote = F, col.names = F, row.names = F)

Extract genotypes:

I wrote a script to pull the doses from the vcf file. Run it with:

 python extractGenotypes.py ../data/apaQTLs/Nuclear_apaQTLs_5fdr_NewUnique.txt ../data/QTLGenotypes/Genotypes_NuclearapaQTLS_newunique.txt
 
  python extractGenotypes.py ../data/apaQTLs/Total_apaQTLs_5fdr_NewUnique.txt ../data/QTLGenotypes/Genotypes_TotalapaQTLS_newunique.txt
  

I also need the header from the VCF to have the individuals:

head -n14 /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf | tail -n1  > ../data/QTLGenotypes/vcfheader.txt

#manually remove # and unneaded columns, keep snp and ind. 

vcfhead=read.table("../data/QTLGenotypes/vcfheader.txt", header = T)

input sample list:

samples=read.table("../data/phenotype/SAMPLE.txt")
samplist=as.vector(samples$V1)

totgeno=read.table("../data/QTLGenotypes/Genotypes_TotalapaQTLS_newunique.txt", col.names = colnames(vcfhead)) %>% select(samplist) %>% t()

totgenosm=totgeno[,1:250]
totgeneCorr=round(cor(totgenosm),2)

heatmap.2(as.matrix(totgeneCorr),trace="none", dendrogram =c("none"), main="Genotype correlation\n for new Total QTL snps")

nucgeno=read.table("../data/QTLGenotypes/Genotypes_NuclearapaQTLS_newunique.txt", col.names = colnames(vcfhead)) %>% select(samplist) %>% t()

nucgenosm=nucgeno[,1:600]
nucgeneCorr=round(cor(nucgenosm),2)
heatmap.2(as.matrix(nucgeneCorr),trace="none", dendrogram =c("none"),main="Genotype correlation \n for new Nuclear QTL snps")

Structure in all QTLs

This was for the new QTLs. I want to see if there is structure more generally.

 python extractGenotypes.py ../data/apaQTLs/Nuclear_apaQTLs_5fdr.txt ../data/QTLGenotypes/Genotypes_NuclearapaQTLS_ALLQTLs.txt

 python extractGenotypes.py ../data/apaQTLs/Total_apaQTLs_5fdr.txt ../data/QTLGenotypes/Genotypes_TotalapaQTLS_ALLQTLs.txt

totgenoall=read.table("../data/QTLGenotypes/Genotypes_TotalapaQTLS_ALLQTLs.txt", col.names = colnames(vcfhead)) %>% select(samplist) %>% t()

totgenoallsm=totgenoall[,1:300]

totgeneallCorr=round(cor(totgenoallsm),2)

heatmap.2(as.matrix(totgeneallCorr),trace="none", dendrogram =c("none"),main="Genotype correlation \n for first 300 Total QTL snps")

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

nucgenoall=read.table("../data/QTLGenotypes/Genotypes_NuclearapaQTLS_ALLQTLs.txt", col.names = colnames(vcfhead)) %>% select(samplist) %>% t()
nucgenoallsm=nucgenoall[,1:400]


nucgeneallCorr=round(cor(nucgenoallsm),2)

heatmap.2(as.matrix(nucgeneallCorr),trace="none", dendrogram =c("none"),main="Genotype correlation \n for first 400 Nuclear QTL snps")

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

Compare beta values in 55 vs 39

I want to make a scatter plot comaparring the new QTL associations in the 55 vs 39 individauls. If the qtls are real we expect a high correlation.

To do this I can recall the qtls with a smaller sample list excluding the 15 new individuals.

I need to make a list of the individuals not in the 4th batch.

batch1.2.3=read.table("../data/MetaDataSequencing.txt", header=T,stringsAsFactors = F)%>% filter(fraction=="total") %>%  select(line, batch) %>% filter(batch != 4)

samplelist=read.table("../data/phenotype/SAMPLE.txt", col.names = c("line"),stringsAsFactors = F)

Make a new directory for the 39ind qtls:

mkdir ../data/ThirtyNineIndQtl_nominal

Filter the sample list

samplelist_39= samplelist %>% semi_join(batch1.2.3, by="line")

write.table(samplelist_39, file="../data/ThirtyNineIndQtl_nominal/samplelist39.txt", col.names = F, row.names = F, quote = F)

Run the QTL code with this sample list

sbatch aAPAqtl_nominal39ind.sh

Concatinate results:

cat APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz.qqnorm_chr* > APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz.qqnorm_AllChr.txt

cat APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz.qqnorm_chr* > APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz.qqnorm_AllChr.txt

I want to filter the results for the new snps in the uniquenewtot. These results are in data/apaQTLs

I need to write a script that makes a dictionary with each of the new QTLs in the format above. Then I can run throguh the nominal values and keep only the values in the dictionary.

I can run this on the 55 and 39 nominal files then combine the files to create the scatterplot.

python selectNominalPvalues.py ../data/apaQTLs/Total_apaQTLs_5fdr_NewUnique.txt ../data/ThirtyNineIndQtl_nominal/APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz.qqnorm_AllChr.txt ../data/ThirtyNineIndQtl_nominal/Total_apaQTLs_NewUniqNom_37ing.txt

python selectNominalPvalues.py ../data/apaQTLs/Total_apaQTLs_5fdr_NewUnique.txt ../data/apaQTLNominal/APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz.qqnorm_AllChr.txt ../data/ThirtyNineIndQtl_nominal/Total_apaQTLs_NewUniqNom_55ind.txt

newin37_tot=read.table("../data/ThirtyNineIndQtl_nominal/Total_apaQTLs_NewUniqNom_37ing.txt",col.names=c("peakID", "snp", "dist", "Nompval39","Beta39"), stringsAsFactors = F)

newin54_tot=read.table("../data/ThirtyNineIndQtl_nominal/Total_apaQTLs_NewUniqNom_55ind.txt",col.names=c("peakID", "snp", "dist", "Nompval54","Beta54"), stringsAsFactors = F) 

newinboth=newin54_tot %>% inner_join(newin37_tot, by=c("peakID", "snp"))

total_qtlind=ggplot(newinboth,aes(x=Beta54, y=Beta39)) + geom_point()  + labs(title="New Total apaQTLs \nin different ind. sets", ylab="Beta 39 ind", xlab="Beta 55 ind")
total_qtlind

python selectNominalPvalues.py ../data/apaQTLs/Nuclear_apaQTLs_5fdr_NewUnique.txt ../data/ThirtyNineIndQtl_nominal/APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz.qqnorm_AllChr.txt ../data/ThirtyNineIndQtl_nominal/Nuclear_apaQTLs_NewUniqNom_37ing.txt

python selectNominalPvalues.py ../data/apaQTLs/Nuclear_apaQTLs_5fdr_NewUnique.txt ../data/apaQTLNominal/APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz.qqnorm_AllChr.txt ../data/ThirtyNineIndQtl_nominal/Nuclear_apaQTLs_NewUniqNom_55ind.txt

newin37_nuc=read.table("../data/ThirtyNineIndQtl_nominal/Nuclear_apaQTLs_NewUniqNom_37ing.txt",col.names=c("peakID", "snp", "dist", "Nompval39","Beta39"), stringsAsFactors = F)

newin54_nuc=read.table("../data/ThirtyNineIndQtl_nominal/Nuclear_apaQTLs_NewUniqNom_55ind.txt",col.names=c("peakID", "snp", "dist", "Nompval54","Beta54"), stringsAsFactors = F)

newinboth_nuc=newin54_nuc %>% inner_join(newin37_nuc, by=c("peakID", "snp"))

nuclear_qtlind=ggplot(newinboth_nuc,aes(x=Beta54, y=Beta39)) + geom_point()  + labs(title="New Nuclear apaQTLs\n in different ind. sets", ylab="Beta 39 ind", xlab="Beta 55 ind")

plot_grid(total_qtlind, nuclear_qtlind)

pvalplot_tot=ggplot(newinboth,aes(x=-log10(Nompval54), y=-log10(Nompval39))) + geom_point()  + labs(title="New Total apaQTLs \nin different ind. sets", y="-log10 pval 39 ind", x="-log10 pval 55 ind")

pvalplot_tot

pvalplot_nuc=ggplot(newinboth_nuc,aes(x=-log10(Nompval54), y=-log10(Nompval39))) + geom_point()  + labs(title="New Nuclear apaQTLs \nin different ind. sets", y="-log10 pval 39 ind", x="-log10 pval 55 ind")

plot_grid(pvalplot_tot,pvalplot_nuc)

newvold=read.table("../data/peaks_5perc/NewVOldPeaks.txt", header = T, stringsAsFactors = F) %>% select(peak, New)

colnames(newvold)=c("Peak", "Set")

UniqueNewTot_set= UniqueNewTot %>% inner_join(newvold, by="Peak")


UniqueNewTot_set$Set=as.factor(UniqueNewTot_set$Set)
summary(UniqueNewTot_set$Set)

     new original 
     209      188 

UniqueNewNuc_set= UniqueNewNuc %>% inner_join(newvold, by="Peak")


UniqueNewNuc_set$Set=as.factor(UniqueNewNuc_set$Set)
summary(UniqueNewNuc_set$Set)

     new original 
     521      349 

Example boxplots

I want to plot these qtls as boxplots. I will plot usage vs. genotype. I want to color the new 15 individuals differently.

peak15390 CREM 10:35411246

less /project2/gilad/briana/YRI_geno_hg19/chr10.dose.filt.vcf.gz  |  grep 10:35411246 > ../data/exampleQTLs/Geno_10:35411246.txt



#get pheno from pheno_5perc
less ../data/phenotype_5perc/APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz | grep peak15390 > ../data/exampleQTLs/Pheno_peak15390.txt

col.names = colnames(vcfhead))

head -n14 /project2/gilad/briana/YRI_geno_hg19/allChrom.dose.filt.vcf | tail -n1  > ../data/exampleQTLs/vcfheader.txt

#get rid of # manually 

#pheno head

 less APApeak_Phenotype_GeneLocAnno.Total.5perc.fc.gz | head -n1 > ../exampleQTLs/phenoHead.txt

vcfheadfull=read.table("../data/exampleQTLs/vcfheader.txt", header = T)
phenohead=read.table("../data/exampleQTLs/phenoHead.txt", header = T)

Try first

geno_peak15390=read.table("../data/exampleQTLs/Geno_10:35411246.txt", col.names = colnames(vcfheadfull), stringsAsFactors = F) %>% select(-CHROM, -POS, -REF, -ALT, -QUAL, -FILTER, -INFO, -FORMAT)
geno_peak15390M=melt(geno_peak15390, id.vars = "ID") %>% separate(value, into=c("genotype", "extra", "extra2"), sep=":") %>% mutate(dose=round(as.integer(extra), digits = 1)) %>% select(-genotype, -extra, -extra2)

geno_peak15390M$variable=as.character(geno_peak15390M$variable)

pheno_peak15390=read.table("../data/exampleQTLs/Pheno_peak15390.txt", col.names = colnames(phenohead), stringsAsFactors = F) 
pheno_peak15390M=melt(pheno_peak15390,id.vars = "chrom") %>% separate(value, into=c("num", "dom"), sep="/") %>% mutate(usage=as.integer(num)/as.integer(dom))

pheno_peak15390M$variable=as.character(pheno_peak15390M$variable)

Join these

genoPhenopeak15390=pheno_peak15390M %>% inner_join(geno_peak15390M, by="variable") %>% mutate(set=ifelse(variable %in% samplelist_39$line, "old39", "new15"))

genoPhenopeak15390$set=as.factor(genoPhenopeak15390$set)
genoPhenopeak15390$dose=as.factor(genoPhenopeak15390$dose)

ggplot(genoPhenopeak15390, aes(x=dose, y=usage)) + geom_boxplot()+ geom_jitter(aes(col=set))

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

Try 1 more: PAFAH1B2 peak25812 11:117048417

#get genotype from vcf
less /project2/gilad/briana/YRI_geno_hg19/chr11.dose.filt.vcf.gz  |  grep 11:117048417 > ../data/exampleQTLs/Geno_11:117048417.txt

#get pheno from pheno_5perc
less ../data/phenotype_5perc/APApeak_Phenotype_GeneLocAnno.Nuclear.5perc.fc.gz | grep peak25812 > ../data/exampleQTLs/Pheno_peak25812.txt

geno_peak25812=read.table("../data/exampleQTLs/Geno_11:117048417.txt", col.names = colnames(vcfheadfull), stringsAsFactors = F) %>% select(-CHROM, -POS, -REF, -ALT, -QUAL, -FILTER, -INFO, -FORMAT)
geno_peak25812M=melt(geno_peak25812, id.vars = "ID") %>% separate(value, into=c("genotype", "extra", "extra2"), sep=":") %>% mutate(dose=round(as.integer(extra), digits = 1)) %>% select(-genotype, -extra, -extra2)

geno_peak25812M$variable=as.character(geno_peak25812M$variable)

pheno_peak25812=read.table("../data/exampleQTLs/Pheno_peak25812.txt", col.names = colnames(phenohead), stringsAsFactors = F) 
pheno_peak25812M=melt(pheno_peak25812,id.vars = "chrom") %>% separate(value, into=c("num", "dom"), sep="/") %>% mutate(usage=as.integer(num)/as.integer(dom))

geno_peak25812M$variable=as.character(geno_peak25812M$variable)

Join these

genoPhenopeak25812=pheno_peak25812M %>% inner_join(geno_peak25812M, by="variable") %>% mutate(set=ifelse(variable %in% samplelist_39$line, "old39", "new15"))

Warning: Column `variable` joining factor and character vector, coercing
into character vector

genoPhenopeak25812$set=as.factor(genoPhenopeak25812$set)
genoPhenopeak25812$dose=as.factor(genoPhenopeak25812$dose)

ggplot(genoPhenopeak25812, aes(x=dose, y=usage)) + geom_boxplot()+ geom_jitter(aes(col=set))

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

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

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

other attached packages:
 [1] cowplot_0.9.4   forcats_0.3.0   stringr_1.3.1   dplyr_0.8.0.1  
 [5] purrr_0.3.2     readr_1.3.1     tidyr_0.8.3     tibble_2.1.1   
 [9] ggplot2_3.1.0   tidyverse_1.2.1 gplots_3.0.1    workflowr_1.3.0
[13] gdata_2.18.0    reshape2_1.4.3 

loaded via a namespace (and not attached):
 [1] gtools_3.8.1       tidyselect_0.2.5   haven_1.1.2       
 [4] lattice_0.20-38    colorspace_1.3-2   generics_0.0.2    
 [7] htmltools_0.3.6    yaml_2.2.0         rlang_0.3.1       
[10] pillar_1.3.1       glue_1.3.0         withr_2.1.2       
[13] modelr_0.1.2       readxl_1.1.0       plyr_1.8.4        
[16] munsell_0.5.0      gtable_0.2.0       cellranger_1.1.0  
[19] rvest_0.3.2        caTools_1.17.1.1   evaluate_0.12     
[22] labeling_0.3       knitr_1.20         broom_0.5.1       
[25] Rcpp_1.0.0         KernSmooth_2.23-15 backports_1.1.2   
[28] scales_1.0.0       jsonlite_1.6       fs_1.2.6          
[31] hms_0.4.2          digest_0.6.18      stringi_1.2.4     
[34] grid_3.5.1         rprojroot_1.3-2    cli_1.0.1         
[37] tools_3.5.1        bitops_1.0-6       magrittr_1.5      
[40] lazyeval_0.2.1     crayon_1.3.4       whisker_0.3-2     
[43] pkgconfig_2.0.2    xml2_1.2.0         lubridate_1.7.4   
[46] assertthat_0.2.0   rmarkdown_1.10     httr_1.3.1        
[49] rstudioapi_0.10    R6_2.3.0           nlme_3.1-137      
[52] git2r_0.23.0       compiler_3.5.1