Dinh/Dinh 2017/NOTES/2017-5-29
Jump to navigation
Jump to search
MONOD reanalysis - part 2[edit]
Data sets[edit]
- Finished all WGBS haploinfo files with rmdup
- Sticking with getHaplo script for now
- Planning to add more Encode RRBS
- Working on reorganizing all the data & making backups.
- generating additional MHL matrices (for the same samples)
- 1) normal MHB, 2) filtered with at least 3 cpgs, 3) merged MHBs from 100 bp apart, 4) new MHBs?
Analysis plan[edit]
- MHL is valuable
- Identify tumor specific blocks using MHL and perform quantitative estimation of tumor load in cell-free DNA
- Identify tissue specific blocks using MHL and perform tissue of origin mapping and deconvolution in cell-free DNA
1. MHL is valuable[edit]
- We want to show that MHL have better signal to noise than AMF
- Fair comparison with AMF means using a published dataset with published differentially methylated regions (DMRs) that we can calculate both MHL and AMF for. Then we will use GSI to compute the top 100 regions using either MHL or AMF for each tissue group.
- So using 36 tissues sample from Schultz et. al. 2015 with DMRs from their website. http://neomorph.salk.edu/SDEC_tissue_methylomes/processed_data/code_data.tar.gz
- I wasn't able to find the tissue specificity labels for the Schultz et. al. DMRs, but using GSI to identify the top 100 tissue specific regions should be fair to both metrics.
- 626,418 total DMR regions in DMR_methylpy_matrix.bed
- merge haploinfo files to DMR regions, example shell command below:
~/WGBS_HaploInfo/scripts/mergeHaploInfo_May22_bed.pl DMR_methylpy_matrix.bed < STL001BL-01.chr22.sorted.clipped.bam.pe.hapInfo.txt >> STL001BL-01.allchrs-mpy.sorted.clipped.bam.pe.dmrs
- from merged haploinfo files, use hapinfo2mhl.pl and hapinfo2amf.pl to calculate data matrices. I computed one chromosome at a time, then use mergeMatrices.pl to merge them together.
- Final data matrices:
amf.dmrs.txt mhl.dmrs.txt
- Analysis code - File:Monod reanalysis 1 functions.R
- Perform analysis using R:
source("1_functions.R") # arguments for running in Rscript mode args <- commandArgs(trailingOnly=T) mhlfile<-args[1] amffile<-args[2] # read in the data matrices mhldata<-read.table(mhlfile, T, row=1) amfdata<-read.table(amffile, T, row=1) # sample 100,000 regions keep.rows<-sample(nrow(amfdata), 100000, replace=F) # compare regions selected by AMF compare.matrices.by.top.gsi(amfdata[keep.rows,], "AMF-dmrs", mhldata, "MHL-dmrs") # Outputs: #[1] "Cur diag: 0.723496010204713" #[1] "Cur off diag: 0.394425152972905" #[1] "Alt diag: 0.541506924678588" #[1] "Alt off diag: 0.256813604981447" # compare regions selected by MHL compare.matrices.by.top.gsi(mhldata[keep.rows,], "MHL-dmrs", amfdata, "AMF-dmrs") # Outputs: #[1] "Cur diag: 0.730646020873588" #[1] "Cur off diag: 0.296063402306549" #[1] "Alt diag: 0.738326381821633" #[1] "Alt off diag: 0.505675275381217"
- Results:
500px