-
Notifications
You must be signed in to change notification settings - Fork 1
/
characterizeTotalApaQtls.Rmd
483 lines (310 loc) · 17.5 KB
/
characterizeTotalApaQtls.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
---
title: "Characterize Total ApaQTLs"
author: "Briana Mittleman"
date: "10/11/2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
This analysis will be used to characterize the total ApaQTLs. I will run the analysis on the total APAqtls in this analysis and will then run a similar analysis on the nuclear APAqtls in another analysis. I would like to study:
* Distance metrics:
+ distance from snp to TSS of gene
+ Distance from snp to peak
* Expression metrics:
+ expression of genes with significant QTLs vs other genes (by rna seq)
+ expression of genes with significant QTLs vs other genes (peak coverage)
* Chrom HMM metrics:
+ look at the chrom HMM interval for the significant QTLs
##Upload Libraries and Data:
Library
```{r}
library(workflowr)
library(reshape2)
library(tidyverse)
library(VennDiagram)
library(data.table)
library(cowplot)
```
Permuted Results from APA:
I will add a column to this dataframe that will tell me if the association is significant at 10% FDR. This will help me plot based on significance later in the analysis. I am also going to seperate the PID into relevant pieces.
```{r}
totalAPA=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"))
totalAPA$sig=as.factor(totalAPA$sig)
print(names(totalAPA))
```
##Distance Metrics
### Distance from snp to TSS
I ran the QTL analysis based on the starting position of the gene.
```{r}
ggplot(totalAPA, aes(x=dist, fill=sig, by=sig)) + geom_density(alpha=.5) + labs(title="Distance from snp to TSS", x="Base Pairs") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + scale_fill_brewer(palette="Paired")
```
It looks like most of the signifcant values are 100,000 bases. This makes sense. I can zoom in on this portion.
```{r}
ggplot(totalAPA, aes(x=dist, fill=sig, by=sig)) + geom_density(alpha=.5)+coord_cartesian(xlim = c(-150000, 150000)) + scale_fill_brewer(palette="Paired")
```
### Distance from snp to peak
To perform this analysis I need to recover the peak positions.
The peak file I used for the QTL analysis is: /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed.bed
```{r}
peaks=read.table("../data/PeaksUsed/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed.bed", col.names = c("chr", "peakStart", "peakEnd", "PeakNum", "PeakScore", "Strand", "Gene")) %>% mutate(peak=paste("peak", PeakNum,sep="")) %>% mutate(PeakCenter=peakStart+ (peakEnd- peakStart))
```
I want to join the peak start to the totalAPA file but the peak column. I will then create a column that is snppos-peakcenter.
```{r}
totalAPA_peakdist= totalAPA %>% inner_join(peaks, by="peak") %>% separate(sid, into=c("snpCHR", "snpLOC"), by=":")
totalAPA_peakdist$snpLOC= as.numeric(totalAPA_peakdist$snpLOC)
totalAPA_peakdist= totalAPA_peakdist %>% mutate(DisttoPeak= snpLOC-PeakCenter)
```
Plot this by significance.
```{r}
ggplot(totalAPA_peakdist, aes(x=DisttoPeak, fill=sig, by=sig)) + geom_density(alpha=.5) + labs(title="Distance from snp peak", x="log10 absolute value Distance to Peak") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + scale_fill_brewer(palette="Paired")
```
Look at the summarys based on significance:
```{r}
totalAPA_peakdist_sig=totalAPA_peakdist %>% filter(sig==1)
totalAPA_peakdist_notsig=totalAPA_peakdist %>% filter(sig==0)
summary(totalAPA_peakdist_sig$DisttoPeak)
summary(totalAPA_peakdist_notsig$DisttoPeak)
```
```{r}
ggplot(totalAPA_peakdist, aes(y=DisttoPeak,x=sig, fill=sig, by=sig)) + geom_boxplot() + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + scale_fill_brewer(palette="Paired")
```
Look like there are some outliers that are really far. I will remove variants greater than 1*10^6th away
```{r}
totalAPA_peakdist_filt=totalAPA_peakdist %>% filter(abs(DisttoPeak) <= 1*(10^6))
ggplot(totalAPA_peakdist_filt, aes(y=DisttoPeak,x=sig, fill=sig, by=sig)) + geom_boxplot() + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + facet_grid(.~strand) + scale_fill_brewer(palette="Paired")
ggplot(totalAPA_peakdist_filt, aes(x=DisttoPeak, fill=sig, by=sig)) + geom_density() + scale_fill_discrete(guide = guide_legend(title = "Significant QTL")) + facet_grid(.~strand)+ scale_fill_brewer(palette="Paired")
```
This gives a similar distribution.
I did snp - peak. This means if the peak is downstream of the snp on the positive strand the number will be negative.
In this case the peak is downstream of the snp.
```{r}
totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="+") %>% filter(DisttoPeak < 0) %>% nrow()
totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="-") %>% filter(DisttoPeak > 0) %>% nrow()
```
Peak is upstream of the snp.
```{r}
totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="+") %>% filter(DisttoPeak > 0) %>% nrow()
totalAPA_peakdist %>% filter(sig==1) %>% filter(strand=="-") %>% filter(DisttoPeak < 0) %>% nrow()
```
This means there is about 50/50 distribution around the peak start.
I am going to plot a violin plot for just the significant ones.
```{r}
ggplot(totalAPA_peakdist_sig, aes(x=DisttoPeak)) + geom_density(fill="darkviolet") + labs(title="Total: Distance from QTL to PAS Peak", x="Distance from SNP to PAS")
```
Within 1000 bases of the peak center.
```{r}
totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 1000) %>% nrow()
totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 10000) %>% nrow()
totalAPA_peakdist_sig %>% filter(abs(DisttoPeak) < 100000) %>% nrow()
```
29 QTLs are within 1000 bp of the peak center, 57 within 10,000bp and 98 within 100,000bp
##Expression metrics
Next I want to pull in the expression values and compare the expression of genes with a total APA qtl in comparison to genes without one. I will also need to pull in the gene names file to add in the gene names from the ensg ID.
Remove the # from the file.
```{r}
expression=read.table("../data/mol_pheno/fastqtl_qqnorm_RNAseq_phase2.fixed.noChr.txt", header = T,stringsAsFactors = F)
expression_mean=apply(expression[,5:73],1,mean,na.rm=TRUE)
expression_var=apply(expression[,5:73],1,var,na.rm=TRUE)
expression$exp.mean= expression_mean
expression$exp.var=expression_var
expression= expression %>% separate(ID, into=c("Gene.stable.ID", "ver"), sep ="[.]")
```
Now I can pull in the names and join the dataframes.
```{r}
geneNames=read.table("../data/ensemble_to_genename.txt", sep="\t", header=T,stringsAsFactors = F)
expression=expression %>% inner_join(geneNames,by="Gene.stable.ID")
expression=expression %>% select(Chr, start, end, Gene.name, exp.mean,exp.var) %>% rename("gene"=Gene.name)
```
Now I can join this with the qtls.
```{r}
totalAPA_wExp=totalAPA %>% inner_join(expression, by="gene")
```
```{r}
ggplot(totalAPA_wExp, aes(x=exp.mean, by=sig, fill=sig)) + geom_density(alpha=.3)+ scale_fill_brewer(palette="Paired")
```
This is not exactly what I want because there are multiple peaks in a gene so some genes are plotted multiple times. I want to group the QTLs by gene and see if there is 1 sig QTL for that gene.
```{r}
gene_wQTL= totalAPA_wExp %>% group_by(gene) %>% summarise(sig_gene=sum(as.numeric(as.character(sig)))) %>% ungroup() %>% inner_join(expression, by="gene") %>% mutate(sigGeneFactor=ifelse(sig_gene>=1, 1,0))
gene_wQTL$sigGeneFactor= as.factor(as.character(gene_wQTL$sigGeneFactor))
```
Therea are 92 genes in this set with a QTL.
```{r}
ggplot(gene_wQTL, aes(x=exp.mean, by=sigGeneFactor, fill=sigGeneFactor)) + geom_density(alpha=.3) +labs(title="Mean in RNA expression by genes with significant QTL", x="Mean in normalized expression") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))+ scale_fill_brewer(palette="Paired")
```
I can do a similar analysis but test the variance in the gene expression.
```{r}
ggplot(gene_wQTL, aes(x=exp.var, by=sigGeneFactor, fill=sigGeneFactor)) + geom_density(alpha=.3) + labs(title="Varriance in RNA expression by genes with significant QTL", x="Variance in normalized expression") + scale_fill_discrete(guide = guide_legend(title = "Significant QTL"))+ scale_fill_brewer(palette="Paired")
```
### Peak coverage for QTLs
I can also look at peak coverage for peaks with QLTs and those without. I will first look at this on peak level then mvoe to gene level. The peak scores come from the coverage in the peaks.
The totalAPA_peakdist data frame has the information I need for this.
```{r}
ggplot(totalAPA_peakdist, aes(x=PeakScore,fill=sig,by=sig)) + geom_density(alpha=.5)+ scale_x_log10() + labs(title="Peak score by significance") + scale_fill_brewer(palette="Paired")
```
This is expected. It makes sense that we have more power to detect qtls in higher expressed peaks. This leads me to believe that filtering out low peaks may add power but will not mitigate the effect.
##Where are the snps:
Download the GM12878 chromHMM annotation. I downleaded this from uscs and put it in:
* /Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.txt
* /project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.txt
Column:
* bin
* chrom
* chromstart
* chromend
* name
* score
* strand
* thich start
* thick end
* item rgb
I can make this a bedfile to use a bedtools pipeline:
* chrom (nochr)
* start
* end
* name (txn hetero ect)
* score
* strand
```{bash,eval=F}
fout = open("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.bed",'w')
for ln in open("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.txt", "r"):
bin, chrom, start, end, name, score, strand, thSt, thE, rgb = ln.split()
chrom=chrom[3:]
name=name.split("_")[0]
fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, score, strand))
fout.close()
```
```{bash,eval=F}
fout = open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.bed",'w')
for ln in open("/Users/bmittleman1/Documents/Gilad_lab/threeprimeseq/data/GM12878.chromHMM.txt", "r"):
bin, chrom, start, end, name, score, strand, thSt, thE, rgb = ln.split()
chrom=chrom[3:]
fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, score, strand))
fout.close()
```
I also need to create a significant QTL snp bedfile for the total qtls. Bed files are 0 bases meaning I want the end to be the position I care about.
chrom, start (pos -1), end (pos), name, score (bh), strand
I can do this in python using the /project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt. I will make the script general to use on the total or nuclear file.
QTLres2SigSNPbed.py
```{bash,eval=F}
def main(inFile, outFile):
fout=open(outFile, "w")
fin=open(inFile, "r")
for num, ln in enumerate(fin):
if num >= 1:
pid, nvar, shape1, shape2, dummy, sid, dist, npval, slope, ppval, bpval, bh = ln.split()
chrom, pos= sid.split(":")
name=sid
start= int(pos)-1
end=int(pos)
strand=pid.split(":")[3].split("_")[1]
bh=float(bh)
if bh <= .1:
fout.write("%s\t%s\t%s\t%s\t%s\t%s\n"%(chrom, start, end, name, bh, strand))
fout.close()
if __name__ == "__main__":
import sys
fraction=sys.argv[1]
inFile = "/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_%s_transcript_permResBH.txt"%(fraction)
outFile= "/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_%s.bed"%(fraction)
main(inFile,outFile)
```
I am going to try to use pybedtools instead for bedtools for this analysis. First I can add it to my conda environment.
Remove header from HMM
```{bash,eval=F}
import pybedtools
sigNuc=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_Nuclear.sort.bed')
sigTot=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/ApaQTLsignificantSnps_10percFDR_Total.sort.bed')
hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")
#map hmm to snps
Tot_overlapHMM=sigTot.map(hmm, c=4)
Nuc_overlapHMM=sigNuc.map(hmm,c=4)
#save results
Tot_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/Tot_overlapHMM.bed")
Nuc_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/perm_APAqtl_trans/sigSnps/Nuc_overlapHMM.bed")
```
I want to make a file that has all of the numbers for the chromatin regions for downstream analysis.
```{bash,eval=F}
cut -f5 GM12878.chromHMM.txt | sort | uniq > chromHMM_regions.txt
```
I then manually seperate the numbers from the name with a tab and remove the name line.
Evaluate results for total:
```{r}
chromHmm=read.table("../data/ChromHmmOverlap/chromHMM_regions.txt", col.names = c("number", "name"), stringsAsFactors = F)
TotalOverlapHMM=read.table("../data/ChromHmmOverlap/Tot_overlapHMM.bed", col.names=c("chrom", "start", "end", "sid", "significance", "strand", "number"))
TotalOverlapHMM_names=TotalOverlapHMM %>% left_join(chromHmm, by="number")
```
```{r}
ggplot(TotalOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Total APAQtls" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))
```
This is the count but I want enrichemnt. I need to randomly choose 188 snps from the ones I tested (nominal res) and get the same inforamtion on where they are.
```{bash,eval=F}
shuf -n 118 /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_NomRes.txt > /project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random118.txt
```
Now I need to make these into the snp bed format.
QTLNOMres2SigSNPbed.py
give this
* total or nuclear
* number of snps
```{bash,eval=F}
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]
number=sys.argv[2]
inFile = "/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_%s_Random%s.txt"%(fraction,number)
outFile= "/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_%s_Random%s.bed"%(fraction,number)
main(inFile,outFile)
```
Sort output
```{bash,eval=F}
import pybedtools
RANDtot=pybedtools.BedTool('/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random118.sort.bed')
hmm=pybedtools.BedTool("/project2/gilad/briana/genome_anotation_data/GM12878.chromHMM.sort.bed")
#map hmm to snps
TotRnad_overlapHMM=RANDtot.map(hmm, c=4)
#save results
TotRnad_overlapHMM.saveas("/project2/gilad/briana/threeprimeseq/data/nominal_APAqtl_trans/randomSnps/ApaQTL_total_Random_overlapHMM.bed")
```
```{r}
TotalRandOverlapHMM=read.table("../data/ChromHmmOverlap/ApaQTL_total_Random_overlapHMM.bed", col.names=c("chrom", "start", "end", "sid", "significance", "strand", "number"))
TotalRandOverlapHMM_names=TotalRandOverlapHMM %>% left_join(chromHmm, by="number")
```
```{r}
ggplot(TotalRandOverlapHMM_names, aes(x=name)) + geom_bar() + labs(title="ChromHMM labels for Total APAQtls (Random)" , y="Number of SNPs", x="Region")+theme(axis.text.x = element_text(angle = 90, hjust = 1))
```
To put this on the same plot I can count the number in each then plot them next to eachother.
```{r}
random_perChromHMM=TotalRandOverlapHMM_names %>% group_by(name) %>% summarise(Random=n())
sig_perChromHMM= TotalOverlapHMM_names %>% group_by(name) %>% summarise(Total_QTLs=n())
perChrommHMM=random_perChromHMM %>% full_join(sig_perChromHMM, by="name", ) %>% replace_na(list(Random=0,Total_QTLs=0))
perChrommHMM_melt=melt(perChrommHMM, id.vars="name")
names(perChrommHMM_melt)=c("Region","Set", "N_Snps" )
```
```{r}
chromenrichTotalplot=ggplot(perChrommHMM_melt, aes(x=Region, y=N_Snps, by=Set, fill=Set)) + geom_bar(position="dodge", stat="identity") +theme(axis.text.x = element_text(angle = 90, hjust = 1)) + labs(title="Enrichment of Total QTLs by chromatin region", y="Number of Snps", x="Chromatin Region") + scale_fill_brewer(palette="Paired")
chromenrichTotalplot
ggsave("../output/plots/ChromHmmEnrich_Total.png", chromenrichTotalplot)
```
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.
```{r}
perChrommHMM$Random= as.integer(perChrommHMM$Random)
perChrommHMM$Total_QTLs= as.integer(perChrommHMM$Total_QTLs)
perChrommHMM_enr=perChrommHMM %>% mutate(Total=Total_QTLs- Random)
```
Write this file so I can put it in the nuclear analysis and compare between groups.
```{r}
write.table(perChrommHMM_enr, file="../data/ChromHmmOverlap/perChrommHMM_Total_enr.txt", quote=F, sep="\t", col.names = T, row.names = F)
```