---
title: "FormatMultiQC"
author: "Shelly Trigg"
date: "4/13/2020"
output: html_document
editor_options: 
  chunk_output_type: console
---

load libraries
```{r}
library(data.table)
library(gtools)
library(ggplot2)
library(knitr)
library(kableExtra)
```

read in data 
```{r}
#cutadapt stats
multiqc_a <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200413/multiqc_data/multiqc_cutadapt.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

#bismark
multiqc_b <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200413/multiqc_data/multiqc_bismark_alignment.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

#bismark meth-extract
multiqc_c <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200413/multiqc_data/multiqc_bismark_methextract.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

#bismark dedup
multiqc_d <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200413/multiqc_data/multiqc_bismark_dedup.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

#general stats
multiqc_g <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200413/multiqc_data/multiqc_general_stats.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

#lambda alignment stats
multiqc_l <- fread("https://gannet.fish.washington.edu/metacarcinus/FROGER_meth_compare/20200424/multiqc_data/multiqc_bismark_alignment.txt", header = TRUE, sep = "\t", stringsAsFactors = FALSE)

meta_data <- read.csv("../../metadata/30-244831635 sample submission form.csv", skip = 1, header = TRUE, quote = "", row.names = NULL)



```


format meta data
```{r}
#remove extraneous columns
meta_data <- meta_data[,c("Sample.Name.","Species.Strain.","Library.Prep.Kit.Used.")]

#sub "Pico Methyl-SeqTM Library Prep Kit" wtih WGBS 
meta_data$Library.Prep.Kit.Used. <- gsub("Pico Methyl-SeqTM Library Prep Kit", "WGBS", meta_data$Library.Prep.Kit.Used.)

meta_data$Library.Prep.Kit.Used. <- gsub("EZ DNA Methylation RRBSTM Library Prep Kit Catalog Nos. D5460S", "RRBS", meta_data$Library.Prep.Kit.Used.)

colnames(meta_data) <- c("Sample", "Species", "Sequencing.Method")
```


format multiqc sample names
```{r}
multiqc_a$Sample <- gsub(" \\| ", "_", multiqc_a$Sample)
multiqc_b$Sample <- gsub(" \\| ", "_", multiqc_b$Sample)
multiqc_c$Sample <- gsub(" \\| ", "_", multiqc_c$Sample)
multiqc_d$Sample <- gsub(" \\| ", "_", multiqc_d$Sample)
multiqc_g$Sample <- gsub(" \\| ", "_", multiqc_g$Sample)
```


format column names so merging won't produce duplicate colnames
```{r}
colnames(multiqc_b)[2] <- "uniq_aligned_reads"
colnames(multiqc_d)[2] <- "uniq_aligned_reads_no_discard"

```

format lambda table
```{r}
#make multiqc_l a data.frame
multiqc_l <- data.frame(multiqc_l)

#simplify sample name to just Meth##
multiqc_l$Sample <- gsub("_R1.*","", multiqc_l$Sample)

#create tot_cpg, tot_chh, tot_chg columns

#calculate the total chg counts by summing unmeth chg counts and meth chg counts
multiqc_l$tot_chg <- multiqc_l$meth_chg + multiqc_l$unmeth_chg

#calculate the total chh counts by summing unmeth chh counts and meth chh counts
multiqc_l$tot_chh <- multiqc_l$meth_chh + multiqc_l$unmeth_chh

#calculate the total cpg counts by summing unmeth cpg counts and meth cpg counts
multiqc_l$tot_cpg <- multiqc_l$meth_cpg + multiqc_l$unmeth_cpg

#calculate the unmeht chg + chh counts by summing unmeth chg counts and unmeth chh counts
multiqc_l$unmeth_chg_chh <- multiqc_l$unmeth_chg + multiqc_l$unmeth_chh

#calculate the total chg + chh counts by summing total chg counts and total chh counts
multiqc_l$tot_chg_chh <- multiqc_l$tot_chg + multiqc_l$tot_chh

#calculate percent conversion efficiency by dividing the unmethylated chg + chh counts by the total (meth + unmeth) chg + chh counts
multiqc_l$conversion_efficiency <- (multiqc_l$unmeth_chg_chh/multiqc_l$tot_chg_chh) *100

#remove columns containing strand info (e.g. strand_ot)
multiqc_l <- multiqc_l[,-grep("strand_",colnames(multiqc_l))]

#write out table
write.csv(multiqc_l, "../../Output/lamda_alignments_descriptive_stats.csv", row.names = FALSE, quote = FALSE)

```



#create table with trimming stats
```{r}
#subset the general stats table for the fastqc data
multiqc_f <- multiqc_g[,1:6]
#remove rows with all NAs
multiqc_f <- multiqc_f[complete.cases(multiqc_f),]

#create file ID column in fastqc table
multiqc_f$File.ID <- gsub(".*_Meth","Meth",multiqc_f$Sample)
multiqc_f$File.ID <- gsub("_val.*","",multiqc_f$File.ID)


#create a file ID column in cutadapt table
multiqc_a$File.ID <- gsub(".*_Meth","Meth",multiqc_a$Sample)


#merge cutadapt and fastqc stats by Sample.ID, exclude path info ("Sample" column)
trim_multiqc <- merge(multiqc_a[,-c("Sample")],multiqc_f[,-c("Sample")], by = "File.ID")

#create a sample column to match to meta data
trim_multiqc$Sample <- gsub("_R.*","",trim_multiqc$File.ID)

#merge with meta data
trim_multiqc <- merge(meta_data, trim_multiqc, by = "Sample")

#remove extra text from colnames
colnames(trim_multiqc) <- gsub("FastQC_mqc-generalstats-","", colnames(trim_multiqc))

#order table by Sample.ID 
trim_multiqc <- trim_multiqc[mixedorder(trim_multiqc$File.ID),]

write.csv(trim_multiqc, "../../Output/Raw_trimmed_data_descriptive_stats.csv", row.names = FALSE, quote = FALSE)

```

create a bismark table
```{r}
#merge dedup and meth extract data
multi_bismark <- merge(multiqc_d, multiqc_c, by = "Sample", all = TRUE)

#remove "dedup" from sample column
multi_bismark$Sample <- gsub("_dedup","",multi_bismark$Sample)
multi_bismark$Sample <- gsub("_nodedup","",multi_bismark$Sample)

#merge with bismark alignment report
#exclude meth calls from bismark alignments because better to include these from meth extract report
colnames(multiqc_b[,c(5:7,10:17,19:21)]) #print columns to exclude
multi_bismark <- merge(multiqc_b[,-c(5:7,10:17,19:21)], multi_bismark, by = "Sample", all = TRUE) #merge

#rename Sample column to sample path
colnames(multi_bismark)[1] <- "Sample.path"

#create a sample column to match to meta data
multi_bismark$Sample <- gsub(".*_Meth","Meth", multi_bismark$Sample.path)
multi_bismark$Sample <- gsub("_R.*","", multi_bismark$Sample)


#merge with meta data
multi_bismark <- merge(meta_data,multi_bismark, by ="Sample")

#split out C1 alignments from other alignments

multi_bismark_C1 <- multi_bismark[grepl("C1", multi_bismark$Sample.path ),]
multi_bismark_C1$Sample.path <- "C1"
colnames(multi_bismark_C1)[4] <- "Reference"

#order by sample
multi_bismark_C1 <- multi_bismark_C1[mixedorder(multi_bismark_C1$Sample),]

#write out C1 alignments table
write.csv(multi_bismark_C1,"../../Output/C1_alignments_descriptive_stats.csv", row.names = FALSE, quote = FALSE)


#remove C1 alignment data from object
multi_bismark <- multi_bismark[-grep("C1", multi_bismark$Sample.path ),]

#remove Sample path column
multi_bismark$Sample.path <- NULL

#order data by sample 
multi_bismark <- multi_bismark[mixedorder(as.character(multi_bismark$Sample)),]

#write out alignments table
write.csv(multi_bismark,"../../Output/Pact_Mcap_alignments_descriptive_stats.csv", row.names = FALSE, quote = FALSE)
```

generate simplified table for manuscript main text
```{r}
multi_bismark_simp <- multi_bismark[,c("Sample","Species", "Sequencing.Method","total_reads","no_alignments", "ambig_reads", "uniq_aligned_reads_no_discard", "dedup_reads")]

#add raw reads column
multi_bismark_simp <- merge(data.frame(unique(trim_multiqc[,c("Sample","Species","Sequencing.Method", "r_processed")])), multi_bismark_simp, by = c("Sample","Species","Sequencing.Method"))

#order by sample no
multi_bismark_simp <- multi_bismark_simp[mixedorder(as.character(multi_bismark_simp$Sample)),]


#fill in RRBS cells with uniq_aligned reads
for(i in 1:nrow(multi_bismark_simp)){
  if(is.na(multi_bismark_simp$uniq_aligned_reads_no_discard[i])){
    multi_bismark_simp$uniq_aligned_reads_no_discard[i] <- multi_bismark$uniq_aligned_reads[i]
  }
  if(is.na(multi_bismark_simp$dedup_reads[i])){
    multi_bismark_simp$dedup_reads[i] <- multi_bismark$uniq_aligned_reads[i]
  }
}

#convert to sci not.
#for(i in 4:ncol(multi_bismark_simp)){
#  multi_bismark_simp[,i] <- formatC(multi_bismark_simp[,i], format = "e", digits = 2)
#}

#put table in terms of M reads
for(i in 4:ncol(multi_bismark_simp)){
  multi_bismark_simp[,i] <- round((multi_bismark_simp[,i]/1000000),2)
}

#remove sample column
multi_bismark_simp$Sample <- NULL

#change MBD ... to just MBD in method column
multi_bismark_simp$Sequencing.Method <- gsub("MBD.*", "MBDBS", multi_bismark_simp$Sequencing.Method)

#rename columns
colnames(multi_bismark_simp) <- c("species", "method", "raw", "trimmed", "unaligned", "ambiguously aligned", "uniquely aligned", "uniquely aligned excluding duplicates")

#pretty format
kable(multi_bismark_simp, row.names = FALSE,align='ccccccc') %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))  %>% row_spec(c(1:3,10:12), color = 'black', background = "#BAE4B3") %>% row_spec(c(4:6,13:15), color = 'black', background = "#CBC9E2") %>% row_spec(c(7:9,16:18), color = 'black', background = '#FDBE85') %>% column_spec (1:8,border_left = T, border_right = T) %>% column_spec(6:8, width = "4cm") %>% as_image(file = "../../Output/Summary_descriptive_statistics_Mcap_Pact.jpg")




```

## plot different stats

percent reads with unique alignment 
```{r}
jpeg("../../Output/Pact_Mcap_boxplot_methodXuniquely_mapped_reads.jpg", width = 8, height= 3, units = "in", res = 300 )
ggplot(multi_bismark) + geom_boxplot(aes(x =Sequencing.Method,y = percent_aligned,group = Sequencing.Method, color = Sequencing.Method,fill = Sequencing.Method), alpha = 0.3) + facet_wrap(~Species) + theme_bw() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x=element_blank()) + ylab("uniquely mapped reads (%)")
dev.off()
```

percent CpG methylation 
```{r}
jpeg("../../Output/Pact_Mcap_boxplot_methodXmethCpG.jpg", width = 8, height= 3, units = "in", res = 300 )
ggplot(multi_bismark) + geom_boxplot(aes(x =Sequencing.Method,y = percent_cpg_meth,group = Sequencing.Method, color = Sequencing.Method, fill = Sequencing.Method), alpha = 0.5) + facet_wrap(~Species) + theme_bw() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x=element_blank()) + ylab("CpG methylation (%)")
dev.off()
```

## Other potential plots
percent reads with any alignment (unique or ambig)
```{r}
ggplot(multi_bismark) + geom_boxplot(aes(x =Sequencing.Method,y = ((uniq_aligned_reads + ambig_reads)/total_reads),group = Sequencing.Method, fill = Sequencing.Method)) + facet_wrap(~Species) + theme_bw() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x=element_blank()) + ylab("mapped reads (%)")

```

percent CHG methylation 
```{r}
ggplot(multi_bismark) + geom_boxplot(aes(x =Sequencing.Method,y = percent_chg_meth,group = Sequencing.Method, color = Sequencing.Method, fill = Sequencing.Method), alpha = 0.5) + facet_wrap(~Species) + theme_bw() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x=element_blank()) + ylab("CHG methylation (%)")
```

percent CHH methylation 
```{r}
ggplot(multi_bismark) + geom_boxplot(aes(x =Sequencing.Method,y = percent_chh_meth,group = Sequencing.Method, color = Sequencing.Method, fill = Sequencing.Method), alpha = 0.5) + facet_wrap(~Species) + theme_bw() + theme(axis.title.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x=element_blank()) + ylab("CHH methylation (%)")
```