TCGA甲基化数据分析

cd /home/train/Liuping/raw_data 
mkdir ProstateCancer_methylation 
cd ProstateCancer_methylation/

#合并下载好的数据
nohup ./Merge.pl metadata.cart.2018-08-01.json > /home/train/Liuping/raw_data/TCGA_Methylation_ProstateCancer/methylayion_rawdata/out_methylation 2>&1 &
#注意记录下normal‘和 tumor的个数
# normal count: 50
# tumor count: 503

mkdir ../result
cp geneMethy.txt ../result/
cd ../result/

#差异甲基化
source("http://bioconductor.org/biocLite.R")
biocLite("limma")
library(limma)
setwd("D:/TCGA_methylation")
inputFile <- "geneMethy.txt"
normalNum=50
tumorNum=503
fdrFilter=0.05
logFCfilter=1

outTab <- data.frame()
grade <- c(rep(1,normalNum),rep(2,tumorNum))
Type <- c(rep("Normal",normalNum),rep("Tumor",tumorNum))
rt <- read.table(inputFile,sep="\t",header=T,check.names=F)
rt <- as.matrix(rt)
rownames(rt) <- rt[,1]
exp <- rt[,2:ncol(rt)] #行名可被直接移植
dimnames <- list(rownames(exp),colnames(exp))
data <- matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
data <- avereps(data)
data <- data[rowMeans(data)>0,]

data <- normalizeBetweenArrays(data)
normalData <- cbind(id=row.names(data),data) #为了让id号写入文件
write.table(normalData,file="normalizeMethy.txt",sep="\t",row.names=F,quote=F)

for (i in row.names(data)){
 rt <- rbind(expression=data[i,],grade=grade)
 rt <- as.matrix(t(rt))
 wilcoxTest <- wilcox.test(expression ~ grade, data=rt)

 normalGeneMeans <- mean(data[i,1:normalNum])
 tumorGeneMeans <- mean(data[i,(normalNum+1):ncol(data)])
 logFC <- log2(tumorGeneMeans) - log2(normalGeneMeans)
 pvalue <- wilcoxTest$p.value
 normalMed <- median(data[i,1:normalNum])
 tumorMed <- median(data[i,(normalNum+1):ncol(data)])
 diffMed <- tumorMed - normalMed
 if((( logFC>0) & ( diffMed >0)) | (( logFC < 0 ) & ( diffMed <0 ))){
 outTab <- rbind(outTab,cbind(gene=i,normalMean=normalGeneMeans,TumorMean=tumorGeneMeans,logFC=logFC,pValue=pvalue))
 }}
pValue <- outTab[,"pValue"]
fdr <- p.adjust(as.numeric(as.vector(pValue)),method="fdr")
outTab <- cbind(outTab,fdr=fdr)
write.table(outTab,file="allGene.xls",sep="\t",row.names=F,quote=F)
outDiff <- outTab[(abs(as.numeric(as.vector(outTab$logFC))) > logFCfilter & as.numeric(as.vector(outTab$fdr)) < fdrFilter),]
write.table(outDiff,file="diff.xls",sep="\t",row.names=F,quote=F)
heatmap <- data[as.vector(outDiff[,1]),]
heatmap <- cbind(id=row.names(heatmap),heatmap)
write.table(heatmap,file="heatmap.txt",sep="\t",row.names=F,quote=F)

#pheatmap图
install.packages("pheatmap")
setwd("D:/TCGA_methylation/pheatmap")
library(pheatmap)

rt <- read.table("heatmap.txt",sep="\t",header=T,row.names=1,check.names=F)
Type <- c(rep("normal",50),rep("tumor",503))
names(Type) <- colnames(rt)
Type <- as.data.frame(Type)

tiff(file="heatmap.tiff",width=45,height=70,units="cm",compression="lzw",bg="white",res=300)
pheatmap(rt,annotation=Type,color=colorRampPalette(c("green","black","red"))(50),cluster_cols=F,fontsize_row=5,fontsize_col=4)
dev.off()

#甲基化与转录组数据融合
setwd("D:/TCGA_methylation/MethyGeneMerge")
methy <- read.table(file="normalizeMethy.txt",row.names=1,header=T,sep="\t",check.names=F)
rna <- read.table(file="normalizeExp.txt",row.name=1,header=T,sep="\t",check.names=F)

colnames(methy) <- gsub("(.*?)-(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3-\\4",colnames(methy))
colnames(rna) <- gsub("(.*?)-(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3-\\4",colnames(rna))
sameSample <- intersect(colnames(methy),colnames(rna))
rownames(methy) <- paste(rownames(methy),"methy",sep="|")
rownames(rna) <- paste(rownames(rna),"exp",sep="|")
merge <- rbind(id=sameSample,methy[,sameSample],rna[,sameSample])
write.table(merge,file="merge.txt",sep="\t",quote=F,col.names=F)

#基因甲基化与转录组数据相关关系
setwd("D:/TCGA_methylation/methycor")
inputFile <- "merge.txt"
methyGene <- "NCAPD3|methy"
expGene <- "NCAPD3|exp"
rt <- read.table(inputFile,sep="\t",header=T,check.names=F,row.names=1)
i <- rt[methyGene,]
j <- rt[expGene,]
x <- as.numeric(i)
y <- log2(as.numeric(j)+1)
corT <- cor.test(x,y)
methyGeneName <- unlist(strsplit(methyGene,"\\|"))[[1]]
expGeneName <- unlist(strsplit(expGene,"\\|"))[[1]]
z <- lm(y~x)
estimate <- corT$estimate
cor <- round(estimate,3)
pvalue <- corT$p.value   #以小数的方式保留位数
pval <- signif(pvalue,4) #以科学计数法的方式保留位数
pval <- format(pval,scientific=T)

outTiff <- "cor.tiff"
tiff(file=outTiff,width=15,height=15,units="cm",compression="lzw",bg="white",res=300)
plot(x,y,type="p",pch=16,main=paste("Cor=",cor,"(p-value=",pval,")",sep=""),xlab=paste(methyGeneName,"methylation"),ylab=paste(expGeneName,"expression"))
lines(x,fitted(z),col=2)
dev.off()

#甲基化位点提取
nohup ./posget.pl metadata.cart.2018-08-01.json NCAPD3 > /home/train/Liuping/raw_data/TCGA_Methylation_ProstateCancer/NCAPD3_Methy_location.log 2>&1 &
#[1] 40635

#甲基化位点与gene表达合并
setwd("D:/TCGA_methylation/posmethygene")
methyFile="posMethy.txt"
expFile="normalizeExp.txt"
methy <- read.table(methyFile,row.names=1,header=T,sep="\t",check.names=F)
rna <- read.table(expFile,row.names=1,header=T,sep="\t",check.names=F)
colnames(methy) <- gsub("(.*?)-(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3-\\4",colnames(methy))
colnames(rna) <- gsub("(.*?)-(.*?)-(.*?)-(.*?)-.*","\\1-\\2-\\3-\\4",colnames(rna))
sameSample <- intersect(colnames(methy),colnames(rna))
rownames(methy) <- paste(rownames(methy),"methy",sep="|")
rownames(rna) <- paste(rownames(rna),"exp",sep="|")
merge <- rbind(id=sameSample,methy[,sameSample],rna[,sameSample])

#甲基化位点与gene表达相关
setwd("D:/TCGA_methylation/cor_posmethygene")
inputFile <- "merge.txt"
posmethy <- "cg14934141|methy"
expGene <- "NCAPD3|exp"
rt <- read.table(inputFile,sep="\t",header=T,check.names=F,row.names=1)
i <- rt[posmethy,]
j <- rt[expGene,]
x <- as.numeric(i)
y <- log2(as.numeric(j)+1)
corT <- cor.test(x,y)
posmethyid <- unlist(strsplit(posmethy,"\\|"))[[1]]
expgenename <- unlist(strsplit(expGene,"\\|"))[[1]]
z <- lm(y~x)
cor <- corT$estimate
cor <- round(cor,3)
pvalue <- corT$p.value
pval <- signif(pvalue,4)
pval <- format(pval,scientific=T)
outTiff <- "posCor.tiff"
tiff(file=outTiff,width=15,height=15,units="cm",compression="lzw",bg="white",res=300)
plot(x,y,type="p",pch=16,main=paste("Cor=",cor,"(p-value=",pval,")",sep=""),xlab=paste(posmethyid,"methylation"),ylab=paste(expgenename,"expression"))
lines(x,fitted(z),col=2)
dev.off()

合并的脚本,类似转录本的合并脚本,需要自己根据需要自行修改
欢迎添加下面这个群提问:

此条目发表在TCGA分类目录。将固定链接加入收藏夹。

TCGA甲基化数据分析》有108条回应

  1. 无聊说:

    请问如何把TCGA下载的文件汇总为一个矩阵

  2. 无聊说:

    你的outtab是不是缺少了几列:logFC、pvalue、normalMed、tumorMed、diffMed

  3. 无聊说:

    想请教一下TCGA下载的甲基化文件如何注释ID?我下载的450K文件每个ID对应多个基因名字,应如何配对?

  4. snowy说:

    甲基化位点提取矩阵,如何获得,谢谢!

  5. 说:

    请问merge之后的文件是第一列为基因名还是探针号,其他列为样本?然后数据都是β值的数据吗? 还是什么样的格式?我没有从TCGA上下载原始数据,是下载的整合后的,但是不知道该整理成什么样,才能作为limma的输入文件

  6. 求教说:

    亲,合并后的矩阵是什么样的格式呢,基因和样本信息,及β值吗

  7. 刘辉说:

    你好,能否帮我合并下甲基化数据,我的电脑合并不起来

    • daizao说:

      最近忙着毕业,不好意思啊

      • 刘辉说:

        没事,我后来合并好了,然后那个甲基化位点的提取,我提取后只有样本的ID,没有甲基化位点的,不知道能否发我个脚本,或者买也可以的,我的扣扣:1098210582,麻烦了!打扰了!感谢!

      • 刘辉说:

        我后来捣鼓了好久,合并好了,但甲基化位点没能提取成功,提取出来只有ID号,没有甲基化位点,不知道能否给我个脚本,或者帮我修改下脚本,有偿的,不好意思,打扰了!感谢!我的扣扣:1098210582

        • daizao说:

          因为实验室原因这个脚本目前还不能给你,而且忙于毕业课题所以木有时间帮你看脚本~
          抱歉啦~

          • 刘辉说:

            没事没事,如果后面有时间的话,方便的话能不能帮我合并下甲基化文件,因为我合并的时候老是有一列样品名全部为0

  8. 说:

    请问合并的矩阵作为输入文件是这个格式的吗?
    GENE SYMBOL NORMAL_1 NORMAL_2 TUMOR_1 TUMOR_2
    gene1 beta_value beta_value beta_value beta_value
    gene2 beta_value beta_value beta_value beta_value
    ….

    正常样本和癌症样本都有多个,往后排列?

  9. 小穆说:

    我在TCGA上下载的甲基化数据,我想确定promoter的甲基化水平。在下载的数据中,每一个cg号探针都对应着一个 position to tss 。按理说,每个探针都对应一个甲基化位点,所以position to tss应该是一个数字,但为什么下载的数据里面position to tss里是好多个数字呢

    • daizao说:

      因为一个基因可以有不同的甲基化位点,所以甲基化位点与TSS的距离,所以可以有多个

      • 小穆说:

        还是有点不明白。一个cg号探针应该是只对应一个甲基化位点吧,所以一个位点对应一个tss应该是一个数字啊。

        • daizao说:

          你说的没问题,但是一个基因对应着不同的转录本,而这些转录本的TSS并不一定是同一个,所以会出现多个数字(原始文件里面相同的gene_symbol有不同的Transcript_ID)
          一个cg号 -> 相同的gene_symbol -> 多个转录本(同一个基因) -> 多个TSS

  10. 谭增琦说:

    非常详细的甲基化数据分析,非常有用!请问你现在有空吗,可否帮忙下载个BRCA甲基化数据呢,我这边网实在是不行,如果可以方便的话,请加我QQ544931709

  11. 谭增琦说:

    非常细致用心的教程!受益匪浅,请问你现在还忙吗?可否麻烦你帮我下载个BRCA的甲基化数据呢?十分感谢!扣544931709

  12. 周瑶瑶说:

    你好,在运行代码
    for (i in row.names(data)){
    rt <- rbind(expression=data[i,],grade=grade)
    rt <- as.matrix(t(rt))
    wilcoxTest <- wilcox.test(expression ~ grade, data=rt)

    normalGeneMeans <- mean(data[i,1:normalNum])
    tumorGeneMeans <- mean(data[i,(normalNum+1):ncol(data)])
    logFC <- log2(tumorGeneMeans) – log2(normalGeneMeans)
    pvalue <- wilcoxTest$p.value
    normalMed <- median(data[i,1:normalNum])
    tumorMed <- median(data[i,(normalNum+1):ncol(data)])
    diffMed 0) & ( diffMed >0)) | (( logFC < 0 ) & ( diffMed <0 ))){
    outTab <- rbind(outTab,cbind(gene=i,normalMean=normalGeneMeans,TumorMean=tumorGeneMeans,logFC=logFC,pValue=pvalue))
    }}
    的时候出现Error in data[i,]:subscript out of bounds.但是不知道哪里错了求解决

  13. gongyue说:

    您好,打扰啦,想请问一下您的基因表达文件“normalizeExp”是怎么生成的呢

  14. 橘子说:

    甲基化差异分析流程按照您的代码跑还是有一些问题可以请教一下吗,不知道怎么联系您,我的QQ2325933960,有偿,麻烦啦

  15. 说:

    outDiff logFCfilter & as.numeric(as.vector(outTab$fdr)) < fdrFilter),]
    请问这里为什么需要挑选logFC的绝对值大于1的行呢?

  16. 宋晓帆说:

    您好,下载的数据里有NA值,这些值您是采取什么方法处理的呢?用什么算法填充一下吗?
    另外很感谢您的分享,对最近的学习有很大帮助。

  17. 新大麦味茶说:

    请问有原始TCGA甲基化数据整合成矩阵的代码,可以分享一下吗,谢谢

  18. le说:

    这是生信自学网的视频脚本吧,我好像用过

  19. xueshengxin说:

    请问怎样把TCGA下载下来的数据合并成矩阵呢?

  20. FLL说:

    您好,我想问一下在第一步甲基化差异分析中,为什么没有把分组信息输入呢?

  21. 马妮籽说:

    真是一个好博主,回答问题这么仔细,你毕业了吗

  22. shaoqing liu说:

    那么多位点,怎样确定一个基因的甲基化水平啊,是把各个位点的b值平均一下吗?

  23. 于永春说:

    不知道博主毕业了没有,有没有时间,我最近也是在合并矩阵上出的问题,文件太大了,电脑运行不下来,停止运行后得出的文件不知道全不全,样本数倒是对,郁闷,可不可以帮个忙,谢谢!

  24. fkxi说:

    ucsc xena上有合并号的矩阵可以直接下载,矩阵差异用什么做的? 查了minfi,ChAMP,输入格式 都不对。

  25. fkxi说:

    差异是怎么做的ttest吗?ucsc xena上有合并好的Illumina Human Methylation 450矩阵.
    查了minfi、ChAMP这些计算甲基化数据的包,格式都不对

  26. Jenson Zhou说:

    博主,你这里对提取后的甲基化数据用normalizeBetweenArrays做标准化看上去有点问题。normalizeBetweenArrays函数是用于标准化不同array之间数据的,它默认一列数据是一个array,而你的data结构中一列是一个样本的各个基因的甲基化β值,那么执行完normalizeBwtweenArrays函数后,不同样本之间的各基因甲基化水平的都被标准化到同样的scale了,而这显然是不合理的,因为偏恶性的肿瘤组织中,其各基因整体的甲基化水平应该比偏良性的肿瘤组织高。

  27. 刘杨说:

    有前列腺TCGA甲基化基因矩阵吗,求分享

  28. 彭振说:

    为什么有的甲基化位点会对应好几个基因啊?比如cg00000029对应的gene_symbol是RBL2;RBL2;RBL2,对于这种位点应该是删除还是做什么处理?

  29. 彭振说:

    为什么有的位点对应好几个gene_symbol?比如cg00000029对应的gene_symbol为RBL2;RBL2;RBL2,对于这类位点应该做什么处理呢?

    • daizao说:

      我取的第一个,因为有的甲基化位点在这个基因的附近或者在基因上,你可以用promoter区域过滤一下,然后再进行合并。

  30. 景程说:

    Heading ##楼主请问TCGA的450k的甲基化矩阵怎么过滤啊,没有idat文件

  31. Wheann说:

    Error in data[i, ] : subscript out of bounds
    您好,依旧是这个问题,请问可以帮忙看看是哪里出问题了吗?

  32. zhiqi Li说:

    您好,我看前面的一些评论包括您都有几次提到了promoter,但是我发现我下载的甲基化文件中并没有包含promoter的行,请问这是为什么?

发表评论

邮箱地址不会被公开。 必填项已用*标注

To create code blocks or other preformatted text, indent by four spaces:

    This will be displayed in a monospaced font. The first four 
    spaces will be stripped off, but all other whitespace
    will be preserved.
    
    Markdown is turned off in code blocks:
     [This is not a link](http://example.com)

To create not a block, but an inline code span, use backticks:

Here is some inline `code`.

For more help see http://daringfireball.net/projects/markdown/syntax

Protected with IP Blacklist CloudIP Blacklist Cloud