Ns126:MONOD

From ZhangLabWiki
Revision as of 19:26, 25 March 2016 by >Shicheng (→‎Computational Script)
Jump to navigation Jump to search

Wiki Index

Results

Library and Public dataset

  • library = non-directional or directional?
  • Single-end or Pired-end
  • Phred score = 33 or 64?
  • Read Length >= 75 or not?
  • WGBS or RRBS

Computational Script

  • bam to haploinfo
* bam2hapinfo-V1.0.pl
* Debug:Update read strain determine function, March 2016: $strand<16 to $strand & 16 or 0x10
* update again: Now, bam2hapinfo.pl works for BAMs created by bismark and bisReadMapper. 
  • haploinfo to MHL
* hapinfo2mhl-V2016.pl

RD10 and 80bp up Region Definition

  • Define WGBS frequently coverage genomic regions
  • Read Dept >= 10
  • Length >= 80bp
Merge all the Bam files: MergeAllBamFilesWithBash.sh[1]
/home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools genomecov -bg -split -ibam N37_10_tissue_pool_chr1.merged.bam>   N37_10_tissue_pooled.chr1.genomecov.bed
awk '$4>9 { print $1"\t"$2"\t"$3}'  N37_10_tissue_pooled.chr1.genomecov.bed | /home/kunzhang/softwares/bedtools-2.17.0/bin/bedtools merge -d 10 -i - > N37_10_tissue_pooled.chr1.RD10.genomecov.bed
awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' N37_10_tissue_pooled.chr1.RD10.genomecov.bed > N37_10_tissue_pooled.chr1.RD10_80up.genomecov.bed

GWBS alignment with Bismark

  • Alignment with Bismark in glean queue
*prepare the sample information manually, therefore, you can alignment single-end and paired-end BS-seq Fastq at same time. Sam_Config_File Format
perl ~/bin/bismark2bamPBS.pl Sam_Config_File submit
*check which alignment is failed and therefore re-submit the job again:
perl ~/bin/bismarkAlignSucceedCheck.pl
*check alignment ratio/coverage:
perl ~/bin/bismarkMappingStatistic.pl

Bam to Hapinfo (Update: March/16/2016)

  • bam and hapinfo files
* N37: /home/shg047/oasis/N37/bam
* SALK: 

  • Blood:
** ln -s /home/k4zhang/my_oasis_tscc/MONOD/whole_blood_WGBS/BAMfiles  bam
** cd /home/shg047/oasis/AgeBlood
** perl ~/bin/mergeBamByChrosome.pl /home/shg047/oasis/AgeBlood/bam /home/shg047/oasis/AgeBlood/mergeBam submit
  • HCC:
  • Plasma:
**Phase  I:
*** where is fastq?
***/oasis/tscc/scratch/ddiep/Working/Rerun_rrbs/BAMfiles 
***/oasis/tscc/scratch/ddiep/Working/Rerun_rrbs/MethylFreq
**Phase II: 
*** where is fastq?
*** /oasis/tscc/scratch/ddiep/Plasma_RRBS_151208/BAMfiles
**Phase III:
***
  • Estellar: /home/shg047/oasis/Estellar2016

Hapinfo to Methylation LD Heatmap(Update: March/16/2016)

Hapinfo to MHB (Update: March/16/2016)

Hapinfo to MHL matrix

  • WB
 cd /home/shg047/oasis/AgeBlood/hapinfo
 echo "perl ~/bin/hapinfo2mhl.pl ./ > ../wb.mhl.march22.txt" | qsub -N wbc.hapinfo2mhl -q pdafm
 WB_hapinfo_mf.pbs
 WB_hapinfo_mhl.pbs
  • Simulation: Figure 2 (manuscript-code and data).
cd /home/shg047/oasis/monod/test/hapinfo2mhl
perl ~/bin/hapinfo2mhl.pl ./
cd /home/shg047/oasis/monod/test/hapinfo2mhb
cat ../hapinfo2mhl/*hapInfo.txt > MergeInfo.txt
perl ~/oasis/monod/bin/hapinfo2R2ByBed.pl rlt chr10:100027918-100027944  < MergeInfo.txt
*File:36CF.tm.png File:6129.tm.pngFile:757E.tm.png
  • Real Data
cd /oasis/tscc/scratch/shg047/AgeBlood/mhl
* Relationship between MF and MHL in WB(Age) data: Rcode.R(Attention: here 1-9 haplotype were used. In Dr. Zhang's idea, only 3-7 haplotype were used in version 1 script)
* smoothscatter.plot.monod.wb.R
* File:Monod.3FD8.tm.pngFile:6229.tm.png


miss ratio: 43.63%: sum(is.na(data))/(nrow(data)*ncol(data))

Estellar2016

* Estellar: cd /home/shg047/oasis/Estellar2016/mergeHapinfo
            Rscript --vanilla ~/bin/matrixPlot.R -i "Estellar2016.mhl.march19.txt" -o "Esterllar2016"

Strategy for genome-wide DNA methylation sequencing for cancer and plasma samples

  • Linkage disequilibrium (R2) of methylation status between pair-wise CpG loci within the methylation blocks

DNA methylation linkage disequilibrium inferred by the methylation status from the single read of the BS-seq were conducted in the present study. The R2 of pair-wise CpG loci within the methylation block were calculated


  • Number of MHBs with R2
scp shg047@genome-miner:/home/kunzhang/CpgMIP/MONOD/Data/WGBS_data/mld_block_hapInfo_July2015/All_chromosomes_combined/* /home/shg047/oasis/monod/mhb/hapinfo

Methylation haplotype blocks were conservative in human tissues


  • how to test the significance between PBMC and TCGA methylation haplotype blocks File:R20150818.txt

Cancer Diagnosis Performance based on Methylation Haplotype Loading of circulating free DNA methylation

In this section, RRBS and Capseq dataset were taken as the discovery dataset and the BSPP dataset were used validate the performance of the biomarkers. For RRBS dataset, 68 cancer and 8 normal samples were included. 56046 regions were obtained by methylation haplotype loading algorithm (MHL) and 17263 regions were 100% un-methylated in 8 normal samples. 248 regions among 56046 were methylated in at least 50% cancer samples while 100% un-methylated in normal samples. The average length of the regions were 103bp (IQR=95bp, SD=126) and were located in the promoter region of 182 genes (2000bp up-stream of TSS). With the help of text mining, 21 genes of them were validated to be methylation relevant cancer related genes (Table). In order to evaluate the distinguish ability of DNA methylation from cancer to the normal samples, Random forest model were applied. In the procedure of RF model, the optimal number of variables tried at each split is131 (mtry) and number of trees (mtrees) makes no difference to the prediction accuracy from 250 to 7000. RF prediction model showed 1592 regions could provide positive ability to distinguish cancer samples from normal samples with sensitivity of 97.06%, specificity of 100% and accuracy of 97.37%. Furthermore, we found the top 206 most importance regions could take account of 80.0% contribution to the accurate prediction. And the random forest model based on top 206 regions with mtry of 14 and mtrees of 500 showed better performance with 100% sensitivity, 100% specificity and 100% accuracy. The result showed these 206 regions were high frequent hypermethylated in cancers (Mean=41.4%, SD=14.5%, IQR=22.4%). (Supplementary Table **) For CapSeq dataset, 31 cancer 24 normal samples were included. 86206 regions were obtained by methylation haplotype loading algorithm (MHL), however, only 2205 regions were 100% un-methylated in 24 normal samples. Parameters of Random forest model were tuned and the optimal number of variables tried at each split is 92 (mtry) and number of trees (mtree) makes no difference to the prediction accuracy from 250 to 7000. Random forest algorithm showed 516 regions could provide positive ability to separate cancer samples from normal samples with sensitivity of 87.1%, specificity of 75.0% and accuracy of 81.26%. We found that the top 21 most importance regions could take account of 80.0% contribution for the classification. Therefore, the second round prediction process of random forest model based on top 21 regions with mtry of 7 and mtrees of 750 showed perfect classification with 96.77% sensitivity, 100% specificity and 98.18% accuracy. The average hyper-methylation frequency in cancer samples for above 516 and 21 regions 90.43% (SD=0.08, IQR=0.096) and 96.7% (SD=0.041, IQR=0.056), respectively. (Supplementary Table ** and supplementary Table **) For BSPP dataset, 16 cancer plasma and 16 normal plasma samples were included. 36281 regions were obtained by methylation haplotype loading algorithm (MHL) and 4194 regions were 100% un-methylated in 16 normal samples. All the regions hypermethylated in cancer and 100% un-methylated in normal samples were collected from RRBS, SeqCap and BSPP dataset and 30 regions were found un-methylation in all the normal samples of the 3 dataset. These 30 regions could explain 93.75% cancers incidence and the specificity is 100%. The methylation frequency of these 30 regions ranged from 6.25% to 25% in free-cell circulating DNA while it is hypermethylated in 1.47%-44% RRBS cancer samples and 6.45%-77% SeqCap cancer samples. These 30 regions were located in the promoter or gene body of 24 genes, including FGF19, MEF2C-AS1, MEF2C, HABP2, AGO3, TAF15, ZWILCH, YTHDF3-AS1, ZNF213, C9orf114, C14orf28, ADAM11, CDKN2AIP, TMEM87B, CDKL3, MSH5, MSH5-SAPCD1, UBC, PCM1, PCK2, NR2F6, JADE2, KANK1 and MCM3. At least 12 genes have been demonstrated to be associated with human cancers with the annotation from NCBI (HABP2, FGF19, MEF2C, AGO3, CDKL3, MCM3, TAF15, PCK2, NR2F6, MSH5, ZWILCH and PCM1). For example, the FGF19-FGFR4 signaling axis has been implicated in the pathogenesis of several cancers in mice and potentially in humans (7) while HABP2 has been found to be hypermethylated in hepatocellular carcinomas cancer (8). AGO3 was found significant decreased in human liver cancer (9). KANK1 re-expression induced by 5-Aza-2'-deoxycytidine could suppresses nasopharyngeal carcinoma cell proliferation and promotes apoptosis (10). Whole-exome sequencing identifies mutated PCK2 associated with carcinoma cell proliferation in a hepatocellular carcinoma patient(11)

Genome-wide DNA methylation profile comparison between solid tissue and circulating free DNA methylation in plasma

All the assumption of we can used DNA methylation in free-cell circulating DNA methylation is based on the methylation status of solid cancers can be released into plasma without any specific selection, thus, the aberrant DNA methylation fragment or biomarkers identified in solid tissues can be used to be the target for cancer non-invasive diagnosis, screening or prognosis surveillance. We then constructed 15 genome-wide DNA methylation profiles for the paired plasmas and solid tissues from same cancer patient with RRBS assay to check whether DNA methylation fragments derived from cancer cells were released into the blood randomly or selectively. 5 paired lung cancer, colon cancer and pancreatic cancer solid tissues and plasma were enrolled in this section and single-base methylome were established by RRBS as mentioned above. Correlation analysis based on 56046 MHL regions showed there is just weak correlation between genome-wide DNA methylation of solid tissues and plasmas. The average correlation between tissue and plasma was 0.50 (95%CI: 0.44-0.56)(Figure). What’s more, only about 28% (95%CI: 25%-31%) regions were simultaneously methylated both in tissues and plasma. Even considering different threshold to define the methylation status and transferring the continuous MHL to binary DNA methylation status, the maximum co-methylation ratio was still less than 63% (Figure). These evidence showed the dramatically significant heterogeneity between solid and plasma DNA methylation profile and provide the importance that we should validate all the biomarkers identified in solid tissues in cell-free circulating plasma DNA to guarantee the validity of the biomarkers. 25204 regions which was methylated in sample pairs for at least one time were selected to evaluate the selectivity of hyper-methylated DNA fragments from solid tissues to plasma. The results showed the process of the releasing were non-random, with at least 1590 fragment were significantly prefer-selected (P<1.98*10-6, binomial test). These fragments were located in the regions of 1190 genes. Function enrichment analysis showed these genes were significantly associated with cancer relevant biological functions (Table), including embryonic morphogenesis, regulation of transcription, neuron differentiation, regionalization, tissue morphogenesis ,transcription factor activity, sequence-specific DNA binding ,transcription regulator activity.

Re-validation of the hyper-methylated fragments of plasma in public GEO dataset.

GSE56044 (124 lung cancer and 12 normal), GSE39279 (444 lung cancer), GSE52401 (244 normal lung), TCGA-lung cancer (), TCGA-colon cancer () and TCGA-pancreatic cancer () dataset were downloaded and summarized.

Discussion

we found the top 289 most importance regions could take account of 80.0% contribution to the accurate prediction. Random forest algorithm showed 1585 regions could provide positive ability to distinguish cancer samples from normal samples while the top 286 most importance regions could take account of about 80.0% contribution, with sensitivity of 97.06%, specificity of 100% and accuracy of 97.37%. The second round prediction process of random forest model based on top 206 regions with mtry of 14 and mtrees of 500 showed 100% sensitivity, 100 specificity and 100% accuracy. One the other side, 248 regions were hypermethylated in at least 50% cancer samples. The average length of the regions were 103bp (IQR=95bp, SD=126) and were located in the promoter region of 182 genes (2000bp up-stream of TSS). With the help of text mining, 21 genes of them were validated to be methylation relevant cancer related genes (Table).

Next-generation methylation sequencing and quality control

Unique mappable reads

In the first step,


In the next step, genome-wide DNA methylation profiles of 30 samples including 15 solid cancer tissues (5 colon cancer, 5 lung cancer and 5 pancreatic cancer) and corresponding plasmas detected with reduced representation bisulfite sequencing (RRBS) were collected to discover the pattern of the shedding for the methylated DNA fragments from tissues to blood.


Task 1.

1, You need to filter all the hypermethylated fragement in cancer solid tissues and corresponding plasma circulating DNA while no methylation signals in health plasma. 2, and then validation these signals in TCGA database.


Task 2. Different method comparison.

Can it be used on RRBS data? We can downsample high coverage data to lower the effect of clonal reads on the analysis. Noi to check the input requirement for performing RRBS, can we perform RRBS in parallel on the test samples.

The optimal number of variables tried at each split is 16 (mtry) and number of trees (mtrees) makes no difference to the prediction accuracy from 500 to 7000. RF prediction model showed 1585 regions could provide positive ability to distinguish cancer samples from normal samples with high specificity of 100%, however, the sensitivity was only 25%, which indicated large number low predictive biomarkers were enrolled into the prediction model. After removed the most 60% lower informative regions in the random forest model the sensitivity was only 78%, (Informative regions from BSPP see supplementary Table *)

Therefore, we merged with the informative biomarkers identified in RRBS and Cap-seq dataset, eventual, there are 2 biomarkers were remained, including:

In the prediction section, top 62 regions could provide the distinguish accuracy of 86.87% which did not indicated Capseq method were worse than RRBS in the prediction ability. First capture, On the another side, when you will top 10 predictive regions, the separate ability of prediction model based on Capseq dataset could come up to the sensitivity of 90.32%, specificity of 95.83% and accuracy of 92.49%. However, small number of predictors would greatly decrease the robust or reproducibility of the prediction model and would bring biased inference in the process of clinical application of Capseq assay.

63.4% MHL distance of BSPP were less than 200bp, indicating they are almost located in the same region/CpG island in human genome, therefore, we need merge these MHL region together to increase the sensitivity of the prediction.


Method and Materials

Clinical sample and DNA collection

  • NSCLC samples and corresponding normal lung tissues for validation study in Chinese population were obtained from 150 patients who underwent pulmonary resection for primary NSCLC at Changhai Hospital, Shanghai, China. The study was approved by Fudan University and Changhai Hospital and Informed consents were obtained from the patients. Exclusion criteria included subjects with a family history of lung cancer, previous radiotherapy, and chemotherapy or adjuvant therapy before surgery. All tissues were immediately frozen at -80℃ after surgical resection. Histological examination and tumor-node-metastasis classification were conducted according to World Health Organization classification criteria (5) and AJCC Cancer Staging Manual, 7th Edition (6), respectively. Age, gender, smoking status, histology type, TNM stage and differentiation status were collected as the covariates when conducting association between DNA methylation and disease status. Smoking status was assigned to binary status: never and ever smoking. TNM stage was assigned to early stage (I and II) or late stage (III and IV) when it is necessary so that the sample size can be big enough to get the efficient statistic power.
  • CTT-frozen and CTT-FFPE RRBS library and sequencing [2]
    • The aim of expriments on TT-frozen and CTT-FFPE
  • N37 samples with WGBS

Bam files were stored in TSCC server:

Mapping, bisulfite conversion rate

  • First batch RRBS alignment statistic[3]
  • Second batch RRBS alignment statistic[4]

Average Methyaltion Level for Specific Genomic Regions

Obviously, we need compare the effect of MHL and average methylation level, especially in the same genomic regions. Here is how to get average methylation level for some genomic postions.

PileOMeth extract -r chr10:123-456 genome.fa alignments.bam

Pairwise R2 (LD) calcualtion within the methylation block from haploinfo files

  • Calculate the pairwise R2 with perl script code
  1. prepare script to calculate LD R2 by genomic region Shicheng:haploinfo2LDR2.pl
  2. collect WGBS(WB,ES,Roadmap) haploinfo files.
  3. R script for LD block plot code
  4. relationship between R2 and distance (absolute and relative) code
  5. Haplotype and MHL

Number of MHBs with different R2 cut-off

cd /home/shg047/oasis/monod/hapinfo/WGBS
qsub hapinfo2mhb.job

Umi-BSPP and Methylation Haplotype Loading (MHL)

Bioinformatics and Statistics

We need to look into each cancer individually on the list, maybe rank the cancer samples from most preferable to least preferable. Once we get optimistic results, we can further ask for samples from other centers. We cannot realistically ask specifically for each stage of cancer, although having samples from multiple stages might be better. prevalence treatment options for early detection (stage 1) current diagnostic methods We need to ask for buffy coat (1 tube) and matched serum samples (2 tubes) per patient We need to ask for primary tumor samples, they might be Formalin-fixed paraffin-embedded (FFPE) samples check Blueprint data to see how successful we were at capturing DNA purified from FFPE samples look into the kit for fixing degraded DNA Noi said we successfully captured with 50 ng converted DNA before, and it is possible to obtain this much from 2 tubes of serum samples, however, we can also look into the amplified DNA from Illumina (can we perform capture on these samples?)

RRBS option Dinh need to check the bayes classification algorithm Can it be used on RRBS data? We can downsample high coverage data to lower the effect of clonal reads on the analysis. Noi to check the input requirement for performing RRBS, can we perform RRBS in parallel on the test samples.


In the biomarker discover stage, the biomarker whose methylation were 100% un-methylated in normal plasmas were enrolled while any regions which were detected to be methylated in normal plasmas were filter out from the candidates. In the procedure of differential methylation test, regions whose variance in total samples were stage at lowest 30% quantile were removed to decrease the burden of multi-test correction. For the random forest prediction, optimal parameters were tuned before the tanning of the model the best number of tries and tress were determined by the grid search method with lowest out of beg prediction error. Intersection analysis of the genome position were conducted by BEDTOOLS and the regions whose distance less than 25bp were taken as the same biomarker region.

R packages of IRanges, Biostrings, stringr, randomForest, impute, rpart, e1071, biclust were used in the the statistic and bioinformatics analysis.

Potential biomarker based on published papers from NCBI were extracted for lung cancer (74 papers), colon cancer (8 papers) and pancreatic cancer (13 papers) with the strategy of title including “methylation” and corresponding cancer symbols as well as abstract including “Diagnosis”.

Gene Ontology Analysis

R code


Cluster and Classification Analysis

1360 tissue specific MHL in tissue cluster

Cancer specific methylation haplotype

Analysis Details

Dr. Zhang hope to visualize some explicit cancer specific methylation haplotype. I checked the wiki page for his previous anlaysis and conducted similar analysis based on paired cancer tissues, cancer plasma and normal plasma.

wiki page:[[5]]


Prepare SRA submission

cd /home/shg047/oasis/monod/hapinfo/SRA
cp ../*-P-* ./
cp ../*-T*  ./
cp ../*N37* ./

Acknowledgement

Author’s Contribution

Can it be used on RRBS data? We can downsample high coverage data to lower the effect of clonal reads on the analysis. Noi to check the input requirement for performing RRBS, can we perform RRBS in parallel on the test samples.

Reference

1.Zhao Y, Sun J, Zhang H, Guo S, Gu J, Wang W, et al. High-frequency aberrantly methylated targets in pancreatic adenocarcinoma identified via global DNA methylation analysis using methylCap-seq. Clinical epigenetics. 2014;6(1):18. Epub 2014/10/03.
2.Siegel RL, Miller KD, Jemal A. Cancer statistics, 2015. CA: a cancer journal for clinicians. 2015;65(1):5-29. Epub 2015/01/07.
3.Hankey BF, Ries LA, Edwards BK. The surveillance, epidemiology, and end results program: a national resource. Cancer Epidemiol Biomarkers Prev. 1999;8(12):1117-21. Epub 1999/12/29.
4.Warton K, Lin V, Navin T, Armstrong NJ, Kaplan W, Ying K, et al. Methylation-capture and Next-Generation Sequencing of free circulating DNA from human plasma. BMC genomics. 2014;15:476. Epub 2014/06/16.
5.Gibbs AR, Thunnissen FB. Histological typing of lung and pleural tumours: third edition. J Clin Pathol. 2001;54(7):498-9. Epub 2001/06/29.
6.Edge SB, Compton CC. The American Joint Committee on Cancer: the 7th edition of the AJCC cancer staging manual and the future of TNM. Ann Surg Oncol. 2010;17(6):1471-4. Epub 2010/02/25.
7.Desnoyers LR, Pai R, Ferrando RE, Hotzel K, Le T, Ross J, et al. Targeting FGF19 inhibits tumor growth in colon cancer xenograft and FGF19 transgenic hepatocellular carcinoma models. Oncogene. 2008;27(1):85-97. Epub 2007/06/30.
8.Revill K, Wang T, Lachenmayer A, Kojima K, Harrington A, Li J, et al. Genome-wide methylation analysis and epigenetic unmasking identify tumor suppressor genes in hepatocellular carcinoma. Gastroenterology. 2013;145(6):1424-35 e1-25. Epub 2013/09/10.
9.Kitagawa N, Ojima H, Shirakihara T, Shimizu H, Kokubu A, Urushidate T, et al. Downregulation of the microRNA biogenesis components and its association with poor prognosis in hepatocellular carcinoma. Cancer science. 2013;104(5):543-51. Epub 2013/02/13.
10.Luo FY, Xiao S, Liu ZH, Zhang PF, Xiao ZQ, Tang CE. Kank1 reexpression induced by 5-Aza-2'-deoxycytidine suppresses nasopharyngeal carcinoma cell proliferation and promotes apoptosis. International journal of clinical and experimental pathology. 2015;8(2):1658-65. Epub 2015/05/15.
11.Liu YX, Zhang SF, Ji YH, Guo SJ, Wang GF, Zhang GW. Whole-exome sequencing identifies mutated PCK2 and HUWE1 associated with carcinoma cell proliferation in a hepatocellular carcinoma patient. Oncology letters. 2012;4(4):847-51. Epub 2012/12/04.