-
Notifications
You must be signed in to change notification settings - Fork 2
/
mapapaQTL.Rmd
90 lines (54 loc) · 2.11 KB
/
mapapaQTL.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
---
title: "apaQTLs"
author: "Briana Mittleman"
date: "4/18/2019"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
In this analysis I will call apaQTls in both fractions. I will start with the phenotype files and normalized the counts using the leafcutter package in order to run the fastq QTL mapper.
##Prepare phenotypes for QTL- phenotype dir
It is best to run this analysis in the data/phenotype directory. I have copied the leafcutter prepare_phenotype_table.py to the code directroy to use here.
```{bash,eval=F}
#!/bin/bash
module load python
gzip APApeak_Phenotype_GeneLocAnno.Total.fc
gzip APApeak_Phenotype_GeneLocAnno.Nuclear.fc
python ../../code/prepare_phenotype_table.py/APApeak_Phenotype_GeneLocAnno.Total.fc
python ../../code/prepare_phenotype_table.py/APApeak_Phenotype_GeneLocAnno.Nuclear.fc
```
This will output bash scripts to run.
```{bash,eval=F}
module load Anaconda3
source activate three-prime-env
sh APApeak_Phenotype_GeneLocAnno.Nuclear.fc.gz_prepare.sh
sh APApeak_Phenotype_GeneLocAnno.Nuclear.fc.gz_prepare.sh
```
Subset the PCs to use the first 2 in the qtl calling:
```{bash,eval=F}
module load Anaconda3
source activate three-prime-env
head -n 3 APApeak_Phenotype_GeneLocAnno.Nuclear.fc.gz.PCs > APApeak_Phenotype_GeneLocAnno.Nuclear.fc.gz.2PCs
head -n 3 APApeak_Phenotype_GeneLocAnno.Total.fc.gz.PCs > APApeak_Phenotype_GeneLocAnno.Total.fc.gz.2PCs
```
##Call QTLs- code dir
Next I will need to make a sample list. From the code directory:
```{bash,eval=F}
python makeSampleList.py
```
Prepare directroy
```{bash,eval=F}
mkdir ../data/apaQTLNominal
mkdir ../data/apaQTLPermuted
```
Run the code to call QTLs within 1mb of each PAS peak. I run both a nominal pass and a permuted pas. The permulted pas chosses the best snp for each peak gene pair.
```{bash,eval=F}
sbatch apaQTL_Nominal.sh
sbatch apaQTL_permuted.sh
```
Concatinate all of the results in the permuted set. I do this so I can account for multiple testing with the benjamini hochberg test.
**Concatinate**
```{bash,eval=F}
Rscripts apaQTLCorrectPvalMakeQQ.R
```