计算生物学期末作业张伟伟综述.docx

1、计算生物学期末作业张伟伟综述利用DNA甲基化450K array数据估计肿瘤细胞纯度与差异甲基化分析计算生物学编程语言考试题目 DNA甲基化是一类重要的表观遗传修饰，它可以通过控制基因表达在多种细胞过程中起作用，比如胚胎干细胞发育 1，基因组印迹化 2, 3，X染色体失活 4 等等。其异常表达与包括癌症在内的多种人类疾病密切相关 5-12。在单核苷酸水平上得到DNA甲基化水平的试验方法一般有两类：一类是基于重亚硫酸盐处理的方法，即“全基因组重亚硫酸盐测序”（WGBS）和“约化表示重亚硫酸盐测序”（RRBS）。这两类方法都是在高通量测序之前用重亚硫酸盐对DNA序列进行处理，没有被甲基化的胞嘧啶（

2、cytosine）被转化为了尿嘧啶（uracil），而甲基化的胞嘧啶保持不变。因此，通过回帖到基因组上并与参考基因组相比较，就可以得知相应位置上的胞嘧啶是否被甲基化。另一类获取DNA甲基化的方法是传统的450K 芯片办法，它可检测人全基因组近450,000个甲基化位点，具有单碱基的分辨率。全面覆盖了96%的CpG岛，并根据需求加入了CpG岛以外的CpG位点、人类干细胞非CpG甲基化位点、正常组织与肿瘤（多种癌症）组织差异甲基化位点、编码区以外的CpG岛、miRNA启动子区域和GWAS疾病相关区域的位点。 在癌症表观遗传分析当中，一个很重要的任务就是获取肿瘤细胞与正常细胞相比发生甲基化差异的位置

3、，即差异甲基化分析。实验上，我们一般采取的方式是在切除肿瘤组织的同时，再切除一部分癌旁正常组织一并测序，然后进行比较。然而，由于手术分辨率的问题，切出的肿瘤组织都或多或少地混有一部分正常细胞，比如TCGA的所有肿瘤样本，正常细胞的含量一般都在30%到70%之间。因此，如果我们不加以处理，而直接把它跟正常组织的甲基化谱相比，势必会出现误差。此外，还有大部分样本是没有正常癌旁组织的，比如脑癌（GBM），这种情况下我们只能把它跟有正常组织的样本相比，这样，一些个体特异性的差异甲基化区域就检测不到了。 为了克服上述困难，我们曾利用极大似然估计和EM算法对WGBS和RRBS数据进行了分解，得到的肿瘤细胞

4、纯度和差异甲基化区域（DMR）跟用“正常-肿瘤”组织相比得到的结果类似，这部分工作发表在8月份的Genome Biology上 13. 然而与WGBS和RRBS数据相比，450K的数据获取更便宜、使用更为广泛，因此，如果我们能提出一种类似的方法估计450K的tumor purity和DMR，无论在理论上还是临床实际上都会有更为重要的意义。当然，思路也不能局限在 GB的paper上，如果是一样的话，我也就顺手做完了。对于这个问题我有如下初步的思路供大家参考：1）对于成对样本，即正常和肿瘤组织都有的样本，我们可以先找到正常和肿瘤组织中有明显差异的甲基化位点，对于tumor purity的估计，这样

5、的位点才是有信息的。对于纯肿瘤细胞，如果我们假设其在每个CpG位点的甲基化水平或者是0，或者是1，那么很容易就可以利用肿瘤组织单个CpG位点的甲基化水平来估计肿瘤细胞的纯度。但很显然，由于测量误差等等因素，利用单个CpG位点来预测肿瘤细胞纯度肯定不会太准确，这里我们可以用所有的有信息的位点一起来估计。TCGA肺腺癌（LUAD）有利用ABSOLUTE软件估计出来的肿瘤细胞纯度，你可以把你估计出来的结果跟ABSOLUTE的结果比较一下，看看是否有很强的相关性？ 2）对于没有癌旁正常组织的肿瘤样本，我们只能退而求其次只利用肿瘤组织的450K数据来估计其tumor purity. 在这里，我们还是假设

6、纯的正常细胞和肿瘤细胞的甲基化水平只取0（全部不甲基化）或者1（全部甲基化）。这样，如果我们知道了肿瘤和正常细胞各自的甲基化水平，那么根据肿瘤样本的450K数据（肿瘤细胞核正常细胞的混合）就可以推算出tumor purity了。但这个信息我们是不知道的，那该怎么办呢？对了，我们可以把这个因素设为隐藏变量，用EM算法就可以解决啦。但即使你成功了，还是有一个问题没有解决，即“相信息”会丢失：你不知道与正常细胞相比，肿瘤细胞在这个位置上是高甲基化了，还是低甲基化了？这个问题我们可以预先分析下TCGA的成对数据，看看是否有一些区域是在所有样本里面都是稳定地“高甲基化”或者“低甲基化”的位点（我想一定会

7、有这样的位点的，而且这些位点应该与肿瘤发生密切相关！）。把这个方法估计出来的结果与ABSOLUTE的结果比较一下，看看相关性有多大？ 3）如果再考虑到肿瘤细胞容易发生拷贝数的变化呢？这个模型应该怎么改进？ 作业要求：1. 把分析过程和计算结果尽可能详细地写出来；2. 把计算的程序附上（python做这个是最合适的）；3. 多画图来展示你的中间结果；4. 放假之前把程序和结果交给我。 1. Li E, Bestor TH, Jaenisch R: Targeted mutation of the DNA methyltransferase gene results in embryonic le

8、thality. Cell 1992, 69(6):915-926.2. Li E, Beard C, Jaenisch R: Role for DNA methylation in genomic imprinting. Nature 1993, 366(6453):362-365.3. Fang F, Hodges E, Molaro A, Dean M, Hannon GJ, Smith AD: Genomic landscape of human allele-specific DNA methylation. Proc Natl Acad Sci U S A 2012, 109(19

9、):7332-7337.4. Panning B, Jaenisch R: RNA and the epigenetic regulation of X chromosome inactivation. Cell 1998, 93(3):305-308.5. Feinberg AP, Cui H, Ohlsson R: DNA methylation and genomic imprinting: insights from cancer into epigenetic mechanisms. Semin Cancer Biol 2002, 12(5):389-398.6. Ehrlich M

10、: DNA methylation in cancer: too much, but also too little. Oncogene 2002, 21(35):5400-5413.7. Jones PA, Baylin SB: The fundamental role of epigenetic events in cancer. Nat Rev Genet 2002, 3(6):415-428.8. Das PM, Singal R: DNA methylation and cancer. J Clin Oncol 2004, 22(22):4632-4642.9. Robertson

11、KD: DNA methylation and human disease. Nat Rev Genet 2005, 6(8):597-610.10. Beck S, Rakyan VK: The methylome: approaches for global DNA methylation profiling. Trends Genet 2008, 24(5):231-237.11. Javierre BM, Fernandez AF, Richter J, Al-Shahrour F, Martin-Subero JI, Rodriguez-Ubreva J, Berdasco M, F

12、raga MF, OHanlon TP, Rider LG et al: Changes in the pattern of DNA methylation associate with twin discordance in systemic lupus erythematosus. Genome Res 2010, 20(2):170-179.12. Hansen KD, Timp W, Bravo HC, Sabunciyan S, Langmead B, McDonald OG, Wen B, Wu H, Liu Y, Diep D et al: Increased methylati

13、on variation in epigenetic domains across cancer types. Nat Genet 2011, 43(8):768-775.13. Zheng X, Zhao Q, Wu H-J, Li W, Wang H, Meyer CA, Qin QA, Xu H, Zang C, Jiang P: MethylPurify: tumor purity deconvolution and differential methylation detection from single tumor DNA methylomes. Genome Biol 2014

14、, 15:419. 问题一：假设在纯的细胞系中，单个位点的甲基化水平或者是0，或者是1对于cancer和normal成对的样本，利用秩和检验选取有信息的cpg位点，建立线性统计模型 ，其中表示肿瘤细胞的纯度，表示cancer中第i个位点的甲基化水平，表示normal中第i个位点的甲基化水平，表示cancer中第i个位点的甲基化水平 ；当时，否则。应用此模型对于每一个位点均得到纯度，再取均值即可到整个样本的纯度，将此方法应用每一个样本，将得到的结果与absolute作比较，得到下图结果：图1：线性统计模型估计肿瘤纯度与absolute估计纯度的相关系数程序：setwd(/mnt/Storage/

15、home/zhengxq/wwzhang/homework/processed/paired_cancer)#read cancer and make a dataframedir-/mnt/Storage/home/zhengxq/wwzhang/homework/luaddata/cancerpair/file.names-list.files(dir)n-length(file.names)cpg-read.delim(paste(dir,file.names1,sep=),header=T),1-1cpg-as.vector.factor(cpg)cancer-c()for(i in

16、1:n) dat-read.delim(paste(dir,file.namesi,sep=),header=T) a-dat,2-1 a-as.vector.factor(a) a-as.numeric(a) cancer-cbind(cancer,a) colnames(cancer)i-strsplit(file.namesi,.,fixed=TRUE)16rownames(cancer)-cpgsave(cancer,file=cancer.RData)# read normal and make a dataframedir1-/mnt/Storage/home/zhengxq/ww

17、zhang/homework/luaddata/normalpair/file-list.files(dir1)m-length(file)normal-c()for(i in 1:m) dat-read.delim(paste(dir1,filei,sep=),header=T) a-dat,2-1 a-as.vector.factor(a) a-as.numeric(a) normal-cbind(normal,a) colnames(normal)i-strsplit(filei,.,fixed=TRUE)16 rownames(normal)-cpgsave(normal,file=n

18、ormal.RData)#选取秩和检验中pvalue10-6的有差异的CPG位点setwd(/home/users/wwzhang/wwzhang/luadprocessed/purity_wilcox_paired_mean)load(/home/users/wwzhang/wwzhang/luadprocessed/normal.RData)load(/home/users/wwzhang/wwzhang/luadprocessed/cancer.RData)tnormal-t(na.omit(normal)tcancer-t(na.omit(cancer)index-intersect(

19、colnames(tcancer),colnames(tnormal)tcancer1-tcancer,indextnormal1-tnormal,indexn-ncol(tcancer1)value-c()for(i in 1:n) value-cbind(value,wilcox.test(tcancer1,i,tnormal1,i,paired=T)$p.value) colnames(value)i-colnames(tcancer1)iindex-which(value10-6)posi-colnames(value)indexmcancer-cancerposi,mnormal-n

20、ormalposi,#利用线性统计模型求cancer的纯度#purity-function(x,y,i) s1-x,i d1-y,i n-length(s1) chundud1j) p-(s1j-d1j)/(1-d1j) else p-(d1j-s1j)/d1j chundu-append(chundu,p) pingjun-mean(chundu) return(pingjun)pur-c()for(i in 1:29) a-purity(mcancer,mnormal,i) pur-cbind(pur,a) colnames(pur)i-colnames(mcancer)i save(pu

21、r,file=purity_wilcox_paired_mean_0.000001.RData)#计算线性模型估计的纯度与absolute估计纯度的相关性#load absolut purityabsolute-read.delim(nature.tumor_purity,header=T,sep=)purity-absolute,3purity-as.data.frame(purity)rownames(purity)-absolute,1load(purity_wilcox_paired_mean_0.000001.RData )pur-t(pur)for(i in 1:29) rowna

22、mes(pur)i-substr(rownames(pur)i,1,12)name-as.vector.factor(absolute,1)index-intersect(rownames(pur),rownames(purity)purnew-as.data.frame(purindex,)puritynew-as.data.frame(purityindex,)rownames(puritynew)-indexcor-cor(purnew,puritynew)cor-round(as.vector(cor),3)tag-paste0(R=,cor)pdf(purity_mean.pdf)m

23、-lm(puritynew,1purnew,1)plot(purnew,1,puritynew,1,xlab=luadpurity,ylab=absolute purity)abline(m)text(0.2, 0.6, tag)dev.off()问题二：对于没有癌旁没有正常组织的样本，即没有对照组的cancer采取如下处理：（1）选取没有癌旁正常组织的所有cancer将其与已有的normal做秩和检验，选取pvalue10-6的位点。（2）将已有的normal取peak作为该位点的甲基化水平（3）将所有的cancer与（2）得到的一个样本应用问题一的线性统计模型，求得没有癌旁正常组织的can

24、cer的纯度，再将其与absolute估计出来的纯度作比较，结果如下：图2：没有癌旁正常组织的cancer纯度与absolute估计纯度的相关系数程序：#load normal samplesetwd(/mnt/Storage/home/zhengxq/wwzhang/homework/cancer_unpaired)load(mnormal_wilcox.test_paired_10-6.RData)#将所有的normal按位点取peak作为该位点的甲基化水平n-nrow(mnormal)methy-c()for(i in 1:n) alpha = density(mnormali,) x

25、= alpha$x y = alpha$y a= xwhich.max(y) methy-cbind(methy,a) colnames(methy)i-rownames(mnormal)i# read cancer sampledir-/mnt/Storage/home/zhengxq/wwzhang/homework/luaddata/cancer/file.names-list.files(dir)m-length(file.names)cpg-read.delim(paste(dir,file.names1,sep=),header=T,skip=1),1cpg-as.characte

26、r(cpg)# read 427 unpaired cancer cancer-c()for(i in 1:m)dat-read.delim(paste(dir,file.namesi,sep=),header=T,skip=1)a-dat,2cancer-cbind(cancer,a)colnames(cancer)i-strsplit(file.namesi,.,fixed=TRUE)16rownames(cancer)-cpg# make 427 cancer sample dataframe match normalcancer-na.omit(cancer)index-interse

27、ct(rownames(cancer),rownames(t(methy)cancerlast-cancerindex,methylast-t(methy)index,# compute cancer purityF.purity-function(x,y,i) s1-x,i d1-y n-length(s1) chundud1j) p-(s1j-d1j)/(1-d1j) else p-(d1j-s1j)/d1j chundu-append(chundu,p) purity-mean(chundu) return(purity)m=ncol(cancerlast)pur-c()for(i in

28、 1:m) a-F.purity(cancerlast,methylast,i) pur-cbind(pur,a) colnames(pur)i-colnames(cancerlast)i save(pur,file=luad-purity_mean.RData)# luad and absolute correlationabsolute-read.delim(nature.tumor_purity,header=T,sep=)purity-absolute,3purity-as.data.frame(purity)rownames(purity)-absolute,1load(luad-purity_mean.RData )pur-t(pur)for(i in 1:427) rownames(pur)i-substr(rownames(pur)i,1,12)name-as.vector.factor(absolute,1)index-intersect(rownames(pur),rownames(purity)purnew-as.data.frame(purindex,)