计算生物学期末作业张伟伟

利用DNA甲基化450Karray数据估计肿瘤细胞纯度与差异甲基化分析——《计算生物学编程语言》考试题目DNA甲基化是一类重要的表观遗传修饰，它可以通过控制基因表达在多种细胞过程中起作用，比如胚胎干细胞发育[1]，基因组印迹化[2,3]，X染色体失活[4]等等。其异常表达与包括癌症在内的多种人类疾病密切相关[5-12]。在单核苷酸水平上得到DNA甲基化水平的试验方法一般有两类：一类是基于重亚硫酸盐处理的方法，即“全基因组重亚硫酸盐测序”（WGBS）和“约化表示重亚硫酸盐测序”（RRBS）。这两类方法都是在高通量测序之前用重亚硫酸盐对DNA序列进行处理，没有被甲基化的胞嘧啶（cytosine）被转化为了尿嘧啶（uracil），而甲基化的胞嘧啶保持不变。因此，通过回帖到基因组上并与参考基因组相比较，就可以得知相应位置上的胞嘧啶是否被甲基化。另一类获取DNA甲基化的方法是传统的450K芯片办法，它可检测人全基因组近450,000个甲基化位点，具有单碱基的分辨率。全面覆盖了96%的CpG岛，并根据需求加入了CpG岛以外的CpG位点、人类干细胞非CpG甲基化位点、正常组织与肿瘤（多种癌症）组织差异甲基化位点、编码区以外的CpG岛、miRNA启动子区域和GWAS疾病相关区域的位点。在癌症表观遗传分析当中，一个很重要的任务就是获取肿瘤细胞与正常细胞相比发生甲基化差异的位置，即差异甲基化分析。实验上，我们一般采取的方式是在切除肿瘤组织的同时，再切除一部分癌旁正常组织一并测序，然后进行比较。然而，由于手术分辨率的问题，切出的肿瘤组织都或多或少地混有一部分正常细胞，比如TCGA的所有肿瘤样本，正常细胞的含量一般都在30%到70%之间。因此，如果我们不加以处理，而直接把它跟正常组织的甲基化谱相比，势必会出现误差。此外，还有大部分样本是没有正常癌旁组织的，比如脑癌（GBM），这种情况下我们只能把它跟有正常组织的样本相比，这样，一些个体特异性的差异甲基化区域就检测不到了。为了克服上述困难，我们曾利用极大似然估计和EM算法对WGBS和RRBS数据进行了分解，得到的肿瘤细胞纯度和差异甲基化区域（DMR）跟用“正常-肿瘤”组织相比得到的结果类似，这部分工作发表在8月份的GenomeBiology上[13].然而与WGBS和RRBS数据相比，450K的数据获取更便宜、使用更为广泛，因此，如果我们能提出一种类似的方法估计450K的tumorpurity和DMR，无论在理论上还是临床实际上都会有更为重要的意义。当然，思路也不能局限在GB的paper上，如果是一样的话，我也就顺手做完了。对于这个问题我有如下初步的思路供大家参考：1）对于成对样本，即正常和肿瘤组织都有的样本，我们可以先找到正常和肿瘤组织中有明显差异的甲基化位点，对于tumorpurity的估计，这样的位点才是有信息的。对于纯肿瘤细胞，如果我们假设其在每个CpG位点的甲基化水平或者是0，或者是1，那么很容易就可以利用肿瘤组织单个CpG位点的甲基化水平来估计肿瘤细胞的纯度。但很显然，由于测量误差等等因素，利用单个CpG位点来预测肿瘤细胞纯度肯定不会太准确，这里我们可以用所有的有信息的位点一起来估计。TCGA肺腺癌（LUAD）有利用ABSOLUTE软件估计出来的肿瘤细胞纯度，你可以把你估计出来的结果跟ABSOLUTE的结果比较一下，看看是否有很强的相关性？2）对于没有癌旁正常组织的肿瘤样本，我们只能退而求其次只利用肿瘤组织的450K数据来估计其tumorpurity.在这里，我们还是假设纯的正常细胞和肿瘤细胞的甲基化水平只取0（全部不甲基化）或者1（全部甲基化）。这样，如果我们知道了肿瘤和正常细胞各自的甲基化水平，那么根据肿瘤样本的450K数据（肿瘤细胞核正常细胞的混合）就可以推算出tumorpurity了。但这个信息我们是不知道的，那该怎么办呢？对了，我们可以把这个因素设为隐藏变量，用EM算法就可以解决啦。但即使你成功了，还是有一个问题没有解决，即“相信息”会丢失：你不知道与正常细胞相比，肿瘤细胞在这个位置上是高甲基化了，还是低甲基化了？这个问题我们可以预先分析下TCGA的成对数据，看看是否有一些区域是在所有样本里面都是稳定地“高甲基化”或者“低甲基化”的位点（我想一定会有这样的位点的，而且这些位点应该与肿瘤发生密切相关！）。把这个方法估计出来的结果与ABSOLUTE的结果比较一下，看看相关性有多大？3）如果再考虑到肿瘤细胞容易发生拷贝数的变化呢？这个模型应该怎么改进？作业要求：1.把分析过程和计算结果尽可能详细地写出来；2.把计算的程序附上（python做这个是最合适的）；3.多画图来展示你的中间结果；4.放假之前把程序和结果交给我。1.LiE,BestorTH,JaenischR:TargetedmutationoftheDNAmethyltransferasegeneresultsinembryoniclethality.Cell1992,69(6):915-926.2.LiE,BeardC,JaenischR:RoleforDNAmethylationingenomicimprinting.Nature1993,366(6453):362-365.3.FangF,HodgesE,MolaroA,DeanM,HannonGJ,SmithAD:Genomiclandscapeofhumanallele-specificDNAmethylation.ProcNatlAcadSciUSA2012,109(19):7332-7337.4.PanningB,JaenischR:RNAandtheepigeneticregulationofXchromosomeinactivation.Cell1998,93(3):305-308.5.FeinbergAP,CuiH,OhlssonR:DNAmethylationandgenomicimprinting:insightsfromcancerintoepigeneticmechanisms.SeminCancerBiol2002,12(5):389-398.6.EhrlichM:DNAmethylationincancer:toomuch,butalsotoolittle.Oncogene2002,21(35):5400-5413.7.JonesPA,BaylinSB:Thefundamentalroleofepigeneticeventsincancer.NatRevGenet2002,3(6):415-428.8.DasPM,SingalR:DNAmethylationandcancer.JClinOncol2004,22(22):4632-4642.9.RobertsonKD:DNAmethylationandhumandisease.NatRevGenet2005,6(8):597-610.10.BeckS,RakyanVK:Themethylome:approachesforglobalDNAmethylationprofiling.TrendsGenet2008,24(5):231-237.11.JavierreBM,FernandezAF,RichterJ,Al-ShahrourF,Martin-SuberoJI,Rodriguez-UbrevaJ,BerdascoM,FragaMF,O'HanlonTP,RiderLGetal:ChangesinthepatternofDNAmethylationassociatewithtwindiscordanceinsystemiclupuserythematosus.GenomeRes2010,20(2):170-179.12.HansenKD,TimpW,BravoHC,SabunciyanS,LangmeadB,McDonaldOG,WenB,WuH,LiuY,DiepDetal:Increasedmethylationvariationinepigeneticdomainsacrosscancertypes.NatGenet2011,43(8):768-775.13.ZhengX,ZhaoQ,WuH-J,LiW,WangH,MeyerCA,QinQA,XuH,ZangC,JiangP:MethylPurify:tumorpuritydeconvolutionanddifferentialmethylationdetectionfromsingletumorDNAmethylomes.GenomeBiol2014,15:419.问题一：假设在纯的细胞系中，单个位点的甲基化水平或者是0，或者是1对于cancer和normal成对的样本，利用秩和检验选取有信息的cpg位点，建立线性统计模型(i)(i)(i)12(1)mmm，其中表示肿瘤细胞的纯度，(i)1m表示cancer中第i个位点的甲基化水平，(i)2m表示normal中第i个位点的甲基化水平，(i)m表示cancer中第i个位点的甲基化水平；当(i)(i)2mm时，(i)11m，否则(i)10m。应用此模型对于每一个位点均得到纯度，再取均值即可到整个样本的纯度，将此方法应用每一个样本，将得到的结果与absolute作比较，得到下图结果：图1：线性统计模型估计肿瘤纯度与absolute估计纯度的相关系数程序：setwd(/mnt/Storage/home/zhengxq/wwzhang/homework/processed/paired_cancer)##readcancerandmakeadataframedir-/mnt/Storage/home/zhengxq/wwzhang/homework/luaddata/cancerpair/file.names-list.files(dir)n-length(file.names)cpg-read.delim(paste(dir,file.names[1],sep=),header=T)[,1][-1]cpg-as.vector.factor(cpg)cancer-c()for(iin1:n){dat-read.delim(paste(dir,file.names[i],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.names[i],.,fixed=TRUE)[[1]][6]}rownames(cancer)-cpgsave(cancer,file=cancer.RData)##readnormalandmakeadataframedir1-/mnt/Storage/home/zhengxq/wwzhang/homework/luaddata/normalpair/file-list.files(dir1)m-length(file)normal-c()for(iin1:m){dat-read.delim(paste(dir1,file[i],sep=),header=T)a-dat[,2][-1]a-as.vector.factor(a)a-as.numeric(a)normal-cbind(normal,a)colnames(normal)[i]-strsplit(file[i],.,fixed=TRUE)[[1]][6]}rownames(normal)-cpgsave(normal,file=normal.RData)###########################################################################选取秩和检验中pvalue10^-6的有差异的CPG位