MeDP-seq count tables at GpG island

... and DESeq2 differential analysis

06_annotationBlocksCount_CGI.R

+library(rtracklayer)
+library(stringr)
+
+# CpG island file, bed to GRanges conversion -----------
+cgi <- import("ucsc_Sscrofa11.2_cpgislands.bed.gz", format = "BED")
+seqlevels(cgi) <- str_remove(seqlevels(cgi), "chr")
+saveRDS(cgi, file = "ucsc_Sscrofa11.2_cpgislands.RDS")
+
+# sarray script -----------
+
+library(dplyr)
+library(purrr)
+library(stringr)
+
+prefix <- "/work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/"
+prefix <- "/home/gdevailly/mnt/genotoul_genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/"
+
+bam_data <- tibble(
+    path = list.files(prefix, pattern = "*.bam$"),
+    id = map_chr(strsplit(path, ".", fixed = TRUE), 1),
+    type = str_sub(id, end = -4)
+)
+
+commands <- paste0(
+    "./annotationBlocksCount.R -i ",
+    prefix,
+    bam_data$path,
+    " -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/",
+    bam_data$id, "_cgi.RDS"
+)
+
+write.table(
+    commands,
+    file = "annotationBlocksCount_CGI_sarray.sh",
+    quote = F, row.names = F, col.names = F
+)
+
+# chmod +x annotationBlocksCount_CGI_sarray.sh
+# module load compiler/gcc-7.2.0
+# module load module load system/R-4.0.2_gcc-7.2.0
+# sarray -J abc -o %j.out -e %j.err -t 01:00:00 --mem=32G --mail-type=FAIL annotationBlocksCount_CGI_sarray.sh
+

07_DEseq2_diff_meth_cgi.R

+library(dplyr)
+library(purrr)
+library(stringr)
+library(DESeq2)
+library(rtracklayer)
+library(ggplot2)
+
+# metadata -----------------
+prefix <- "/work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/"
+# prefix <- "/home/gdevailly/mnt/genotoul_genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/"
+
+bam_data <- tibble(
+    path = list.files(prefix, pattern = "*.bam$"),
+    id = map_chr(strsplit(path, ".", fixed = TRUE), 1),
+    type = str_sub(id, end = -4)
+)
+bam_data$type2 <- str_replace(bam_data$type, "\\+", "_plus")
+bam_data$type2 <- str_replace(bam_data$type2, "\\-", "_minus")
+bam_data$line <- str_extract(bam_data$type2, "^G11_[:alpha:]*")
+bam_data$cond <- str_extract(bam_data$type2, "[:alpha:]*$")
+
+
+my_rds <- paste0("rds_cgi/", bam_data$id, "_cgi.RDS")
+names(my_rds) <- bam_data$id
+
+cgi <- import("ucsc_Sscrofa11.2_cpgislands.bed.gz", format = "BED")
+
+# reading data --------------
+cgi_full <- do.call(cbind, lapply(my_rds, readRDS))
+colnames(cgi_full) <- bam_data$id
+row.names(cgi_full) <- paste0(seq_along(cgi$name), "_", cgi$name)
+cgi$name <- paste0(seq_along(cgi$name), "_", cgi$name)
+
+# DEseq2 -------
+
+my_exp <- grepl("G11", bam_data$type)
+
+dds <- DESeqDataSetFromMatrix(
+    countData = cgi_full[, my_exp],
+    colData = bam_data[my_exp, c("id", "type2", "line", "cond")],
+    design = ~ line + cond + line:cond
+)
+
+keep <- rowSums(counts(dds)) >= 25
+summary(keep)
+# Mode   FALSE    TRUE 
+# logical    2055   44164 
+dds <- dds[keep,]
+
+dds <- DESeq(dds)
+
+rline <- results(dds, c("line", "G11_plus", "G11_minus"))
+rcond <- results(dds, c("cond", "fed", "fasted"))
+rinte <- results(dds, list("lineG11_plus.condfed", character()))
+
+plotMA(rline)
+plotMA(rcond)
+plotMA(rinte)
+
+# results exploration under igv ----------------
+
+plotCounts(dds, gene=which.min(rline$padj), intgroup="line")
+
+signif <- data.frame(rline)[order(rline$padj, na.last = NA), ]
+signif <- signif[signif$padj <= 0.01, ]
+
+cgi[cgi$name == "10327_CpG:_75", ]
+names(cgi) <- cgi$name
+mr <- cgi[rownames(signif), ]
+
+mr$baseMean <- signif$baseMean
+mr$log2FoldChange <- signif$log2FoldChange
+mr$padj <- signif$padj
+
+options("showHeadLines"= 25)
+
+mr

annotationBlocksCount_CGI_sarray.sh

+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R4_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R5.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R5_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fasted_R6.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fasted_R6_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fed_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fed_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fed_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fed_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fed_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fed_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fed_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fed_R4_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11-_fed_R5.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11-_fed_R5_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R4_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R5.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R5_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fasted_R6.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fasted_R6_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R4_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R5.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R5_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/G11+_fed_R6.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/G11+_fed_R6_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-input_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-input_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-input_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-input_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-input_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-input_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-input_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-input_R4_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-methyl_R1.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-methyl_R1_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-methyl_R2.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-methyl_R2_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-methyl_R3.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-methyl_R3_cgi.RDS
+./annotationBlocksCount.R -i /work2/genphyse/genepi/data/porc/medpseq/ROSEpigs_processed2/results/bwa/mergedLibrary/Tem-methyl_R4.mLb.clN.sorted.bam -w /work2/genphyse/genepi/guillaume/rosepigs/ucsc_Sscrofa11.2_cpgislands.RDS -o rds_cgi/Tem-methyl_R4_cgi.RDS