生物芯片—ChAMP包学习笔记

一、甲基化芯片简介

一张芯片包括12个array，也就是一张芯片可以做12个sample，一台机子可以跑8张芯片，也就是一共96个sample，每个样本可以测到超过450，000个CpG位点的甲基化信息（大概人所有的1%，但是覆盖了多数CpG岛和启动子区），芯片本身包含一些控制探针可以做质控。

芯片注释信息：

#查看450k芯片注释信息
library(IlluminaHumanMethylation450kmanifest)
show(IlluminaHumanMethylation450kmanifest)

#查看850k芯片信息使用另一个R包
library(IlluminaHumanMethylationEPICkmanifest)

minfi包中有特定的函数getManifest()可以查看指定对象对应芯片的注释信息

二、ChAMP包简介

ChAMP可以用来处理450k和EPIC的数据，包中整合了质量控制、数据标准化、识别DMR、识别CNA等函数，使甲基化数据的处理形成一个pipeline。

如果你的整个流程中参数都设置为默认值，那么你只需运行函数：champ.process(directory = testDir)，参数directory为甲基化原始数据.IDAT所在文件目录。

当然，ChAMP中的每个函数也可以单独运行，自定义参数。

三、具体处理流程

（一）数据准备

ChAMP包读取数据不仅要.idat原始文件，还需要一个pd(phenotype)文件，名为sample_sheet.csv（当然sample_sheet这个名字不重要，重要的是这是一个csv文件），具体信息如下：

其中比较重要的属性列是Sample_Name、Sample_Group、Slide(Sentrix_ID)、Array(Sentrix_position)

Array（Sentrix_position）：一张甲基化芯片上最多有12个样本，每个样本根据Sentrix_position标识；
Slide（Sentrix_ID）：当样本个数大于12个时，必然需要另外一张芯片，对于每张芯片，使用Sentrix_ID标识；
minfi就是通过Sentrix_ID和Sentrix_position这两个字段来查找样本的原始数据

需要根据各样本信息自行创建一个sample_sheet.csv，如何构建？Sentrix_ID 和 Sentrix_Position来自IDAT文件名。
注意：GEO下载的数据命名是提交者自己命名，难以构建该csv文件。TCGA下载数据，文件名包含所需信息，可以构建该csv文件

1.使用ChAMPdata包中的数据

library(ChAMP)
testdir<-system.file("extdata",package="ChAMPdata")

去上述路径查看数据：

testDir[1]=”/public/workspace/shaojf/anaconda3/lib/R/library/ChAMPdata/extdata”

包含8个样本，4个肿瘤样本，4个对照样本

2.数据库下载，以GEO为例说明数据准备

（二）数据读取与过滤champ.load

myLoad <-champ.load(testDir,arraytype="450K")
#芯片类型，可设置为"EPIC"用于加载850K芯片的数据

champ.load()参数见champ.import()和champ.filter()，champ.load()利用minfi包来读入数据
champ.load()==champ.import()+champ.filter()
返回一个列表包含3个元素：“beta”,“intensity”,“pd”

beta是一个矩阵，列代表样本，行代表探针
pd记录样本信息，即自己编写的sample_sheet.csv文件，重要的属性列”Slide”,“Array”

读入数据后有6种数据过滤:

① filterDetP、detPcut：默认过滤掉<0.01的探针，使用detPcut来修改阈值

②filterBeads、beadCutoff：探针与DNA序列杂交，beadcounts就是与探针结合的reads数，若结合得较少，则认为该种探针信号不可靠，实际处理中，默认如果这个探针在5%的样本中，beadcount<3这个探针就会被过滤掉，使用beadCutoff来修改阈值

③ filterNoCG：过滤非CpG位点的探针，对甲基化的C而言，有CpG\CHG\CHH等甲基化的C，若不想过滤这些探针可以设置filterNoCG=FALSE

④ filterSNPs：过滤SNP附近的CpG位点的探针，SNP位点附近的CpG位点不稳定，设filterSNPs=FALSE不进行过滤

⑤filterMultiHit：过滤非特异性的探针，Nordlund等研究人员使用bwa将探针序列与基因组进行比对，若唯一比对，说明探针特异性好，根据比对结果，去除比对到基因组上多个位置的探针，修改参数filterMultiHit=FALSE不进行过滤

⑥filterXY：过滤性染色体上的探针，性染色体上的差异CpG是不准确的，无法解释这个差异是否是实验造成，若性别相同时，可以不用过滤，修改参数filterXY=FALSE不进行过滤

（三）可视化查看质量控制情况champ.QC

champ.QC(beta = myLoad$beta,
		 #---产生的β矩阵
		 pheno=myLoad$pd$Sample_Group,
		 #---样本分类信息
		 mdsPlot=TRUE,
		 densityPlot=TRUE,
		 dendrogram=TRUE,
         PDFplot=TRUE,Rplot=TRUE,
		 Feature.sel="None",
		 #---绘制dendrogram时所使用的计算方法
		 resultsDir="./CHAMP_QCimages/")
QC.GUI()

（四）数据标准化 champ.norm

对探针的矫正即消除两类探针因技术差异带来的误差

myNorm<-champ.norm(beta=myLoad$beta,rgSet=myLoad$rgSet,mset=myLoad$mset,
			resultsDir=".CHAMP_Normalization/",
			method="BMIQ",
			plotBMIQ=FALSE,
			arraytype="450K",
			cores=3)

（五）批次效应

1.查看是否有批次效应champ.SVD

champ.SVD(beta = myNorm,
		   rgSet=NULL,
		   pd=myLoad$pd,
		   RGEffect=FALSE,
		   #---If Green and Red color control probes would be calculated
          PDFplot=TRUE,Rplot=TRUE,
		  resultsDir="./CHAMP_SVDimages/")

产生2个图像：Singular Value Decomposition Analysis (SVD)；Scree Plot

①Singular Value Decomposition Analysis (SVD)
颜色越深，p-value越小，变异越显著，表明变异成分越大
从图像可以看出变异主要来自样本间(control组和tumor组)，即我们关注的生物学因素，但变异显著性不大，经过过滤后得到的数据比较clear，适合后续分析。
如果发现在array或slide等有较显著的p-value，则需要检查实验设计以及使用其他normalization methods来进行后续的批次效应矫正

②scree plot：

各主成分可解释的变异占总变异的百分比

2.校正批次效应champ.runCombat

mycombat <- champ.runCombat(beta=myNorm,
				pd=myLoad$pd,
                variablename="Sample_Group",
				#要进行矫正的协变量，参考sva包的ComBat函数
				batchname=c("Slide"),
				#批次，来自pd文件
                logitTrans=TRUE)
				#对beta是否进行logit转化，利用原始beta进行矫正时设为T，利用M值进行计算时设为F 

ComBat如果直接用beta值进行矫正，输出值可能不在0-1之间，所以在计算之前会进行一个logit转化成M-value；
若用M值进行矫正则没必要进行转化，因此设置logitTrans=F

（六）甲基化分析

1.champ.DMP

myDMP <- champ.DMP(beta = myNorm,
	  				pheno = myLoad$pd$Sample_Group,
      				compare.group = NULL,
	 	 			#当group>2时设定该参数为一个列表，其中每个元素为两个Sample_Group，进而指定要进行比较的group，否则每两个group之间都会进行比较
     				adjPVal = 0.05,
      				#设定的p-value检验水平
	  				adjust.method = "BH",
	  				#p-value矫正，使用的method参考函数p.adjust()
      				arraytype = "450K")
#直接通过基于类别group的探针比较计算得到DMP

# 可视化查看分布
DMP.GUI(DMP=myDMP[[1]],beta=myNorm,pheno=myLoad$pd$Sample_Group)
#使用GUI需要shiny包

2.champ.DMR

myDMR<-champ.DMR(beta=myNorm,
          pheno=myLoad$pd$Sample_Group,
          compare.group=NULL,
          arraytype="450K",
          method = "Bumphunter",
          minProbes=7,
          #---DMR包含很多CpG位点,每个位点都对应一个探针，minProbes设置每个DMR包含的最小探针数
          adjPvalDmr=0.05,
          #---DMR的p-value是用Stouffler's method计算得到，再进行矫正后得到adjPvalDmr
          cores=3,
          ## following parameters are specifically for Bumphunter method.
          maxGap=300,
          cutoff=NULL,
          pickCutoff=TRUE,
          smooth=TRUE,
          smoothFunction=loessByCluster,
          useWeights=FALSE,
          permutations=NULL,
          B=250,
          nullMethod="bootstrap",
          ## following parameters are specifically for probe ProbeLasso method.
          meanLassoRadius=375,
          minDmrSep=1000,
          minDmrSize=50,
          adjPvalProbe=0.05,
          Rplot=T,
          PDFplot=T,
          resultsDir="./CHAMP_ProbeLasso/",
          ## following parameters are specifically for DMRcate method.
          rmSNPCH=T,
          fdr=0.05,
          dist=2,
          mafcut=0.05,
          lambda=1000,
          C=2)

# 可视化查看分布
DMR.GUI(DMR=myDMR)

（七）找拷贝数变异CNAchamp.CNA

（八）GSEA分析champ.GSEA

今天的文章生物芯片—ChAMP包学习笔记分享到此就结束了，感谢您的阅读。