• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • (未测试)TCGA数据库的癌症甲基化芯片数据重分析

      (未测试)TCGA数据库的癌症甲基化芯片数据重分析

      • 发布者 weinfoauthor
      • 分类 未分类
      • 日期 2020年3月4日
      • 评论 0评论

      最近生信技能树更新了一系列甲基化数据分析教程,本周我们一起来完成其中的一个作业:TCGA数据库的各个癌症甲基化芯片数据重新分析

      作业来自文献:Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma,下载TCGA数据库中口腔癌的甲基化芯片信号值矩阵,然后挑选有N-T配对的32个病人的数据进行差异分析,画火山图和热图。

      准备数据

      Step1. 从UCSC xena(https://xenabrowser.net/datapages/)中下载头颈癌HNSC的甲基化数据HumanMethylation450.gz(数据很大,如果下载失败,那应该是网络的问题),表型数据HNSC_clinicalMatrix.gz,生存数据HNSC_survival.txt.gz。

      Step2. 挑选有Normal-Tumor配对的32个OSCC病人

      # UCSC的xena浏览器下载指定HNSC的phenotype数据
      pd.all <- read.delim("./raw_data/HNSC_clinicalMatrix", header = T, stringsAsFactors = F)
      colnames(pd.all)
      pd <- pd.all[,c("sampleID","anatomic_neoplasm_subdivision","sample_type")]
      table(pd$anatomic_neoplasm_subdivision)
      # 挑出OSCC样本(参考文章methods)
      oscc <- c("Oral Cavity","Oral Tongue","Buccal Mucosa","Lip","Alveolar Ridge","Hard Palate","Floor of mouth")
      pd <- pd[pd$anatomic_neoplasm_subdivision %in% oscc,]
      # 挑出normal-tumor配对样本
      pd$patient <- substr(pd$sampleID,1,12)
      pd <- pd[pd$patient %in% pd$patient[pd$sample_type=="Solid Tissue Normal"],]
      rownames(pd) <- pd$sampleID
      pd$sample_type <- ifelse(pd$sample_type=="Primary Tumor","Tumor","Normal")
      # 保存配对样本(此时有48对)
      write.table(pd$sampleID, file = "./raw_data/samples.txt",quote = F, row.names = F, col.names = F)

      Step3.缩减甲基化信号矩阵:因为整个HNSC甲基化矩阵有580个样本,读进R会很慢,而且没有必要,所以先用linux命令根据sample.txt筛选出N-T配对的OSCC样本。

      # 解压甲基化信号矩阵文件
      $ gunzip HumanMethylation450.gz
      
      # 文件内容是这样的
      $ head -5 samples.txt 
      TCGA-CV-5436-01
      TCGA-CV-5436-11
      TCGA-CV-5442-01
      TCGA-CV-5442-11
      TCGA-CV-5966-01
      $ cat HumanMethylation450 | cut -f 1-6 |head -5
      sample    TCGA-CN-6022-01 TCGA-CQ-5327-01 TCGA-CV-5435-11 TCGA-HD-7753-01 TCGA-CQ-5325-01
      cg13332474    0.0220  0.0264  0.0307  0.0263  0.0218
      cg00651829    0.0209  0.4514  0.0214  0.0231  0.0248
      cg17027195    0.0443  0.2414  0.0330  0.0361  0.0416
      cg09868354    0.0490  0.0546  0.0469  0.0754  0.0516
      
      ## 筛选N-T配对的OSCC样本
      # 1.获得OSCC样本所在的列
      $ cols=($(sed '1!d;s/\t/\n/g' HumanMethylation450 | grep -nf samples.txt | sed 's/:.*$//'))
      # 2.筛选出OSCC样本的甲基化信号矩阵
      $ cut -f 1$(printf ",%s" "${cols[@]}") HumanMethylation450 > OSCC_NT_methy450k.txt

      载入数据为champ对象

      ## 安装ChAMP包
      if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
      BiocManager::install("ChAMP") #有点漫长,可能会报错,把报错的包手动安装一下应该没什么问题
      library("ChAMP")
      
      ## 读取具有N-T配对的OSCC甲基化信号矩阵
      library(data.table)
      a <- fread("./raw_data/OSCC_NT_methy450k.txt", data.table = F )
      a[1:4,1:4]
      rownames(a)=a[,1]
      a=a[,-1]
      beta=as.matrix(a)
      # beta信号值矩阵里面不能有NA值
      beta=impute.knn(beta) 
      betaData=beta$data
      betaData=betaData+0.00001
      sum(is.na(betaData))
      
      #表型数据
      pd <- pd[colnames(betaData),]
      table(pd$patient) #此时有些normal样本没有配对的tumor
      NT.s <- pd$patient[duplicated(pd$patient)] #得到32对样本,和文章一致了
      pd <- pd[pd$patient %in% NT.s,]
      betaData <- betaData[,pd$sampleID]
      dim(betaData)
      
      # 载入为ChAMP对象
      library(ChAMP)
      myLoad=champ.filter(beta = betaData ,pd = pd) #这一步已经自动完成了过滤
      dim(myLoad$beta)
      save(myLoad,file = './Rdata/step1-output.Rdata')

      甲基化矩阵质控

      rm(list = ls())
      load('./Rdata/step1-output.Rdata')
      # 数据归一化
      myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=5)
      dim(myNorm)
      QC.GUI(beta=myNorm,arraytype="450K") # 显示有NA值
      num.na <- apply(myNorm,2,function(x)(sum(is.na(x))))
      hist(num.na) # 看来是有样本有问题,个别样本校正后NA值过高
      table(num.na) # 61个样本没有NA,3个样本含有超过250k个NA
      myNorm <- myNorm[,which(num.na < 250000)] # 删掉异常样本
      save(myNorm,file = './Rdata/step2-champ_myNorm.Rdata')
      
      # 主成分分析
      library("FactoMineR")
      library("factoextra") 
      dat <- t(myNorm)
      
      pd <- myLoad$pd[colnames(myNorm),] #去掉异常样本
      group_list=pd$sample_type
      table(group_list)
      
      dat.pca <- PCA(dat, graph = FALSE) 
      fviz_pca_ind(dat.pca,
                   geom.ind = "point", 
                   col.ind = group_list, 
                   addEllipses = TRUE, 
                   legend.title = "Groups")
      
      # 热图
      cg=names(tail(sort(apply(myNorm,1,sd)),1000))
      library(pheatmap)
      ac=data.frame(group=group_list)
      rownames(ac)=colnames(myNorm)  
      pheatmap(myNorm[cg,],show_colnames =F,show_rownames = F,
               annotation_col=ac,filename = 'heatmap_top1000_sd.png')
      dev.off()
      
      # 相关关系矩阵热图
      # 组内的样本的相似性应该是要高于组间的!
      pheatmap::pheatmap(cor(myNorm[cg,]),
                         annotation_col = ac,
                         show_rownames = F,
                         show_colnames = F)
      
      # 去除异常样本
      # 图中看出TCGA-CV-5971-01,TCGA-CV-6953-11,TCGA-CV-6955-11这3个样本有点异常
      myNorm <- myNorm[,!(colnames(myNorm) %in% c("TCGA-CV-5971-01","TCGA-CV-6953-11","TCGA-CV-6955-11"))]
      pd <- pd[colnames(myNorm),]
      save(pd,myNorm,file = "./Rdata/filtered.Rdata")

      从主成分分析图可以看出,其实甲基化信号并没有将tumor和normal分得很开,top1000sd的信号信号矩阵和相关关系矩阵热图显示有3个样本看着很奇怪,我把它们删除了。

      差异分析

      rm(list=ls())
      load("./Rdata/filtered.Rdata")
      # 差异分析
      group_list <- pd$sample_type
      myDMP <- champ.DMP(beta = myNorm,pheno=group_list)
      head(myDMP[[1]])
      save(myDMP,file = './Rdata/step3-output-myDMP.Rdata')

      画图

      根据差异分析结果画火山图,因为作者的差异分析用的是paired t test,和我的结果差异有点大,我自己选了差异甲基化的cutoff。

      df_DMP <- myDMP$Tumor_to_Normal[,1:5]
      logFC_cutoff <- 0.45
      df_DMP$change <- ifelse(df_DMP$adj.P.Val < 10^-15 & abs(df_DMP$logFC) > logFC_cutoff,
                              ifelse(df_DMP$logFC > logFC_cutoff ,’UP’,’DOWN’),’NOT’) 
      table(df_DMP$change) 

      this_tile <- paste0(‘Cutoff for logFC is ‘,round(logFC_cutoff,3),
                          ‘\nThe number of up gene is ‘,nrow(df_DMP[df_DMP$change ==’UP’,]) ,
                          ‘\nThe number of down gene is ‘,nrow(df_DMP[df_DMP$change ==’DOWN’,]))

      library(ggplot2)
      g <- ggplot(data=df_DMP, 
                  aes(x=logFC, y=-log10(adj.P.Val), 
                      color=change)) +
        geom_point(alpha=0.4, size=1) +
        theme_set(theme_set(theme_bw(base_size=10)))+
        xlab(“log2 fold change”) + ylab(“-log10 p-value”) +
        ggtitle( this_tile ) + theme(plot.title = element_text(size=10,hjust = 0.5))+
        scale_colour_manual(values = c(‘blue’,’black’,’red’))
      print(g)

      作者用了Cox回归从差异分析结果中获得了影响生存的CpGs,画热图

      ### 挑选影响生存的CpG位点,画热图
      # 提取差异甲基化CpGs
      choose_gene <- rownames(df_DMP[df_DMP$change != "NOT",])
      choose_matrix <- myNorm[choose_gene,]
      
      ## 批量生存分析
      # 载入生存信息
      phe <- read.delim("./raw_data/HNSC_survival.txt.gz",header = T, stringsAsFactors = F)
      phe <- phe[phe$sample %in% colnames(choose_matrix),]
      phe <- phe[substr(phe$sample,14,15)=="01",]
      choose_matrix <- choose_matrix[,phe$sample]
      
      library(survival)
      logrankP <- apply(choose_matrix, 1, function(x){
        #x <- choose_matrix[1,]
        phe$group <- ifelse(x>mean(x),"High","Low")
        res <- coxph(Surv(OS.time, OS)~group, data=phe)
        beta <- coef(res)
        se <- sqrt(diag(vcov(res)))
        p.val <- 1 - pchisq((beta/se)^2, 1)
      })
      names(logrankP) <- rownames(choose_matrix)
      table(logrankP<0.05) #151个CpG位点
      surv_gene <- names(sort(logrankP))[1:20] #挑了20个生存最显著的DMP,用于画图
      
      # 热图
      plot_matrix <- myNorm[surv_gene,]
      annotation_col <- data.frame(Sample=myLoad$pd$sample_type ) #热图数据的bar
      rownames(annotation_col) <- colnames(plot_matrix)
      ann_colors = list(Sample = c(Normal="green", Tumor="red"))
      library(pheatmap)
      pheatmap(plot_matrix,show_colnames = F,annotation_col = annotation_col,
               border_color=NA,
               color = colorRampPalette(colors = c("white","navy"))(50),
               annotation_colors = ann_colors)

      本次作业主要是学习用ChAMP包处理TCGA甲基化数据的流程,用的方法和原文不同,结果差异有点大。原文筛选配对样本的目的是为了进行“配对样本t检验”,我们这里其实可以保留全部OSCC肿瘤样本的。

      总结:

      1. TCGA的甲基化数据一般都很大,可以的话,在读入R之前先进行一下样本筛选
      2. 差异分析之前一定要进行质控检验,本例中就发现了异常样本
      3. 更完整的甲基化芯片信号值矩阵差异分析代码见Jimmy老师的github:https://github.com/jmzeng1314/methy_array/

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      weinfoauthor
      weinfoauthor

      1233

      上一篇文章

      (未测试)3+肿瘤EMT纯生信,热点分析来袭
      2020年3月4日

      下一篇文章

      (未测试)一文看懂主成分分析
      2020年3月4日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?