七步走纯R代码通过数据挖掘复现一篇实验文章(第1到6步)

文章里面是自己测了8个TNBC病人的转录组然后分析,但是我们有TCGA数据库,所以可以复现,这就是为什么标题是七步走纯R代码通过数据挖掘复现一篇实验文章!

主要流程:

  • 下载数据
  • 数据清洗
  • 质量控制
  • 差异分析
  • 注释mRNA,lncRNA
  • 富集分析
  • WGCNA(因为排版限制,内容见本期第3条推文

Step1. 数据下载

这里参考去年的学徒数据挖掘:送你一篇TCGA数据挖掘文章

  • 表达矩阵下载及临床信息下载
  • GTF文件

Step2. 数据清洗


提高数据清洗的能力,将会很大程度的提高你做分析数据的速度,可能有的人还是习惯用Excel来清洗数据,但是我建议能用代码的尽量用代码解决,数据清洗思路也很重要,一定要清楚你的目标,然后思考可能实现的途径。


# 挑选三阴性乳腺癌的样本
# 
## FALSE  TRUE 
##  1165   118
tnbc_sample = phenotype_file[phenotype_file$breast_carcinoma_estrogen_receptor_status == 'Negative' & 
                    phenotype_file$breast_carcinoma_progesterone_receptor_status == 'Negative' & 
                    phenotype_file$lab_proc_her2_neu_immunohistochemistry_receptor_status == 'Negative', ]

save(tnbc_sample,file = '../02_data/tnbc_sample.Rdata')
  • 然后在118个三阴性乳腺癌中找到有N-T配对样本的病人

在TCGA中第14,15位的数字01~09代表肿瘤样本,10以上则为正常样本


# 把肿瘤样本提取出来,把正常样本提取出来然后根据前十二字符merged到的样本就是属于配对样本
library(stringr)
b = tnbc_sample[1:2]

tumor_sample <- b[substr(b$submitter_id.samples,14,15) < 10,]
tumor_sample$TCGA_ID <- str_sub(tumor_sample$submitter_id.samples, 1, 12)

normal_sample <- b[!substr(b$submitter_id.samples,14,15) < 10,]
normal_sample$TCGA_ID <- str_sub(normal_sample$submitter_id.samples, 1, 12)

tmp = merge(normal_sample, tumor_sample, by = "TCGA_ID")
# 以下是为了方便后续提取数据
a = tmp[,2:3]
colnames(a)
## [1] "submitter_id.samples.x"             
## [2] "additional_pharmaceutical_therapy.x"
names(a) = c("submitter_id.samples", "additional_pharmaceutical_therapy")

b = tmp[,4:5]
names(b) = c("submitter_id.samples", "additional_pharmaceutical_therapy")

TNBC_pair_sample = rbind(a,b)
head(TNBC_pair_sample)
save(TNBC_pair_sample, file = "../02_data/TNBC_pair_sample.Rdata")
  • 在配对样本中过滤掉并非同时有正常和肿瘤组织测序的样本
rm(list = ls())
library(readr)
# 提取表达矩阵
load(file = "../02_data/TNBC_pair_sample.Rdata")
rawdata <- read_tsv("../02_data/TCGA-BRCA.htseq_counts.tsv")
rawdata[1:3,1:3]
## # A tibble: 3 x 3
##   Ensembl_ID         `TCGA-3C-AAAU-01A` `TCGA-3C-AALI-01A`
##   <chr>                           <dbl>              <dbl>
## 1 ENSG00000000003.13               9.35               8.71
## 2 ENSG00000000005.5                1.58               1.58
## 3 ENSG00000000419.11              10.9               10.8
rawdata = as.data.frame(rawdata)
rownames(rawdata) = rawdata[,1]
rawdata = rawdata[, -1]
rawdata[1:3,1:3]
##                    TCGA-3C-AAAU-01A TCGA-3C-AALI-01A TCGA-3C-AALJ-01A
## ENSG00000000003.13         9.348728         8.714246         10.35645
## ENSG00000000005.5          1.584963         1.584963          5.72792
## ENSG00000000419.11        10.874981        10.834471         10.32980
table(colnames(rawdata) %in% TNBC_pair_sample$submitter_id.samples)
## 
## FALSE  TRUE 
##  1191    26
expr = rawdata[,colnames(rawdata) %in% TNBC_pair_sample$submitter_id.samples]
expr = as.data.frame(t(expr))
expr[1:3,1:3]
##                  ENSG00000000003.13 ENSG00000000005.5 ENSG00000000419.11
## TCGA-A7-A4SE-01A           10.89482          4.906891          10.903882
## TCGA-AC-A2BK-01A           13.06811          0.000000          11.758640
## TCGA-AC-A2QJ-01A           10.27496          4.643856           9.763212
# 还要继续过滤掉不符合要求的样本

b = expr[1:2]
b$submitter_id.samples = rownames(b)
# 为了去掉其中一个样本测了两次tumor
tumor_sample <- b[substr(b$submitter_id.samples,14,16) == "01A",]
tumor_sample$TCGA_ID <- str_sub(tumor_sample$submitter_id.samples, 1, 12)

normal_sample <- b[!substr(b$submitter_id.samples,14,15) < 10,]
normal_sample$TCGA_ID <- str_sub(normal_sample$submitter_id.samples, 1, 12)

tmp = merge(normal_sample, tumor_sample, by = "TCGA_ID")
###-----------------------------------------------------------------------------
a = tmp[,2:4]
colnames(a)
## [1] "ENSG00000000003.13.x"   "ENSG00000000005.5.x"   
## [3] "submitter_id.samples.x"
names(a) = c("ENSG00000000003.13", "ENSG00000000005.5", "submitter_id.samples")

b = tmp[,5:7]
names(b) = c("ENSG00000000003.13", "ENSG00000000005.5", "submitter_id.samples")

TNBC_pair_sample = rbind(a,b)
table(table(str_sub(TNBC_pair_sample$submitter_id.samples,1,12)))
## 
##  2 
## 10
# ---------------------------------------------------------------------------------
table(rownames(expr) %in% TNBC_pair_sample$submitter_id.samples)
## 
## FALSE  TRUE 
##     6    20
expr = expr[rownames(expr) %in% TNBC_pair_sample$submitter_id.samples,]
dim(expr)
## [1]    20 60483
save(expr,file = "../02_data/TNBC_pair_expr.Rdata")

Step3. 质量控制


表达矩阵质量控制还有很多其他的图和方式,例如PCA等等,需要大家了解


rm(list = ls())
load(file = "../02_data/TNBC_pair_expr.Rdata")       plot(hclust(dist(expr)))

img

rownames(expr) # 学徒在这里考虑到模仿的文章里面是8个病人,所以这里剔除了2个病人。# 考虑的是hclust的离群点
exprSet = expr[c(-5,-6,-13,-14),]
save(exprSet,file = "../02_data/TCGA_TNBC_QC_exprSet.Rdata")

step4. 差异分析

  • 差异分析可选edger,limma,DESeq2,这里使用DESeq2
rm(list = ls())
library(stringr)
library(DESeq2)
library(DT)
load(file = "../02_data/TCGA_TNBC_QC_exprSet.Rdata")
raw_data = exprSet
raw_data <- as.data.frame( t(raw_data) )
raw_data[1:5, 1:6] 
##                    TCGA-BH-A0B3-01A TCGA-BH-A0B3-11B TCGA-BH-A18V-01A
## ENSG00000000003.13        11.287135        12.728346        10.100662
## ENSG00000000005.5          3.584963         8.864186         4.643856
raw_data <- 2^raw_data - 1
raw_data <- ceiling( raw_data )
raw_data[1:5, 1:6]
##                    TCGA-BH-A0B3-01A TCGA-BH-A0B3-11B TCGA-BH-A18V-01A
## ENSG00000000003.13             2499             6785             1097
## ENSG00000000005.5                11              465               24
pick_row <- apply( raw_data, 1, function(x){
  sum(x == 0) < 16
})
raw_data <- raw_data[pick_row, ]
dim(raw_data  )
## [1] 48236    16
# 使用DESeq2需要表达矩阵和group_list
a = as.data.frame(rownames(exprSet))
names(a) = "sample_id"
group_list=factor(ifelse(as.numeric(substr(colnames(raw_data),14,15)) < 10,'tumor','normal'))
table(group_list)
## group_list
## normal  tumor 
##      8      8
colData <- data.frame(row.names=colnames(raw_data),
                       group_list=group_list)

## 第一列有名称时,tidy=TRUE
raw_data[1:3,1:3]
##                    TCGA-BH-A0B3-01A TCGA-BH-A0B3-11B TCGA-BH-A18V-01A
## ENSG00000000003.13             2499             6785             1097
## ENSG00000000005.5                11              465               24
## ENSG00000000419.11             1853             1921             3564
head(colData)
##                  group_list
## TCGA-BH-A0B3-01A      tumor
## TCGA-BH-A0B3-11B     normal
dds <-DESeqDataSetFromMatrix(countData=raw_data, 
                             colData=colData, 
                             design=~group_list,
                             tidy=FALSE)

dds <- DESeq(dds)
vsd <- vst(dds, blind = FALSE)
# 制作表达矩阵小心结果相反,这里是tumor vs. normal
contrast <- c("group_list","tumor","normal")
dd1 <- results(dds, contrast=contrast, alpha = 0.05)
plotMA(dd1, ylim=c(-2,2))

img

# lfcShrink矫正
dd2 <- lfcShrink(dds, contrast=contrast, res=dd1)
plotMA(dd2, ylim=c(-2,2))

img

save(res,vsd, file = "../02_data/TCGA_TNBC_different_analysis.Rdata")

step5. 基因注释


lncRNA注释还是有挺多方法的,gencode上的文件同样可以可以来注释,至于哪些类型的注释是属于lncRNA需要自己查看相关的资料


# 加载包
# ------------------------------------------------------------------------------- 
rm(list = ls())
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
library(rtracklayer)
library(dplyr)
library(tidyr)
library(pheatmap)

# 构建基因注释的GTF文件
# --------------------------------------------------------------------------------
gtf1 <- rtracklayer::import('../02_data/Homo_sapiens.GRCh38.97.chr.gtf')
gtf_df <- as.data.frame(gtf1)
colnames(gtf_df)
##  [1] "seqnames"                 "start"                   
##  [3] "end"                      "width"                   
##  [5] "strand"                   "source"                  
##  [7] "type"                     "score"                   
##  [9] "phase"                    "gene_id"                 
## [11] "gene_version"             "gene_name"               
## [13] "gene_source"              "gene_biotype" 
gtf = gtf_df[,c(10,14)]
head(gtf)
##           gene_id                       gene_biotype
## 1 ENSG00000223972 transcribed_unprocessed_pseudogene
## 2 ENSG00000223972 transcribed_unprocessed_pseudogenesave(gtf,file = "../02_data/Homo_sapiens.GRCh38.97.chr.Rdata")

# 注释mRNA并找出显著差异的mRNA# ----------------------------------------------------------------------------------
load(file = "../02_data/TCGA_TNBC_different_analysis.Rdata")
res_nopoint <- res %>% 
  tidyr::separate(gene_id,into = c("gene_id"),sep="\\.") 
res_nopoint[1:3,1:3]##           gene_id  baseMean log2FoldChange## 1 ENSG00000000003 3910.1490   -0.008576929## 2 ENSG00000000005  428.1032   -5.759778527## 3 ENSG00000000419 2000.6141    0.618374851mRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #筛选gene,和编码指标
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>% 
  dplyr::inner_join(res_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = "|")
dim(mRNA_exprSet)## [1] 18949     7
diffSig_mRNA <- mRNA_exprSet[with(mRNA_exprSet, (abs(log2FoldChange)>2 & padj < 0.05 )), ]
dim(diffSig_mRNA)## [1] 2308    7# 注释lncRNA并找出显著差异的lncRNA# ----------------------------------------------------------------------------------
table(gtf_df$gene_biotype)LncRNA_exprSet <- gtf_df %>% 

  dplyr::filter(type=="gene",gene_biotype %in% "lncRNA") %>% #注意这里是transcript_biotype

  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>% 

  dplyr::distinct() %>% #删除多余行

  dplyr::inner_join(res_nopoint,by ="gene_id") %>% 

  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = "|")dim(LncRNA_exprSet)## [1] 12065     7
diffSig_lncRNA <- LncRNA_exprSet[with(LncRNA_exprSet, (abs(log2FoldChange)>2 & padj < 0.05 )), ]
dim(diffSig_lncRNA)## [1] 1127    7# 提取显著差异的mRNA,lncRNA表达矩阵,并进行ID转换# -----------------------------------------------------------------------
vst = as.data.frame(assay(vsd))vst$gene_id = rownames(vst)
vst_nopoint <- vst %>% 
  tidyr::separate(gene_id,into = c("gene_id"),sep="\\.") 
vst_nopoint[1:3,1:3]rownames(vst_nopoint) = vst_nopoint$gene_id
vst_nopoint[1:3,1:3]vst_nopoint = vst_nopoint[, -ncol(vst_nopoint)]
vst_nopoint[1:3,1:3]# 构建lncRNA_ids#-----------------------------------------------------------------------------------------
a = as.data.frame(diffSig_lncRNA$gene_id %>%
                    str_split("\\|",simplify = T))
names(a) = c("symbol", "probe_id","gene_type")
# diffSig_lncRNA ID转换#-----------------------------------------------------------------------------------------

ids = a[,1:2]
exprSet_lncRNA = vst_nopoint
exprSet_lncRNA[1:3,1:3] class(ids$symbol)## [1] "character"
table(rownames(exprSet_lncRNA) %in% ids$probe_id)## ## FALSE  TRUE ## 47197  1039
exprSet_lncRNA = exprSet_lncRNA[rownames(exprSet_lncRNA) %in% ids$probe_id,]
dim(exprSet_lncRNA)## [1] 1039   16
ids = ids[match(rownames(exprSet_lncRNA), ids$probe_id),]
tmp = by(exprSet_lncRNA,
         ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x)
         )])
probes = as.character(tmp)
dim(tmp)## [1] 1039
exprSet_lncRNA = exprSet_lncRNA[rownames(exprSet_lncRNA) %in% probes, ]
dim(exprSet_lncRNA)## [1] 1039   16
rownames(exprSet_lncRNA)=ids[match(rownames(exprSet_lncRNA),ids$probe_id),1]
exprSet_lncRNA[1:3,1:3]# 构建mRNA_ids#----------------------------------------------------------------------------------------
a = as.data.frame(diffSig_mRNA$gene_id %>%
                    str_split("\\|",simplify = T))
names(a) = c("symbol", "probe_id","gene_type")
# diffSig_mRNA ID转换# --------------------------------------------------------------------------------------
ids = a[,1:2]
exprSet_mRNA = vst_nopointexprSet_mRNA = exprSet_mRNA[rownames(exprSet_mRNA) %in% ids$probe_id,]
dim(exprSet_mRNA)## [1] 2273   16
ids = ids[match(rownames(exprSet_mRNA), ids$probe_id),]tmp = by(exprSet_mRNA,
         ids$symbol,
         function(x) rownames(x)[which.max(rowMeans(x)
         )])
probes = as.character(tmp)
dim(tmp)## [1] 2272
exprSet_mRNA = exprSet_mRNA[rownames(exprSet_mRNA) %in% probes, ]
exprSet_mRNA[1:3,1:3]##                 TCGA-BH-A0B3-01A TCGA-BH-A0B3-11B TCGA-BH-A18V-01A## ENSG00000000005         4.443206         8.357049         5.006992## ENSG00000003249        11.148636         8.672442        12.093123## ENSG00000003436         8.780552        12.019202        10.205594
dim(exprSet_mRNA)## [1] 2272   16
rownames(exprSet_mRNA)=ids[match(rownames(exprSet_mRNA),ids$probe_id),1]
exprSet_mRNA[1:3,1:3]##        TCGA-BH-A0B3-01A TCGA-BH-A0B3-11B TCGA-BH-A18V-01A## TNMD           4.443206         8.357049         5.006992## DBNDD1        11.148636         8.672442        12.093123## TFPI           8.780552        12.019202        10.205594
  • 画个显著差异mRNA的热图和火山图
 # diffSig_mRNA热图
  #--------------------------------------------------------------------------------
  group_list=factor(ifelse(as.numeric(substr(colnames(exprSet_lncRNA),14,15)) < 10,'tumor','normal'))
  table(group_list)
  (colData <- data.frame(row.names=colnames(exprSet_lncRNA),
                         group_list=group_list))
  pheatmap(exprSet_lncRNA,scale = "row",show_rownames = FALSE, show_colnames = FALSE,
           annotation_col = colData)

img

  ## mRNA火山图
  # -----------------------------------------------------------------------------------


  logFC_cutoff = 2
  DEG = mRNA_exprSet
  DEG = na.omit(DEG)

  DEG$change = as.factor( ifelse( DEG$padj< 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
                                  ifelse( DEG$log2FoldChange > logFC_cutoff , 'UP', 'DOWN' ), 'NOT' ) )
  table(DEG$change)
## 
##  DOWN   NOT    UP 
##  1228 15548  1045
  this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                      '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                      '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
  )
  library(ggplot2)
  g = ggplot(data=DEG, 
             aes(x=log2FoldChange, y=-log10(padj), 
                 color=change)) +
    geom_point(alpha=0.4, size=1.75) +
    theme_set(theme_set(theme_bw(base_size=20)))+
    xlab("log2 fold change") + ylab("-log10 p-value") +
    ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
    scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
  print(g)

img

  • lncRNA同样可以画出相应的图
 # diffSig_lncRNA热图
  # ----------------------------------------------------------------------------------------
  group_list=factor(ifelse(as.numeric(substr(colnames(exprSet_mRNA),14,15)) < 10,'tumor','normal'))
  table(group_list)
## group_list
## normal  tumor 
##      8      8
  (colData <- data.frame(row.names=colnames(exprSet_mRNA),
                         group_list=group_list))
##                  group_list
## TCGA-BH-A0B3-01A      tumor
## TCGA-BH-A0B3-11B     normal
## TCGA-BH-A18V-01A      tumor
## TCGA-BH-A18V-11A     normal
## TCGA-BH-A1F6-01A      tumor
## TCGA-BH-A1F6-11B     normal
## TCGA-BH-A1FC-01A      tumor
## TCGA-BH-A1FC-11A     normal
## TCGA-E2-A158-01A      tumor
## TCGA-E2-A158-11A     normal
## TCGA-E2-A1LH-01A      tumor
## TCGA-E2-A1LH-11A     normal
## TCGA-E2-A1LS-01A      tumor
## TCGA-E2-A1LS-11A     normal
## TCGA-GI-A2C9-01A      tumor
## TCGA-GI-A2C9-11A     normal
  pheatmap(exprSet_mRNA,scale = "row",show_rownames = FALSE, show_colnames = FALSE,
           annotation_col = colData)

img

  # lncRNA 火山图
  # -------------------------------------------------------------------------------------
  logFC_cutoff = 2
  DEG = LncRNA_exprSet
  DEG = na.omit(DEG)

  DEG$change = as.factor( ifelse( DEG$padj< 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
                                  ifelse( DEG$log2FoldChange > logFC_cutoff , 'UP', 'DOWN' ), 'NOT' ) )
  table(DEG$change)
## 
## DOWN  NOT   UP 
##  611 7067  428
  this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                      '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                      '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
  )
  library(ggplot2)

  g = ggplot(data=DEG, 
             aes(x=log2FoldChange, y=-log10(padj), 
                 color=change)) +
    geom_point(alpha=0.4, size=1.75) +
    theme_set(theme_set(theme_bw(base_size=20)))+
    xlab("log2 fold change") + ylab("-log10 p-value") +
    ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
    scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
  print(g)

img

  save(mRNA_exprSet,LncRNA_exprSet,diffSig_lncRNA, diffSig_mRNA, exprSet_mRNA,exprSet_lncRNA,
       file = "../02_data/TCGA_TNBC_pair_annotation_result.Rdata")

step6.富集分析

rm(list = ls())
library("clusterProfiler")
library("org.Hs.eg.db")
library(ggplot2)
library(RColorBrewer)
library(gridExtra)
library(enrichplot)
#library(plyr)
library(ggrepel)
library(stringr)
load(file = "../02_data/TCGA_TNBC_pair_annotation_result.Rdata")
mRNA = as.data.frame(diffSig_mRNA$gene_id %>%
                         str_split("\\|",simplify = T))
names(mRNA) = c("symbol", "probe_id","gene_type")
mRNA = na.omit(mRNA)
gene.df <- bitr(mRNA$symbol, fromType="SYMBOL",
                toType="ENTREZID", 
                OrgDb = "org.Hs.eg.db")
ego_BP <- enrichGO(gene          = gene.df$ENTREZID,
                   OrgDb         = org.Hs.eg.db,
                   keyType       = "ENTREZID",
                   ont           = "BP",
                   pvalueCutoff  = 0.05,
                   qvalueCutoff  = 0.05,
                   readable      = TRUE)

dotplot(ego_BP, showCategory = 20,font.size = 8)

img

# ?dotplot
ego_CC <- enrichGO(gene          = gene.df$ENTREZID,
                   OrgDb         = org.Hs.eg.db,
                   keyType       = "ENTREZID",
                   ont           = "CC",
                   pvalueCutoff  = 0.05,
                   qvalueCutoff  = 0.05,
                   readable      = TRUE)

dotplot(ego_CC, showCategory = 20,font.size = 8)

img

ego_MF <- enrichGO(gene          = gene.df$ENTREZID,
                   OrgDb         = org.Hs.eg.db,
                   keyType       = "ENTREZID",
                   ont           = "MF",
                   pvalueCutoff  = 0.05,
                   qvalueCutoff  = 0.05,
                   readable      = TRUE)

dotplot(ego_MF, showCategory = 20,font.size = 8)

img

ego_KEGG <- enrichKEGG(gene = gene.df$ENTREZID, organism = "hsa", 
                       keyType = "kegg",
                       pvalueCutoff = 0.05,
                       pAdjustMethod = "BH",
                       minGSSize = 10, maxGSSize = 500, 
                       qvalueCutoff = 0.05,
                       use_internal_data = FALSE)

dotplot(ego_KEGG, showCategory = 20,font.size = 8)

img

save(ego_BP, ego_CC,ego_MF,ego_KEGG, file = "../02_data/TCGA_TNBC_pair_enrichment.Rdata")

富集分析本身的难度不大,但是如何利用富集到的通路去讲你的生物学故事,则十分重要

因为内容太多,微信推文排版限制,第七步WGCNA需要大家移步到本期第3条查看,谢谢!

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2019-08-30

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券