---
title: "mqtl_analysis"
author: "Katherine Silliman"
date: "2/21/2022"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## R Markdown

```{r}
library(MatrixEQTL)
library(tidyverse)
library(reshape2)
library(LEA)
```

```{r}
pos <- read.table("../data/Cv_10x_fxmerge_05.012.pos",header = F,sep = "\t",stringsAsFactors = F)
head(pos)
row.names(pos) <- paste(pos$V1,pos$V2,sep = "_")
id <- as.character(row.names(pos))

pos <- cbind(id,pos)
colnames(pos) <- c('snp','chr','pos')
write.table(pos,file="../data/Cv_10x_fxmerge_05.pos",quote = F,sep = "\t",)

inds <- scan("../data/Cv_10x_fxmerge_05.012.indv", character())
inds <- str_remove(inds, "/home/sr320/github/2018_L18-adult-methylation/bg_data/")
inds <- str_remove(inds,"_R1_val_1_bismark_bt2_pe.deduplicated.sorted.bam")

geno <- t(read.table("../data/Cv_10x_fxmerge_05.012",header = F,sep = "\t",stringsAsFactors = F,row.names = 1))

geno <- cbind(id,geno)
geno <- as.data.frame(geno,stringsAsFactors = F)
rownames(geno) <- NULL
# change sample names to match epigenetics
colnames(geno) <- c("id",inds)
geno[geno < 0] <- NA
head(geno)
write.table(geno,"../data/Cv_10x_fxmerge_05.geno",sep="\t",quote = F,row.names = F)
```

	chrom	start	end	12	13	16	19	22	23	29	31	35	36	39	3	41	44	48	50	52	53	54	59	64	6	76	77	7	9	total

```{r}
mf <- read.table("../data/union_10x.bedgraph",header=T,stringsAsFactors = F)
#add sample names

head(mf)
dim(mf)

```

Remove CpGs where the difference in methylation between the samples with the highest and lowest values is less than 10%.
```{r}
mf1 <- mf %>% mutate_at(vars(contains("X")), as.numeric) %>%
  rowwise() %>%
  mutate(max=max(c_across(contains("X"))), min = min(c_across((contains("X"))))) 
  
  
  
  #rowwise() %>% mutate(max = max(as.diff = pmax(contains("X"))-pmin(contains("X")))
```

```{r}
mf1 <- mf1 %>% mutate(diff = max-min, loci = paste(chrom,start,end,sep="_"))
```

Plot range per locus before filtering
```{r}
ggplot(mf1,aes(x=diff))+ geom_histogram(binwidth = 1)
```


```{r}


mf_filt <- mf1 %>%
  filter(diff>10) %>%
  select(loci, contains("X"),-max) %>% 
  rename_with(~str_remove(.,"X")) %>%
  rowwise() %>%
   mutate(count_na = sum(is.na(cur_data()))) %>%
   ungroup %>% filter(count_na/26 <0.1)
dim(mf_filt)
```

```{r}
save(mf_filt,file="../data/mf_filt.Robj")
```


Plot range per locus after filtering

```{r}
mf_filtdiff <- mf_filt %>% 
  mutate_at(vars(contains("M"),contains("F")), as.numeric) %>%
  rowwise() %>%
  mutate(max=max(c_across(contains(c("M","F")))), min = min(c_across(contains(c("M","F"))))) %>%
  mutate(diff = max-min)
```

```{r}

ggplot(mf_filtdiff,aes(x=diff))+ geom_histogram(binwidth = 1)
```

```{r}
perc_meth <- melt(mf_filt,variable.name = "sample",value.name = "percMeth",id.vars = "loci") 
head(perc_meth)
```

```{r}
# check a couple loci to make sure % methylation is calculated separately for each sample and gene region 
perc_meth %>% filter(loci=="NC_035780.1_1882_1884") %>% ggplot(aes(x=sample, y=percMeth)) + geom_bar(stat="identity")
```
## Rank-based normalization
```{r}
myRN <- function(x){
    
    y <- x[2:27][!is.na(x[2:27])]
    n <- length(y)
    k = 0.375  
  # Ranks.
      r <- rank(y)
  
  # Apply transformation.
    out <- qnorm((r - k) / (n - 2 * k + 1))
    x[2:27][!is.na(x[2:27])] <- out
    return(x)
}
```

```{r}
mf_filt.df <- as.data.frame(select(mf_filt,-count_na))
df_list <- lapply(1:nrow(mf_filt.df), function(r_num) { myRN(mf_filt.df[r_num, ]) })
dcN10<- do.call(rbind, df_list)
```

```{r}
write.table(dcN10,"../data/meth_filter_10_RankNorm.tab",sep="\t",row.names = F,quote = F)
head(dcN10)
```

```{r}
# get meth position file
dcN10 <- as.data.frame(dcN10,stringsAsFactors = F)
chr1 = sapply(strsplit(dcN10$loci,"_"), '[',1)
chr2 = sapply(strsplit(dcN10$loci,"_"), '[',2)
chr = paste(chr1,chr2,sep="_")
s1 = as.numeric(sapply(strsplit(dcN10$loci,"_"), '[',3))
s2 = as.numeric(sapply(strsplit(dcN10$loci,"_"), '[',4))
methpos <- as.data.frame(cbind(dcN10$loci,chr,s1,s2),stringsAsFactors = F)
methpos$s1 <- as.numeric(methpos$s1)
methpos$s2 <- as.numeric(methpos$s2)
colnames(methpos) <- c('geneid','chr','s1','s2')
head(methpos)
write.table(methpos,file="../data/methsitespos_10_mqtl.tab",quote = F,row.names = F,sep = "\t")
```

Look at Lauras code to get cov
```{r}
read.csv("../RAnalysis/data/adult-meta.csv", header=T, na.strings = NA) %>%
  select(Sample.ID, Treatment, Sex) %>%
  mutate_at(vars(Treatment, Sex), as.factor) %>%
  mutate_at(vars(Treatment, Sex), as.numeric) %>%
  column_to_rownames("Sample.ID") %>% t() %>% as.data.frame() %>% rownames_to_column(var = "id") %>%
  write.table("../data/covariates.txt", sep="\t",col.names = TRUE, row.names = FALSE, quote = FALSE) # save to file

```

# Run methQTL
```{r}
# set parameters
useModel = modelLINEAR
snpfile = "../data/Cv_10x_fxmerge_05.geno"
methfile = "../data/meth_filter_10_RankNorm.tab"
covfile = "../data/covariate_mqtl.tab"

snppos = read.table("../data/Cv_10x_fxmerge_05.pos",header=T)
methpos = read.table("../data/methsitespos_10_mqtl.tab",header=T)

pvOutputThresholdcis = 0.05
pvOutputThresholdtrans = 0.01
output_file_name_cis = "../output/mQTL_sites_cis50kb_12_p05_maf05.txt"
output_file_name_tra = "../output/mQTL_sites_trans50kb_12_p01_maf05.txt"
#errorCovariance = rab2;
errorCovariance = numeric()
# Distance for local gene-SNP pairs
cisDist = 50000  #50kb
```

```{r}
snps = SlicedData$new();
snps$fileDelimiter = "\t";      # the TAB character
snps$fileOmitCharacters = "NA"; # denote missing values;
snps$fileSkipRows = 1;          # one row of column labels
snps$fileSkipColumns = 1;       # one column of row labels
snps$fileSliceSize = 2000;      # read file in pieces of 2,000 rows
snps$LoadFile( snpfile);
```

```{r}
meth = SlicedData$new();
meth$fileDelimiter = "\t";      # the TAB character
meth$fileOmitCharacters = "NA"; # denote missing values;
meth$fileSkipRows = 1;          # one row of column labels
meth$fileSkipColumns = 1;       # one column of row labels
meth$fileSliceSize = 5000;      # read file in pieces of 2,000 rows
meth$LoadFile( methfile )
```

```{r}
cvrt = SlicedData$new();
cvrt$fileDelimiter = "\t";      # the TAB character
cvrt$fileOmitCharacters = "NA"; # denote missing values;
cvrt$fileSkipRows = 1;          # one row of column labels
cvrt$fileSkipColumns = 1;       # one column of row labels
cvrt$fileSliceSize = 2000;      # read file in pieces of 2,000 rows
cvrt$LoadFile( covfile )
```

```{r}
me = Matrix_eQTL_engine(
    snps = snps,
    gene = gene,
    cvrt = cvrt,
    output_file_name = output_file_name_tra,
    output_file_name.cis = output_file_name_cis,
    pvOutputThreshold = pvOutputThresholdtrans,
    pvOutputThreshold.cis = pvOutputThresholdcis,
    useModel = useModel, 
    errorCovariance = errorCovariance, 
    snppos = snppos,
    genepos = genepos,
    cisDist = cisDist,
    verbose = TRUE,
    pvalue.hist = TRUE,
    min.pv.by.genesnp = FALSE,
    noFDRsaveMemory = FALSE);
```


|Sample.ID|OldSample.ID|Treatment|Sex|TreatmentN|Parent.ID|
|---------|------------|---------|---|----------|---------|
|12M      |S12M        |Exposed  |M  |3         |EM05     |
|13M      |S13M        |Control  |M  |1         |CM04     |
|16F      |S16F        |Control  |F  |2         |CF05     |
|19F      |S19F        |Control  |F  |2         |CF08     |
|22F      |S22F        |Exposed  |F  |4         |EF02     |
|23M      |S23M        |Exposed  |M  |3         |EM04     |
|29F      |S29F        |Exposed  |F  |4         |EF07     |
|31M      |S31M        |Exposed  |M  |3         |EM06     |
|35F      |S35F        |Exposed  |F  |4         |EF08     |
|36F      |S36F        |Exposed  |F  |4         |EF05     |
|39F      |S39F        |Control  |F  |2         |CF06     |
|3F       |S3F         |Exposed  |F  |4         |EF06     |
|41F      |S41F        |Exposed  |F  |4         |EF03     |
|44F      |S44F        |Control  |F  |2         |CF03     |
|48M      |S48M        |Exposed  |M  |3         |EM03     |
|50F      |S50F        |Exposed  |F  |4         |EF01     |
|52F      |S52F        |Control  |F  |2         |CF07     |
|53F      |S53F        |Control  |F  |2         |CF02     |
|54F      |S54F        |Control  |F  |2         |CF01     |
|59M      |S59M        |Exposed  |M  |3         |EM01     |
|64M      |S64M        |Control  |M  |1         |CM05     |
|6M       |S6M         |Control  |M  |1         |CM02     |
|76F      |S76F        |Control  |F  |2         |CF04     |
|77F      |S77F        |Exposed  |F  |4         |EF04     |
|7M       |S7M         |Control  |M  |1         |CM01     |
|9M       |S9M         |Exposed  |M  |3         |EM02     |

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see <http://rmarkdown.rstudio.com>.

When you click the **Knit** button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

```{r cars}
summary(cars)
```

## Including Plots

You can also embed plots, for example:

```{r pressure, echo=FALSE}
plot(pressure)
```

Note that the `echo = FALSE` parameter was added to the code chunk to prevent printing of the R code that generated the plot.
