# Guinea Pig methylation analysis This contains the analysis of RRBS methylation data from 12 guinea pig livers. These were from 6 controls (term animals) and 6 cases (premature animals). * **sample:** RRBS methylation sequence data from liver * **design:** 6 vs 6 case:control (term vs prem) ## `methylkit` processing ### data loading ```{r, eval=FALSE} require(methylkit) # methylkit processing # create list of bismark bam files fileList <- as.list(list.files(path = "./bismark", pattern = "*_bismark_bt2_pe.sorted.bam$", full.names = TRUE)) # create a list of sample IDs, same order as above samples <- as.list(list.files(path = "./bismark", pattern = "*_bismark_bt2_pe.sorted.bam$") %>% strtrim(., 11)) # CpG context # create methylkit files as well as load the data into memory methCalled_CpG <- processBismarkAln(location = fileList, sample.id = samples, treatment = c(rep(0,6), rep(1,6)), assembly="cavPor", save.context= c("CpG"), read.context = "CpG", mincov = 3) ``` ### QC ```{r, eval=FALSE} # filter on coverage; low/min = 5, high/max = 99.9 percentile meth.filtered <- filterByCoverage(methCalled_CpG, lo.count = 5, lo.perc = NULL, hi.count = 100, hi.perc = NULL) # normalise the counts, 'median' meth.filtered <- normalizeCoverage(meth.filtered, method = "median") # check methyltion (beta) distribution getMethylationStats(methCalled_CpG[[2]], plot = TRUE, both.strands = FALSE) getMethylationStats(meth.filtered[[2]], plot = TRUE, both.strands = FALSE) # check coverage getCoverageStats(methCalled_CpG[[2]], plot = TRUE, both.strands = FALSE) getCoverageStats(meth.filtered[[2]], plot = TRUE, both.strands = FALSE) ``` ```{r, eval=FALSE} ## Get bases covered by all samples and cluster samples # merge all samples to one table by using base-pair locations that are covered in all samples # setting destrand=TRUE, will merge reads on both strands of a CpG dinucleotide. This provides better # coverage, but only advised when looking at CpG methylation meth <- methylKit::unite(meth.filtered, destrand = TRUE, mc.cores = 20, min.per.group = 5L) ### Annotation information ```{r, eval=FALSE} require("biomaRt") head(listMarts()) ensembl <- useMart("ensembl") ensembl head(listDatasets(ensembl), n = 25) # ensembl <- useMart("ensembl", dataset="cporcellus_gene_ensembl") ensembl # head(listAttributes(ensembl), n = 25) head(getBM(attributes="chromosome_name", mart=ensembl)) head(getBM(attributes="ensembl_gene_id", mart=ensembl)) listFilters(ensembl) # cavpor mapping information # grab information to map scaffold IDs cavporAnno <- read.delim("http://mirrors.vbi.vt.edu/mirrors/ftp.ncbi.nih.gov/genomes/refseq/vertebrate_mammalian/Cavia_porcellus/representative/GCF_000151735.1_Cavpor3.0/GCF_000151735.1_Cavpor3.0_assembly_report.txt", head = T, skip = 33) # cavporAnnoMap <- cavporAnno[c(5,10)] colnames(cavporAnnoMap)[2] <- 'chr' ``` #### create promoter bed file ```{r, message=FALSE, warning=FALSE} cavporPromoterBed <- data.frame(GenBank.Accn = res$chr, start = res$start - 2000, end = res$start + 200, hgnc_symbol = res$hgnc_symbol) # # map scaffold info cavporPromoterBed <- merge(cavporPromoterBed, cavporAnnoMap, by = "GenBank.Accn") cavporPromoterBed <- cavporPromoterBed[c(5,2:4,1)] cavporPromoterBed <- cavporPromoterBed[!duplicated(cavporPromoterBed),] cavporPromoterBed[cavporPromoterBed$start < 0,]$start <- 1 ``` ### Differential methylation analysis #### regression with overdispersion ```{r, eval=FALSE} myDiff <- calculateDiffMeth(meth, mc.cores = 20, overdispersion = "MN", test = "Chisq") ``` Create a percent methylation matrix: ```{r, eval=FALSE} percent.meth <- percMethylation(meth) percent.meth <- as.data.frame(percent.meth) per.meth <- percent.meth # per.meth <- as.matrix(percent.meth) rownames(per.meth) <- paste0(meth$chr, ':', meth$start) per.meth <- as.data.frame(per.meth) per.meth <- per.meth/100 ``` Generate an annotation GRanges object: ```{r, eval=FALSE, message=FALSE, warning=FALSE} # create a GRanges object methRes <- myDiff class(methRes) <- "data.frame" methRes$end <- methRes$start + 1 methRes <- merge(methRes, cavporAnnoMap, by = 'chr', sort = F) methRes$scaffold <- methRes$chr methRes$chr <- methRes$GenBank.Accn methRes <- methRes[-8] # make GRanges methRes.obj <- makeGRangesFromDataFrame(methRes, keep.extra.columns = TRUE) # get unique scaffolds to find genes diffScaffold <- as.character(unique(methRes.obj@seqnames)) # grab all gene annotations for above scaffolds res <- getBM(attributes=c("chromosome_name", "hgnc_symbol", "entrezgene_id", "strand", "start_position", "end_position", "transcription_start_site", "description"), filters = "chromosome_name", values = diffScaffold, mart = ensembl) head(res, n = 5) colnames(res)[c(1,5,6)] <- c("chr", "start", "end") res$strand <- "*" # geneAnno <- makeGRangesFromDataFrame(res, keep.extra.columns = T) ``` ##### differential methylation analysis ```{r, eval=TRUE} # get differentially methylated bases/regions with specific cutoffs myDiff5 = getMethylDiff(myDiff, difference = 5, qvalue = 0.05, type="all") myDiff5 <- data.frame(myDiff5) myDiff5$absDiff <- abs(myDiff5$meth.diff) myDiff5$ID <- paste(myDiff5$chr, myDiff10$start, sep = ':') ``` Results ```{r, table_caption = "Table: all significant traditional regression sites (abs 5% methylation difference and q<0.05)"} myDiff10[order(myDiff10$pvalue),] ``` Annotate genes into the results: ```{r, eval=TRUE} # create a GRanges object from results methRes.reg <- myDiff5 class(methRes.reg) <- "data.frame" methRes.reg$end <- methRes.reg$start + 1 methRes.reg <- merge(methRes.reg, cavporAnnoMap, by = 'chr', sort = F) methRes.reg$scaffold <- methRes.reg$chr methRes.reg$chr <- methRes.reg$GenBank.Accn methRes.reg <- methRes.reg[-8] head(methRes.reg) # make GRanges methRes.reg.obj <- makeGRangesFromDataFrame(methRes.reg, keep.extra.columns = TRUE) head(methRes.reg.obj) # what regions overlap what genes? overlapGenes.reg <- findOverlaps(methRes.reg.obj, geneAnno) # Return any genes with an overlap. # Convert the resulting "Hits" object to a data frame # and use it as an index overlapGenes.reg.df <- as.data.frame(overlapGenes.reg) head(overlapGenes.reg.df) # got queryHits and subjectHits # extract the regions that hit genes regionsWithHits.reg <- as.data.frame(methRes.reg.obj[overlapGenes.reg.df$queryHits]) # add the names of the genes regionsWithHits.reg$genes <- geneAnno$hgnc_symbol[overlapGenes.reg.df$subjectHits] regionsWithHits.reg$entrezID <- geneAnno$entrezgene_id[overlapGenes.reg.df$subjectHits] regionsWithHits.reg$description <- geneAnno$description[overlapGenes.reg.df$subjectHits] # look for duplication and remove regionsWithHitsNoDup.reg <- regionsWithHits.reg[!duplicated(regionsWithHits.reg), ] regionsWithHitsNoDup.reg <- regionsWithHitsNoDup.reg[order(regionsWithHitsNoDup.reg$pvalue),] regionsWithHitsNoDup.reg$absDiff <- abs(regionsWithHitsNoDup.reg$meth.diff) regionsWithHitsNoDup.reg <- regionsWithHitsNoDup.reg[order(-regionsWithHitsNoDup.reg$absDiff),] # grab TSS information head(res) tssInfo <- res[,c(1,2,7,8)] head(tssInfo) colnames(tssInfo) <- c("seqnames", "genes", "transcription_start_site", "description") # remove non-annotated genes tssInfo <- tssInfo[tssInfo$genes != "",] # merge into results from DM analysis resTSS.reg <- merge(regionsWithHitsNoDup.reg, tssInfo, by = c("seqnames", "genes"), sort = F) # check if TSS is close to start outTSS.reg <- resTSS.reg %>% mutate(ok = start >= (transcription_start_site-2000) & start <= (transcription_start_site + 2000)) table(outTSS.reg$ok) outTSS.reg <- outTSS.reg[outTSS.reg$ok == TRUE,] outTSS.reg <- outTSS.reg %>% mutate(flag = absDiff > 10 & qvalue < 0.05) table(outTSS.reg$flag) finalTSSRes.reg <- outTSS.reg %>% filter(flag == TRUE) ``` Differential methylation for regions DMR analysis metilene ```{r, engine='bash', eval=F} wget http://www.bioinf.uni-leipzig.de/Software/metilene/metilene_v02-8.tar.gz tar -xzvf metilene_v02-8.tar.gz /data/software/metilene_v0.2-8/./metilene_linux64 ``` Need to format the data accordingly (see: https://www.bioinf.uni-leipzig.de/Software/metilene/Start/) ```{r, eval=FALSE} DMRdata <- per.meth # add prefix for metilene analysis colnames(DMRdata) <- paste0(group, "_", colnames(per.meth)) %>% gsub("Sample_", "", .) # create the correct input format DMRdata <- data.frame(chr = gsub(":.*", "", rownames(DMRdata)), pos = gsub(".*:", "", rownames(DMRdata)), DMRdata) # need to sort but faster doing this in bash # output to tab delim file write.table(DMRdata, file = "guinea_cpg_betavals.txt", sep = "\t", row.names = F, quote = FALSE) ``` Sort data in bash: ```{r, engine='bash', eval=F} # sort sort -V -k1,1 -k2,2n guinea_cpg_betavals.txt > guinea_cpg_betavals_sorted.txt # replace NA with - sed -i 's/NA/-/g' guinea_cpg_betavals_sorted.txt # this allows metilene approx the missing values from the beta distribution ``` ##### denovo DMR analysis Run DMR analysis: ```{r, engine='bash', eval=F} /data/software/metilene_v0.2-8/./metilene_linux64 --threads 24 --maxdist 1000 --mincpgs 5 --mtc 2 -a control -b case guinea_cpg_betavals_sorted.txt | sort -V -k1,1 -k2,2n > guinea_metilene_DMRs_sorted.txt ``` > The output for the de-novo DMR annotation mode consists of a bed-like format: > > | chr | start | stop | q-value | mean methylation difference | #CpGs | p (MWU) | p (2D KS) | mean g1 | mean g2 | > > While "mean g1" and "mean g2" refer to the absolute mean methylation level for the corresponding segment in both groups, the difference is given in the 5th column. Single CpGs are not tested using the 2D KS-test. Here, q-values are based on MWU-test p-values. ```{bash, eval=F} head guinea_metilene_DMRs_sorted.txt scaffold_0 339053 339233 0.002386 0.105686 14 3.3789e-07 7.3401e-08 0.23247 0.12678 scaffold_0 1122286 1122390 1 0.127995 5 0.0081351 0.027077 0.8372 0.7092 scaffold_0 1338259 1338293 1 -0.104127 5 0.011467 0.074191 0.48724 0.59136 scaffold_0 1518742 1519088 1 -0.100654 5 0.19836 0.12539 0.632 0.73266 scaffold_0 1759363 1759812 1 0.195237 5 0.0096739 0.034687 0.69949 0.50425 scaffold_0 1788372 1788745 1 0.126111 5 0.0039395 0.011637 0.81717 0.69106 scaffold_0 2258580 2258754 1 0.172655 5 0.0041284 0.026717 0.64867 0.47602 scaffold_0 2512427 2512617 1 0.106062 5 0.16915 0.44325 0.74746 0.6414 scaffold_0 2516037 2516363 1 0.128659 7 0.003065 0.018726 0.70907 0.58041 scaffold_0 2762952 2765018 1 0.119801 12 3.1363e-05 0.0097513 0.79305 0.67325 ``` Can do some basic stats and plotting with their perl script: ```{bash, eval=F} perl /data/software/metilene_v0.2-8/metilene_output.pl -q guinea_metilene_DMRs_sorted.txt -o prelim_results -p 0.05 -c 5 -d 0.1 -a control -b case ``` Load the results into R: ```{r} # load data DMRresults <- read.delim('guinea_metilene_DMRs_sorted.txt', head = F, as.is = F) # label columns colnames(DMRresults) <- c("chr", "start", "end", "qval", "meanDiff", "noCpGs", "pMWU", "p2DKS", "meanControl", "meanCase") # add DMR size DMRresults$DMRsize <- DMRresults$end - DMRresults$start DMRresults <- DMRresults[c(1:3,11,4:10)] head(DMRresults) ``` Annotating ```{r} # create a GRanges object from results methRes.overlap <- DMRresults class(methRes.overlap) <- "data.frame" # methRes.overlap$end <- methRes.overlap$start + 1 methRes.overlap <- merge(methRes.overlap, cavporAnnoMap, by = 'chr', sort = F) methRes.overlap$scaffold <- methRes.overlap$chr methRes.overlap$chr <- methRes.overlap$GenBank.Accn # methRes.overlap <- methRes.overlap[-8] # make GRanges methRes.overlap.obj <- makeGRangesFromDataFrame(methRes.overlap, keep.extra.columns = TRUE) # what regions overlap what genes? overlapGenes.overlap <- findOverlaps(methRes.overlap.obj, geneAnno) # Return any genes with an overlap. # Convert the resulting "Hits" object to a data frame # and use it as an index overlapGenes.overlap.df <- as.data.frame(overlapGenes.overlap) # geneAnno$hgnc_symbol[overlapGenes.reg.df$subjectHits] # extract the regions that hit genes regionsWithHits.overlap <- as.data.frame(methRes.overlap.obj[overlapGenes.overlap.df$queryHits]) # add the names of the genes regionsWithHits.overlap$genes <- geneAnno$hgnc_symbol[overlapGenes.overlap.df$subjectHits] regionsWithHits.overlap$entrezID <- geneAnno$entrezgene_id[overlapGenes.overlap.df$subjectHits] regionsWithHits.overlap$description <- geneAnno$description[overlapGenes.overlap.df$subjectHits] # look for duplication and remove regionsWithHitsNoDup.overlap <- regionsWithHits.overlap[!duplicated(regionsWithHits.overlap), ] # Filter for number of CpGs >=5 and nominal p (pMWU) < 0.05 DMRresultsFilter <- regionsWithHitsNoDup.overlap[order(regionsWithHitsNoDup.overlap$qval),] %>% filter(noCpGs >= 5 & pMWU < 0.05) ``` Adding TSS info: ```{r} tssINFO <- res[,c(1,2,3,7,8)] colnames(tssInfo) <-c("seqnames", "genes", "entrezID", "transcription_start_site", "description") DMRTSS <- merge(DMRresultsFilter, tssINFO, by=c("seqnames", "genes", "entrezID"), sort=F) #define TSS 5000bp upstream and 2000bp downstream outTSS.DMR <- DMRTSS %>% mutate(ok=end >= (transcription_start_site-5000) & end <= (transcription_start_site+2000)) #just DMR in TSS (as defined above) Final_DMR <- outTSS.DMR %>% filter(ok == TRUE) ``` ## Session info ```{r} sessionInfo()