代码编织梦想

多形性成胶质细胞瘤(GBM)甲基化区域的计算鉴别

 

目的:找出胶质细胞瘤特异性甲基化区域,为临床诊断提供理论依据

 

步骤:

1、  查找数据:下载TCGA中GBM的RNA-seq和甲基化数据

2、  甲基化数据分析,正常肿瘤对比,进行差异甲基化分析,找出肿瘤样本中高甲基化区域

3、  对RNA-seq数据进行分析,正常肿瘤对比,差异表达基因的筛选,找出肿瘤样本中低表达基因。

4、  结合甲基化和RNA-seq数据,将高甲基化和低表达基因取交集,这些基因很可能属于抑癌基因,与抑癌基因取交集,再结合promoter区域的CpG整合分析,寻找候选靶标。

5、  对找出的靶标进行验证,利用pubmed以及其他数据库,反向验证靶标的可靠性

整体流程

整体流程

一、数据下载

首先进入TCGA下载数据GBM的RNA-seq和甲基化数据,从下表可见GBM共有172套RNA-seq数据以及437套DNA甲基化数据,由于TCGA提供Infinium HumanMethylation27 BeadChip和Infinium HumanMethylation450 BeadChip两种芯片平台的数据,为了避免后续不同芯片平台间数据合并的困难,仅下载HumanMethylation450的芯片数据,共计154套。

 

图表 1TCGA数据汇总

二、初步整理数据

使用TCGA-Assembler.2.0.5进行GBM数据批量下载与初步整理,并且绘制RNA-seq基因表达量盒型图以及甲基化芯片数据盒型图,由于数据量较大,此处不贴图。

三、整体可视化

首先对于甲基化数据,选取ID为TCGA.06.AABW.11A.31D.A368.05的数据,查看总体甲基化程度。由于每个位点真实情况只存在:甲基化/非甲基化两种,所以对全部位点甲基化程度进行统计,也应该是大部分位点处于“完全甲基化”(Methylation state=1)和“完全非甲基化”(Methylation state=0)两种状态,下图绘制了数据的频数柱状图,可以明显看出形状处于“两头高,中间低”,反向说明芯片数据质量较好。

 

图表 2单个样本CpG甲基化程度统计

接下来,对多个样本绘制CpG甲基化程度小提琴图,同一行是同一个病人,左边样本来源于Primary Solid Tumor,右边样本来源于Recurrent Solid Tumor,除了甲基化程度大部分分布于0和1附近外,还能看出来源于同一病人肿瘤的甲基化程度依旧会有略微差异。

TCGA barcode:https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode

 

图表 3小提琴图

 

同样的,对于RNA-seq数据也可以进行一些初步可视化,除了数据下载后绘制的盒型图,亦可以进行PCA初步查看数据分布,下图左为PCA陡坡图,反映了第一主成分、第二主成分…等等所拥有信息量的比例,下图右为使用PCA1和PCA2绘制的散点图,可以发现5个正常样本距离较近,从侧面反映数据可信度较好。

 

最后,对于RNA-seq表达谱数据,使用系统聚类方法,绘制树状图,可以发现5个正常样本距离也是很近,数据质量还行。

 

 

四、差异甲基化区域筛选

为了更加科学高效地筛选差异甲基化位点,参考bioconductor中甲基化芯片的分析流程,使用minfi包进行差异甲基化分析,得到差异甲基化位点。

http://www.bioconductor.org/help/workflows/methylationArrayAnalysis/

在检测的526733个CpG位点中,共有4927个CpG位点P值<0.01,且在肿瘤样本中保持着甲基化程度高于0.7,对应2054个基因。

五、差异表达基因筛选

由于数据源自RNA-seq,最主流的分析方法当然是基于负二项分布模型的DESeq2包。

先用MA-plot查看差异表达基因大致分布。意外的是,图形左侧有大概七条线状条纹,最初我怀疑这是sample之间有batch effect导致,需要用其他更好normalize的方法,后来用identify方法挨个找出每条线上的基因名及其对应的表达量,发现这些基因在172套样本中表达量几乎全为0,仅有一两个样本有一点点表达,这种数据的存在导致这些线状条纹的产生。

 

图表 4 MA-plot

然后, 选取p值最小的差异表达基因,绘制其在不同组间表达量,确实差异很显著。

 

图表 5表达量散点图

接着,绘制差异表达基因在不同组间的表达量热图,正常样本是图片最左边的五列,当然如果需要解释具体的生物学问题,需要将聚类出来的每一类,将差异表达基因进行GO以及KEGG注释,结合有关的生物学表型,探讨其分子机制及意义

 

图表 6表达量热图

最后,选取筛选条件为p值小于0.01且log2FoldChange<-2的差异表达基因,在肿瘤样本低表达的基因,共计1657个

六、抑癌基因的获取

在pubmed中查询研究人员整理的tumor suppressor genes,果然在Nucleic Acids Res发表过TSGene数据库,存储了抑癌基因的列表。https://bioinfo.uth.edu/TSGene/download.cgi

下载全部1217个人类抑癌基因的列表。

All the 1217 human tumor suppressor genes with basic annotations

 

图表 7 TSGene database

七、数据整合

对于甲基化数据中,肿瘤样本高甲基化CpG附近的基因,RNA-seq中肿瘤样本低表达的基因,以及TSGene数据库中下载的抑癌基因列表,三者做overlap,找出特异性的候选靶标,为后续分析做准备,下图为三者overlap的韦恩图。

 

图表 8数据整合韦恩图

共计找出12个候选靶标基因。

八、靶标筛选

之前筛选选择的单个CpG的差异甲基化,而实际临床检测应用时候,可能需要多个CpG作为对照,因此统计了12个候选靶标基因TSS前1.5kb内所有CpG的甲基化程度,然后绘制热图,可以明显发现,虽然当初用CpG的差异甲基化位点筛出来的基因都是肿瘤样本高甲基化的,可是统计TSS前1.5kb内所有CpG的甲基化程度,这些基因却有很多在所有样本中都是低甲基化状态,而看上去很靠谱的是NUAK1基因,其正常样本在TSS前1.5kb内低甲基化,肿瘤样本中对应区域高甲基化。

 

图表 9 methylation heatmap

NUAK1基因TSS前1.5kb内共检测了7个CpG,这7个CpG在154个样本中检测出来的甲基化程度如下图,可以明显看出来这7个CpG在Tumor组织中甲基化程度都相对高,而在Normal组织中甲基化程度相对较低。

 

图表 10 NUAK1的TSS区CpG甲基化程度

使用Cpgplot预测CpG island位置,这7个CpG在promoter5’到3’序列图上相对位置如下

1035 1094 1408 1413 1444 1448 1471

 

图表 11 CpG island预测

参数使用:

NUAK1promoter from 1 to 1500

     Observed/Expected ratio > 0.40

     Percent C + Percent G > 40.00

     Length > 100

CpG island详细信息: Length 101 (1086..1186)   Length 105 (1366..1470)

这七个CpG基本都在CpGisland中,具体序列见附录

 

九、靶标基因相关讨论

进入Gene数据库搜索NUAK1相关内容,可以发现基因全称NUAK family kinase 1,还是个激酶,激酶的话就对调控会有很大作用了,而在HPA RNA-seq normal tissues项目中,又看出来这个激酶在脑中表达量明显高于其他组织,这又与发生在脑部的GBM不谋而合。

 

图表 12 NUAK1相关讨论

 

十、分子机制探讨

对于肿瘤组织中高甲基化CpG附近的,并且在肿瘤样本中低表达的intersect共计274个基因,使用Gene Ontology进行富集分析,可以明显发现在GO biological process生物学过程中的“神经系统发育”、“化学性突触传递”和“细胞膜的组织”等部分里面有着富集,特别是“中枢神经系统的髓鞘形成”,富集程度达到26.95倍,这又与研究的多发生于脑补的GBM有着密切的联系,反向验证实验结果的正确性。

 

图表 13 GO富集分析

十一、FurtherMore

根据生物学知识可以得到,CpG的甲基化会调控基因的转录,因此,Transcript Start Site(TSS)附近的甲基化程度值得进行一番深入研究,选用人类基因组hg19版本,对23056基因共计46489个转录起始位点,进行转录起始位点富集甲基化程度统计。

统计TSS前后5000bp内CpG甲基化程度,并且使用曲线进行拟合,可以发现TSS处的CpG Methylation水平明显降低,这也与科学常识相吻合。

 

图表 14 TSS附近甲基化程度

附录

NUAK1promoter区CpG island用蓝色标注,检测的CpG用加粗横线下标标注。

>NUAK1promoter

TATGAAAGGAGAAGGGGGAGCTTTGGAACTGGACAGGTAGGGTTTAAATCCTGGTCCTGC

CATTTACAAACTGTGTAACCTCTGGGAAATTACTTAACCTTTCTGATACGGTTTCTTCAA

TTGAAAATAGGGATTGTAAACAGCTACTTTACAGAGGAGGGTTTACTGTCATAAAACAGT

ACCAGCCTATGGTAGATGATGCTGTTGATTCAATAGATACTGATGAAGTCCCACATATCT

GGGAGTAACACTATCAGCCAGAATAAGCCAGGTTATGCTGCTTTAACAAATCCCACTGAC

TTAACATAAATAAAGTTTAATATTTGCTGACACTACTTTTCCAATGCAGATTGGCTGGGA

TTTTCTGCCATCTTTACACAGAGACTCAGGCTGGCTTGGGGAGGCTCCAACTTTGTACCA

CCATGATTGTCAAGATAGGAAAAAGAAGACATGGGGATTTGCTCACTGGCTCCTAAAGGT

TCCAAGTGAAAGCAACACAATCACTTCTGCTCTCATTCCATTGGCCTAAGCAAGTCCGAT

GGTCTCATCTAACTTCAAGAGGATGGGGTAGTGGATTTCTACCACATCGCAGGGGAGGGG

AAGAAGTAGAAAATATTTTAAAATACTACTGTAGACACAATGTGTTTATCTCTACAATAG

ATCTGCTAAATCCATATCTTAAGAGAAAACCGAAGCTCTGAGAGGTTATGTCATTTGACC

AAGGCCACACAGCCAATCCTCTGGGAAGCCAGGACTCAACCTGCATGCCACTCATCCTGA

CTGTGGGATCTATAGGCACAACGTCCACAAACTGTATAGAAGCAGAACAGTTTCAGGATG

GGGGTGGGTGGCAGGGAGGCCCTGCAAATGCGGTGGACAGAATAGGAATTGAATAGGCAT

GATGCGCCTTTGTGTCTGTGTTCTGTCATTGCTCCTGGGGAAGGGAAGAAGGGGCAGGGA

AGTTTGTGTGAAATATCAGTAGAAGGAAAGTTTGGACAGGTAAGAATATCACGCGTCAGA

GTACAGCAGAGACACGTGTGGAGGATGAGGGCAGTGTTTCAGGCCATTACTCTGGCAGCA

GTGAGGAGGCTCCCGGGGAGGGTTGGGGAGAGCGGGCTGTTGCTGGAGTCCGGGGTAAGG

TGACCAGGGGTAGACAGGAATGAAGGCGGCGGGAATGTAGAGGAAGGGCTGCACCTGAGA

GCTGCAGAGGAGGCCCCAGTGAGGTCCACTGGTGGAAAAGAACCGCCATGCGATCTGCAG

AACACCAGACCTTCTTCAGCCTTCACTCTTCCCACCCTAGTCTGGTACATTGCACCACTT

GTTAAAAAAAAAAATTTCCCTAAGACCTCTTTTTCTCTAGCCTTCTTCCTTATTTTTCAT

GGTCTCTTCTTAGAACGGCGGCAGCCACGCCGCGGTGGGAGGCCCTTCCTGCCTGACCCT

TACCGTGCGGGGTACCGTTCCTGTCACCATCGCCAGGATCTGGCCCTTTCAGTGAAAGGA

 

diffFind.R

 

setwd("G:/AllShare/SkillTrainHomework/BasicDataProcessingResult")
load("GBM__methylation_450.rda")
setwd("G:/AllShare/SkillTrainHomework")
DATAFRAME<- data.frame(Data)

TCGASampleName<- colnames(DATAFRAME)
TCGASamplelist <- strsplit(TCGASampleName,split=".",fixed = TRUE)
TCGASampletpye <- c()
TCGASampleID <- c()
for(i in 1:length(TCGASamplelist)){
  TCGASampletpye <- c(TCGASampletpye,TCGASamplelist[[i]][4])
  TCGASampleID <- c(TCGASampleID,paste(TCGASamplelist[[i]][2],TCGASamplelist[[i]][3],sep = "."))
}
TCGASampletpyeResult <- c()
for(i in 1:length(TCGASampletpye)){
  if(as.integer(substr(TCGASampletpye[i],1,2))<10){
    TCGASampletpyeResult <- c(TCGASampletpyeResult,"Tumor")
  }else{
    TCGASampletpyeResult <- c(TCGASampletpyeResult,"Normal")
  }
}

library(ggplot2)
#h<- DATAFRAME$TCGA.06.0125.01A.01D.A45W.05
h<- DATAFRAME$TCGA.06.AABW.11A.31D.A368.05
pdf("CpGmethlation.pdf")
ggplot(NULL,aes(x=h))+
  geom_histogram(binwidth = 0.02)+
  labs(x = "% methylation per base",title = "Histogram of % CpG methylation\nSampleID:TCGA.06.AABW.11A.31D.A368.05")
dev.off()


library(reshape2)
tmp <- Data[,1:8]
tmp <- data.frame(tmp)
colnames(tmp) <- TCGASampleID[1:8]
tmp2 <- stack(tmp)
tmp2 <- data.frame(tmp2)
colnames(tmp2) <- c("CpG methylation","SampleID")
pdf("MethlationViolin.pdf",width =12,height = 10 )
ggplot(tmp2,aes(x=SampleID,y=`CpG methylation`,fill=SampleID))+
  geom_violin()+
  facet_wrap(~SampleID,ncol=2)
dev.off()


library(minfi)
dmp <- dmpFinder(Data, pheno = TCGASampletpyeResult, type = "categorical")
save(dmp,file = "diff.Rdata")
load("diff.Rdata")

meanTumorData<- apply(Data[,TCGASampletpyeResult=="Tumor"],1,mean)

#&!is.na(meanTumorData)

# head(meanTumorData)
# head(dmp)
# head(signifidmp,10)
# summary(dmp)
# rownames(dmp)
dmp$genelistID <- as.integer(rownames(dmp))
o<- order(dmp[,"genelistID"])
dmp<- dmp[o,]
#筛选出p<0.01且无空值的CpG,并且正常样本甲基化程度<0.3,即筛选肿瘤中高甲基化的基因
signifidmp <- dmp[meanTumorData>0.7 & !is.na(meanTumorData) & dmp$pval<0.01
                  &!is.na(dmp$pval) & !is.na(dmp$qval) & !is.na(dmp$intercept),]
#signifidmp <- signifidmp[signifidmp$intercept<0.3,]dmp <- 
signifiData<- Data[as.integer(rownames(signifidmp)),]
signifiDes<- Des[as.integer(rownames(signifidmp)),]

# Data[Des[,1]=="ch.3.2438620R",]
# meanTumorData[Des[,1]=="ch.3.2438620R"]>0.7
# mean(Data[103429,TCGASampletpyeResult=="Normal"])
# mean(Data[103429,TCGASampletpyeResult=="Tumor"])
# summary(Data[Des[,1]=="ch.3.2438620R",])

# apply((head(signifiData,10)[,TCGASampletpyeResult=="Tumor"]),1,mean)

signifiGene <- signifiDes[,2]

length(unique(signifiGene))

head(signifidmp)

#intercept正常样本甲基化程度
# Data[362831,118]
# Data[362831,48]
# Data[362831,]
# 
# Data[310369,118]
# Data[310369,48]
# Data[310369,]



#以下是关于CpG的一些计算


#因为Des有部分空缺,取出非空部分,生成Desfull
Desfull<- Des[!is.na(Des[,3]),]
#同样取出Data里非空部分,计算mean
methtstat<- apply(Data[!is.na(Des[,3]),],1,mean,na.rm=TRUE)
#去掉计算mean值以后NA的部分对应的Desfull
Desfull<- Desfull[!is.na(methtstat),]
#再去掉methtstat的mean值中NA的部分
methtstatresult<- methtstat[!is.na(methtstat)]
grCpG <- GRanges(seqnames = paste("chr",Desfull[,3],sep = ""),
              ranges = IRanges(start = as.integer(Desfull[,4]), width = 1))
grCpG$value <- methtstatresult


library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb_hg19 <- TxDb.Hsapiens.UCSC.hg19.knownGene
trans <- as.data.frame(transcripts(txdb_hg19))
trans<- trans[trans$seqnames %in% c("chrX","chrY",paste("chr",1:22,sep="")),]
grTSS <- GRanges(seqnames = trans$seqnames,
                 ranges = IRanges(start = trans$start-5000,end = trans$start+5000))

grCpG
grTSS
hitObj<- findOverlaps(grTSS,grCpG)
CpGRelativeSite <- c()
CpGRelativeMeth <- c()
for(i in 1:length(hitObj)){
  #取出对应CpG真实位置
  CpGsite <- grCpG[hitObj[i]@to]@ranges@start
  #计算回TSS真实位置
  TSSsite<- mean(grTSS[hitObj[i]@from]@ranges)
  #计算CpG相对位置并且存储
  CpGRelativeSite <- c(CpGRelativeSite,TSSsite - CpGsite)
  #取出CpG平均甲基化程度
  CpGRelativeMeth <- c(CpGRelativeMeth,grCpG[hitObj[i]@to]$value)
}
save.image(file = "diffresult.Rdata")

load("diffresult.Rdata")

CpGResult <- c()
for(i in -5000:5000){
  CpGResult <- c(CpGResult,mean(CpGframe[CpGframe$CpGRelativeSite==i,"CpGRelativeMeth"]))
}
CpGframe <- data.frame(-5000:5000,CpGResult)
colnames(CpGframe) <- c("CpGRelativeSite","MethState")
ggplot(CpGframe, aes(x=CpGRelativeSite, y=MethState))+
  #geom_point()+
  geom_smooth()
ggsave(filename="TSS附近甲基化程度.pdf",width = 12,height=8)





library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb_hg19 <- TxDb.Hsapiens.UCSC.hg19.knownGene
trans <- as.data.frame(transcripts(txdb_hg19))
grTSS <- GRanges(seqnames = trans$seqnames,
                 ranges = IRanges(start = trans$start-5000,end = trans$start+5000))


genesall<- genes(txdb_hg19)
SERPIND1<- genesall[genesall$gene_id=="3053"]
SERPIND1promoter<- promoters(SERPIND1,upstream = 1500,downstream = 0)
SERPIND1promoter
SERPIND1hit<- findOverlaps(SERPIND1promoter,grCpG)
SERPIND1hit@to
grCpG[SERPIND1hit@to]
#将Data里对应数据取出
Datafull<- Data[!is.na(Des[,3]),]
Datafull<- Datafull[!is.na(methtstat),]

NUAK1<- genesall[genesall$gene_id=="9891"]
NUAK1promoter<- promoters(NUAK1,upstream = 1500,downstream = 0)
NUAK1promoter

NUAK1hit<- findOverlaps(NUAK1promoter,grCpG)
NUAK1hit@to
grCpG[NUAK1hit@to]
#将Data里对应数据取出
Datafull<- Data[!is.na(Des[,3]),]
Datafull<- Datafull[!is.na(methtstat),]


apply(Datafull[NUAK1hit@to,TCGASampletpyeResult=="Normal"],1,mean)
apply(Datafull[NUAK1hit@to,TCGASampletpyeResult=="Tumor"],1,mean)
library(ggplot2)
normalNUAK1<- data.frame(stack(Datafull[NUAK1hit@to,TCGASampletpyeResult=="Normal"]))
tumorNUAK1 <- data.frame(stack(Datafull[NUAK1hit@to,TCGASampletpyeResult=="Tumor"]))
normalNUAK1$sampletype <- "Normal"
normalNUAK1$site <- 1:7
tumorNUAK1$sampletype <- "Tumor"
tumorNUAK1$site <- 1:7
frame1<- data.frame(normalNUAK1$NA..1,normalNUAK1$sampletype,normalNUAK1$site)
colnames(frame1) <- c("MethState","sampletype","site")
frame2<- data.frame(tumorNUAK1$NA..1,tumorNUAK1$sampletype,tumorNUAK1$site)
colnames(frame2) <- c("MethState","sampletype","site")
resultframe <- rbind(frame1,frame2)
ggplot(resultframe, aes(x=site, y=MethState,colour=sampletype)) + 
  scale_x_discrete(limits=1:7)+
  geom_point(position="jitter")
ggsave(filename="TargetCpGNUAK1.pdf",width = 8,height=6)


apply(Datafull[SERPIND1hit@to,TCGASampletpyeResult=="Normal"],1,mean)
apply(Datafull[SERPIND1hit@to,TCGASampletpyeResult=="Tumor"],1,mean)
library(ggplot2)
normalSERPIND1<- data.frame(stack(Datafull[SERPIND1hit@to,TCGASampletpyeResult=="Normal"]))
tumorSERPIND1 <- data.frame(stack(Datafull[SERPIND1hit@to,TCGASampletpyeResult=="Tumor"]))
normalSERPIND1$sampletype <- "Normal"
normalSERPIND1$site <- 1:2
tumorSERPIND1$sampletype <- "Tumor"
tumorSERPIND1$site <- 1:2
frame1<- data.frame(normalSERPIND1$NA..1,normalSERPIND1$sampletype,normalSERPIND1$site)
colnames(frame1) <- c("MethState","sampletype","site")
frame2<- data.frame(tumorSERPIND1$NA..1,tumorSERPIND1$sampletype,tumorSERPIND1$site)
colnames(frame2) <- c("MethState","sampletype","site")
resultframe <- rbind(frame1,frame2)
ggplot(resultframe, aes(x=site, y=MethState,colour=sampletype)) + 
  scale_x_discrete(limits=1:2)+
  geom_point(position="jitter")
ggsave(filename="TargetCpGSERPIND1.pdf",width = 8,height=6)

library("BSgenome.Hsapiens.UCSC.hg19")
library(seqinr)
genome <- BSgenome.Hsapiens.UCSC.hg19
grCpG[NUAK1hit@to]
promoterSeq<- getSeq(genome,NUAK1promoter)
write.fasta(promoterSeq$`9891`,"NUAK1promoter","promoterSeq.txt")


normalNUAK1$NA..1
tumorNUAK1$NA..1

dim(Datafull)
dim(Desfull)
# 尝试用mice包补缺失值,由于数据量太大而取消
# library(mice)
# library(VIM)
# #md.pattern(Data)
# tmp <- Data[,1:3]
# aggr_plot <- aggr(tmp, col = c('navyblue', 'red'), numbers=TRUE, sortVars=TRUE, 
#                   labels=names(tmp),cex.axis=.7, gap=2, 
#                   ylab=c("Histogram of missing data", "Pattern"))
# tempData <- mice(Data,m=1,maxit=50,meth='pmm',seed=500)
# 
# 尝试T检验,由于有缺失值而取消
# Ttest<- t.test(Data[1,1:3],Data[1,c(48,118)])
# 
# Data[1,118]
# Data[1,48]
# tmp<- Data[1,]
# myfun <-function(x){
#   Ttest<- t.test(x[c(1:47,49:117,119:154)],x[c(48,118)])
#   return(Ttest$p.value)
# }
# myfun(Data[1,])

# result<- apply(Data,1,myfun)
# group <- factor(TCGASampletpyeResult,levels=c("Normal","Tumor"))
# design <- model.matrix(~-1+group)
# design
# fit.reduced <- lmFit(Data,design)
# fit.reduced <- eBayes(fit.reduced)
# summary(decideTests(fit.reduced))
# 
# 尝试使用missMethyl包,最后designMatrix设置有错误而差异区域识别错误,差异不明显
# top<-topTable(fit.reduced,coef=1)
# top
# cpgs <- as.integer(rownames(top))
# Data[cpgs[10],]
# Data[cpgs[10],48]
# Data[cpgs[10],118]
# 
# par(mfrow=c(2,2))
# for(i in 1:4){
#   stripchart(Data[rownames(Data)==cpgs[i],]~design[,4],method="jitter",
#              group.names=c("Normal","Tumor"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
#              vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
#   title(cpgs[i],cex.main=1.5)
# }


Downloaded.R

 

 

 

 

setwd("G:/AllShare/SkillTrainHomework/TCGA-Assembler.2.0.5/TCGA-Assembler")
#' Load functions
source("Module_A.R")
source("Module_B.R")
setwd("G:/AllShare/SkillTrainHomework")
#' choose a cancer type
#' 可查看网址https://tcga-data.nci.nih.gov/docs/publications/tcga/
sCancer <- "GBM"
sPath1 <- "./DownloadedData"
sPath2 <- "./BasicDataProcessingResult"
sPath3 <- "./AdvancedDataProcessingResult"
#下载RNA-Seq数据
path_geneExp <- DownloadRNASeqData(cancerType = sCancer,
                        assayPlatform = "gene.normalized_RNAseq",
                        saveFolderName = sPath1)
#' Download DNA methylation 450 data
#' 使用Illumina的甲基化分析芯片测出来的甲基化数据
path_methylation_450 <- DownloadMethylationData(cancerType = sCancer,
                        assayPlatform = "methylation_450",
                        saveFolderName = sPath1)
# 下载出来的格式说明:
# TCGA-3C-AAAU-01A-11R-A41B-07
# 前三位TCGA-3C-AAAU表示病人ID
# 第四位01A表示肿瘤类型Tumor types range from 01 - 09,
# normal types from 10 - 19 and control samples from 20 - 29.
# 第五位Order of portion in a sequence of 100 - 120 mg sample portions
# 可查看网址https://wiki.nci.nih.gov/display/TCGA/TCGA+barcode

#' Process gene expression data
#' 对RNA-seq下载结果进行处理,将基因名进行处理
list_geneExp <-
  ProcessRNASeqData(inputFilePath = path_geneExp[1],
                    outputFileName = paste(sCancer,
                                           "geneExp",
                                           sep = "__"),
                    dataType = "geneExp",
                    outputFileFolder = sPath2)
#处理methylation数据
list_methylation_450 <-
  ProcessMethylation450Data(inputFilePath =
                    path_methylation_450[1], outputFileName = paste(sCancer,
                    "methylation_450", sep = "__"), outputFileFolder = sPath2)

#'Perform advanced data processing using Module B functions
#'Advanced进一处理methylation数据
list_methylation_450_OverallAverage <-
  CalculateSingleValueMethylationData(input = list_methylation_450,
                    regionOption = c("TSS1500", "TSS200"), DHSOption = "Both",
                    outputFileName = paste(sCancer, "methylation_450", sep = "__"),
                    outputFileFolder = sPath3,
                    chipAnnotationFile = "./SupportingFiles/MethylationChipAnnotation.rda")
save.image("tmpdata.Rdata")
load("tmpdata.Rdata")

 

 

RNAdiff.R

 

#setwd("G:/AllShare/SkillTrainHomework/BasicDataProcessingResult")
#load("GBM__geneExp.rda")
setwd("G:/AllShare/SkillTrainHomework")
load("RNAdiff.RData")
DATAFRAME<- data.frame(Data)
TCGASampleName<- colnames(DATAFRAME)
TCGASamplelist <- strsplit(TCGASampleName,split=".",fixed = TRUE)
TCGASampletpye <- c()
TCGASampleID <- c()
for(i in 1:length(TCGASamplelist)){
  TCGASampletpye <- c(TCGASampletpye,TCGASamplelist[[i]][4])
  TCGASampleID <- c(TCGASampleID,paste(TCGASamplelist[[i]][2],TCGASamplelist[[i]][3],sep = "."))
}
TCGASampletpyeResult <- c()
for(i in 1:length(TCGASampletpye)){
  if(as.integer(substr(TCGASampletpye[i],1,2))<10){
    TCGASampletpyeResult <- c(TCGASampletpyeResult,"Tumor")
  }else{
    TCGASampletpyeResult <- c(TCGASampletpyeResult,"Normal")
  }
}
phenotype <- data.frame(TCGASampletpyeResult)
rownames(phenotype) <- colnames(Data)
rownames(Data) <- Des[,2]
library(DESeq2)
Data<- floor(Data)
dds <- DESeqDataSetFromMatrix(countData = Data,
                              colData = phenotype,
                              design = ~ TCGASampletpyeResult)
dds <- DESeq(dds)
res <- results(dds)
resOrdered <- res[order(res$padj),]
resSig <- subset(resOrdered, padj < 0.01)
res001 <- results(dds, alpha=0.01)
pdf(file = "MA-plot.pdf",width = 15,height = 9)
plotMA(res001, ylim=c(-3,3))
idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]
res[idx,]
Data[c("100131561","677811","414899","677846"),]
Data[c("9271","140893"),]
Data[c("390077","144742"),]
dev.off()
d <- plotCounts(dds, gene=which.min(res$padj), intgroup="TCGASampletpyeResult", 
                returnData=TRUE)

library("ggplot2")
ggplot(d, aes(x=TCGASampletpyeResult, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))+
  labs(title=paste("EntrezID:",res@rownames[which.min(res$padj)]))+
  theme(plot.title = element_text(hjust = 0.5))
ggsave(filename="p值最小的基因不同表达量.pdf",width = 8,height=6)

ntd <- normTransform(dds)
library("pheatmap")
select <- rownames(resOrdered)[1:20]
df <- as.data.frame(colData(dds)[,"TCGASampletpyeResult"])
pdf(file = "热图.pdf",width = 20,height = 9)
pheatmap(assay(ntd)[select,],show_rownames=TRUE)
dev.off()

RnaSigGeneID<- rownames(resSig)
RnaSigGeneName<- Des[Des[,2]%in%RnaSigGeneID,1]

intersect(RnaSigGeneName,signifiGene)

#log2FoldChange为负数,正常高表达;log2FoldChange为正数,肿瘤高表达
#筛选log2FoldChange为负数,即肿瘤低表达的部分
resSigdown<- resSig[resSig$log2FoldChange<(-2),]
#dim(resSigdown)
RnaSigGeneIDdown<- rownames(resSigdown)
RnaSigGeneNamedown<- Des[Des[,2]%in%RnaSigGeneIDdown,1]

# Data["2498",c(38:41,76)]
# Data["2498",]
# Data["7153",c(38:41,76)]
# Data["7153",]

#与抑癌基因取交集
Human_TSGs <- read.table("Human_TSGs.txt",header = TRUE,sep = "\t")
#TSGene <- read.table("TSGene-LOFdataset.txt",header = TRUE,sep = "\t")
#sig_exp <- read.table("sig_exp.txt",header = TRUE,sep = "\t")
TSGs <- Human_TSGs$GeneSymbol
#TSGs <- TSGene$GeneName
#TSGs <- sig_exp$Symbol
result1<- intersect(RnaSigGeneNamedown,TSGs)
result2<- intersect(result1,signifiGene)#三者交集
resultmuch <- intersect(RnaSigGeneNamedown,signifiGene)
write.table(resultmuch,"候选靶标基因多多多.txt",quote = FALSE,row.names = FALSE,col.name = FALSE)

#"CCHCR1"%in%signifiGene




#此段代码需要methylation的TCGASampletpyeResult对象
setwd("G:/AllShare/SkillTrainHomework/AdvancedDataProcessingResult")
load("GBM__methylation_450__TSS1500-TSS200__Both.rda")
MethTssFrame<- data.frame(Data)
rownames(MethTssFrame) <- Des[,1]
targetFrame<- MethTssFrame[result2,]
colnames(targetFrame) <- TCGASampletpyeResult
targetFrame<-na.omit(targetFrame)
library(pheatmap)
pdf(file = "target热图.pdf",width = 20,height = 6)
pheatmap(targetFrame,show_rownames=TRUE)
dev.off()

library (VennDiagram)
draw.triple.venn(area1=5, area2=5, area3=5
                 ,n12=3, n23=3, n13=3, n123=3
                 ,category = c('A','B','C'))
pdf("交集.pdf",width = 10,height = 10)
T<-venn.diagram(list(A=RnaSigGeneNamedown,B=TSGs,C=signifiGene),filename=NULL
                ,lwd=1,lty=2,category = c('RNA-seq down','Tumor suppressor gene','Hypermethylation'))
grid.draw(T)
dev.off()


#save.image("RNAdiff.RData")


PCAdata<-t(Data)
#作主成分分析
PCAdata.pr<-prcomp(PCAdata,scale=FALSE)
#作预测
PCA_eset<- predict(PCAdata.pr)
colnames(PCA_eset)
pdf(file = "陡坡图.pdf",width = 14,height = 7)
screeplot(PCAdata.pr,type="lines")
dev.off()

data.hc <- hclust( dist(PCAdata))
pdf(file = "树状图.pdf",width = 22,height = 12)
plot(data.hc, hang = -1)
dev.off()

#plot(PCA_eset[,1:2])

library("ggplot2")
ggplot(NULL, aes(x=PCA_eset[,1], y=PCA_eset[,2], colour=TCGASampletpyeResult)) + 
  geom_point() + 
  guides(color=guide_legend(title=NULL)) +
  labs(x = "PCA1 34.7%",y = "PCA2 15.3%",title = "RNA-seq Principal components analysis")
ggsave(filename="PCA RNA-seq.pdf",width = 8,height=6)


GOterm <- read.table("analysis.txt",header = TRUE,sep = "\t")
library(ggplot2)
ggplot(GOterm,aes(x=GO.biological.process,y=upload_1..over.under.,fill = P.value))+
  geom_bar(stat="identity")+
  labs(x = "GO terms",y = "Fold Enrichment",title = "GO biological process analysis")+
  coord_flip()
ggsave(filename="GO analysis.pdf",width = 8,height=6)

 

 

 

 

 

 

 

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/qq_29300341/article/details/78816189

一文读懂:dna甲基化的作用及各种高通量检测方法比较-爱代码爱编程

DNA甲基化相关定义 DNA甲基化一直以来都是表观遗传学领域研究的重点之一。DNA甲基化(Methylation)是指在DNA甲基转移酶(DNA methyltransferase, 缩写DNMT)的作用下,基因组DNA序

DNA甲基化在重头甲基转移酶远古丢失后数百万年的进化持久性-爱代码爱编程

胞嘧啶甲基化是DNA中广泛存在的一种修饰,发挥重要的作用。在酵母新型隐球菌(Cryptococcus neoformans),CG甲基化发生在富含转座子的重复序列中,而且需要DNA甲基转移酶Dnmt5。科学家发现Dnmt5在体内和体外都表现出精确的甲基化维持特异性,并且在体内显示出与后生动物维持甲基化的Dnmt1酶利用相似的辅酶因子。值得注意的是,系统发育

精准医学: 尿液DNA甲基化检测有助于膀胱癌早期和复发监测|早期筛查专题-爱代码爱编程

易点评 本研究是膀胱癌预防、检测方面医学研究的典范。作者巧妙运用了中山纪念医院(SYSMH)、癌症基因组图谱(TCGA)和基因表达总表(GEO)的研究数据,以及MassARRAY检测技术,确定了膀胱癌特异性甲基化标记物,并针对其开发出的utMeMA监测方法进行了鉴定,从而证实了对肿瘤患者尿液DNA甲基化进行检测可以诊断出早期、微小以及残留的肿瘤这一结果,与

DeepRMethylSite:一种基于深度学习的蛋白质精氨酸甲基化位点预测方法-爱代码爱编程

DeepRMethylSite:一种基于深度学习的蛋白质精氨酸甲基化位点预测方法 https://www.researchgate.net/publication/341890599_DeepRMethylSite_A_Deep_Learning_based_approach_for_Prediction_of_Arginine_Methylati

学科前沿:基因启动子甲基化与宫颈癌发展的关系 | 文献科普-爱代码爱编程

宫颈癌作为四大常见癌症和世界上第二大最常见的妇科恶性肿瘤,据环球观察称,2018年全球约有57万例新发病例,超过31.1万万例死亡。宫颈癌起源于正常的上皮下包,有明显的形态学改变,包括两种主要的组织类型:鳞状细胞癌(SCC)和腺癌(AC),分别占85-90%和10-25%。 宫颈癌流行病学和分子研究表明,人乳头瘤病毒(HPV)的持续存在是宫颈癌发生的主要原

DNA甲基化数据分析专题-爱代码爱编程

欢迎关注”生信修炼手册”! DNA 甲基化作为重要的表观遗传学的标记,能够调控基因的表达,在生长发育和疾病相关研究领域都有着重要意义。测定甲基化的手段有很多,芯片作为一种成熟的手段,其稳定性,可重复性以及性价比,使得在DNA甲基化研究领域芯片占据了半壁江山。 对于human来说,目前主流的DNA甲基化芯片有450K 和 850K 两种,都是illumi

易基因 | 文献速递:重亚硫酸盐扩增子测序研究通过DNA甲基化监测急性髓系白血病MRD-爱代码爱编程

大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。 关于DNA甲基化测序方法,可以点击查看:常用的6种DNA甲基化测序方法,你知道几个? 今天要分享解读的是一篇发表于Leukemia(IF 11.528)的DNA甲基化研究文章,本文通过重亚硫酸盐扩增子测序,分析了通过研究DNA甲基化模式来监测急性髓系白血病(acute myeloid le

易基因 | 文献速递:RRBS方法绘制1538例乳腺癌甲基化图谱并预测癌症发生/预后-爱代码爱编程

大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。 错过RRBS技术在人和小鼠疾病表观遗传特征研究的可点:Mol Biol Evol | 利用RRBS技术多维度分析人和小鼠的疾病表观遗传特征 今天一起来看看发表于Nat Commun期刊的利用RRBS测序分析技术绘制了1538例乳腺癌组织样本和244个相邻正常组织样本的DNA甲基化图谱,并揭

综述科普:单细胞测序技术下的小鼠脑部DNA甲基化图谱-爱代码爱编程

表观基因组的动态变化与哺乳动物大脑中的细胞分化和成熟有关,在调节神经元功能和动物行为方面具有重要作用。胞嘧啶DNA甲基化(5mC)是一种稳定的共价修饰,对基因调控至关重要。在哺乳动物基因组中,5mC主要发生在CpG位点(mCG),在具有组织和细胞类型特异性的调控元件上表现出动态模式,调节转录因子的结合亲和力,控制基因转录。非CpG胞嘧啶在小鼠和人脑中也有丰

易基因|精准医学: TERT启动子DNA甲基化在癌症中的双重作用-爱代码爱编程

DNA甲基化是基因表达的重要调控机制,加拿大多伦多大学的Donghyun D. Lee团队从生物学和临床角度讨论TERT启动子DNA甲基化在肿瘤中双重作用的最新进展,相关结果发表在2021年11月的J Clin Invest(14.808)期刊上。 缩略词: TERT:端粒酶逆转录酶(telomerase reverse transcriptase

易基因:cfDNA甲基化分析揭示造血细胞移植的所有主要并发症|前沿技术-爱代码爱编程

大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。 异体造血细胞移植(HCT)为血液恶性肿瘤和免疫疾病提供了有效的治疗方案,然而频繁的移植并发症限制了HCT后的长期获益。本研究发现循环游离DNA(cfDNA)是一种可用于监测HCT后主要并发症的多功能分析物,通过cfDNA重亚硫酸盐测序分析(cfDNA-WGBS)和疾病特异性生物信息学分析

易基因|基于cfDNA甲基化的肿瘤液体活检如何开展?-爱代码爱编程

大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。 液体活检是指通过微创或无创的方式收集血液或其他体液(包括尿液、腹水、胸腔积液等),并对其中肿瘤相关的物质进行分析的一种检验方式。以血液为例,除了正常的红细胞、白细胞、来自正常细胞的游离 DNA (cell free DNA, cfDNA),还可以检测到一些来自肿瘤的物质,包括循环肿瘤细胞 (

NmRF:从RNA序列中鉴定多物种RNA2‘-o-甲基化修饰位点(假尿苷位点)-爱代码爱编程

文章目录 前言一、简介二、方法和材料2.1 数据收集和准备2.2 序列特征向量的构造2.3 核苷酸二元剖面2.4 二核苷酸2.5 核苷酸化学性质2.6 特征选择2.7 机器学习方法2.8 评价指标三、结果3.1 不同分类器的性能比较3.2 特征选择方法的比较3.3 特征分析和跨物种性能评价3.4 比较现有的最先进的方法3.5 Web服务器实现总结

易基因|综述:单细胞DNA甲基化分析方法全介绍及未来发展前景预测-爱代码爱编程

​大家好,这是专注表观组学十余年,领跑多组学科研服务的易基因。 2021年07月10日,《Biomolecules》杂志上发表一篇关于单细胞表观测序的综述文章,详细介绍了单细胞DNA甲基化的实验策略、分析方法、数据分析以及相关科研应用并讨论了单细胞甲基化测序的未来发展前景。以下为原文总结分享: Introduction 在这篇综述中,作者

易基因|干货:m6a rna甲基化merip-seq测序分析实验全流程解析_易基因科技的博客-爱代码爱编程

易基因|干货:m6A RNA甲基化MeRIP-seq测序分析实验全流程解析 大家好,这是专注表观组学十余年,领跑多组学科研服务的易基因。 本期,我们讲讲甲基化RNA免疫共沉淀(MeRIP-seq/m6A-seq)实验怎么做,从技术原理、建库测序流程、信息分析流程和研究套路等四方面详细介绍。 一、甲基化RNA免疫共沉淀(MeRIP-seq/

喜报!易基因“同源基因特异性甲基化时序数据分析方法”获专利授权_易基因科技的博客-爱代码爱编程

大家好,这里是专注表观组学十余年,领跑多组学科研服务的易基因。 喜报!深圳市易基因科技有限公司发明的《一种同源基因特异性甲基化时序数据的分析方法》获专利授权。该方法首先利用DNA甲基化测序数据序列特征和突变信息判定ASM区域,筛选出每个时期或者组别中特有的所述ASM区域,之后结合不同突变型对应的甲基化信息,利用某个时期的甲基化信息和其他不同时期样品进

易基因|dna甲基化方法全解析:方法发展、技术应用、优缺点_技术应用的优缺点-爱代码爱编程

大家好,这是专注表观组学十余年,领跑多组学科研服务的易基因。 2021年03月,《Methods》杂志以“DNA methylation methods: global DNA methylation and methylomic analyses”为题发表了关于DNA甲基化分析方法的综述文章,详细介绍了DNA甲基化分析方法的发展变化、DNA甲基化分析方