【原】你研究的基因凭什么重要（这才是数据挖掘的用武之地）

不少小伙伴找到我们教学团队想学生物信息学数据挖掘的时候总是喜欢问一下，能否发文章，更有甚者希望我们的教学承诺包发文章，我就呵呵了，首先这个承诺包发文章是否违规违法不说起码也是一个灰色地带，我们是不会做的。我们直播互动教学是希望授人以渔，在数据挖掘课程方面，总体上来说就是希望大家有可以证明你研究的目标基因的重要性。通俗一点就是你的文章的figure1.

比如2022年6月的文章：《Restoring the epigenetically silenced lncRNA COL18A1-AS1 represses ccRCC progression by lipid browning via miR-1286/KLF12 axis》，其生物学故事的落脚点是一个lncRNA，名字是COL18A1-AS1 ，文章用figure1说明了它在肿瘤里面表达量下降了，而且具体到每个病人的配对正常组织和肿瘤组织都可以看到这个规律。好好的一个 lncRNA COL18A1-AS1为什么在肿瘤里面表达量就降低了呢？然后看了看肿瘤里面这个基因也确实是一个保护因子，它表达量越高的话病人生存越好，而且可以看到无论是TNM分期，还是stage和grade这样的分级，都是肿瘤病人情况越糟糕，它表达量越低。

figure1

事实上，经过了我们马拉松基础课程一个月勤学苦练的小伙伴基本上都可以完成，我把这个作业布置给了十几个学员，完成率90%，任意选择了一个“新乡医学院”的本科刚刚毕业的小伙伴作业给大家点评一下。

小作业

Fig. 1: COL18A1-AS1 was downregulated and could be a biomarker for ccRCC. 里面的 A到 J的9个图 文献链接：https://www./articles/s41419-022-04996-2

TAGC-KIRC

0.运行环境与输入数据的准备

### 1.清空环境，加载包
rm(list=ls())
library(tidyverse)
library(data.table) 

### 2.自定义去重基因函数，用于后续基因去重
deDuplicateGene <- function(data){
  index <- order(rowMeans(data[,-1]),decreasing = T)
  data_ordered <- data[index,]
  keep <- !duplicated(data_ordered[1])
  data <- data_ordered[keep,] 
  row.names(data) <- data[,1]
  data <- data[,-1]
}

### 3.下载数据
proj = "TCGA-KIRC"
if(!file.exists(paste0(proj, ".htseq_counts.tsv.gz"))){
 download.file(url = paste0("https://gdc./download/",proj, ".htseq_counts.tsv.gz"),destfile = paste0(proj,".htseq_counts.tsv.gz"))
 download.file(url = paste0("https://gdc./download/",proj, ".GDC_phenotype.tsv.gz"),destfile = paste0(proj,".GDC_phenotype.tsv.gz"))
 download.file(url = paste0("https://gdc./download/",proj, ".survival.tsv"),destfile = paste0(proj,".survival.tsv"))
  }
  
### 4.整理数据（获取后续分析的基因cpm表达矩阵）
data_meta <- fread(paste0(proj,".GDC_phenotype.tsv.gz"), check.names = F) %>% 
  select(., contains(c('submitter','pathologic_T', 'pathologic_N',
  'pathologic_M', 'histologic_grade','tumor_stage'))) %>% 
  set_names(., c('sample','simple.ID','T','N','M','grade','stage')) %>% 
  mutate( T = case_when( T == 'T1'  ~ 'T1',
 T == 'T1a' ~ 'T1',
 T == 'T1b' ~ 'T1',
 T == 'T2'  ~ 'T2',
 T == 'T2a' ~ 'T2',
 T == 'T2b' ~ 'T2',
 T == 'T3'  ~ 'T3',
 T == 'T3a' ~ 'T3',
 T == 'T3b' ~ 'T3',
 T == 'T3c' ~  'T3')) %>% 
  left_join(.,(fread(paste0(proj,".survival.tsv"),check.names = F) %>% select(-3)))

data_exp <- fread(paste0(proj,".htseq_counts.tsv.gz"),
   check.names = F, data.table = F) %>% 
  mutate( Ensembl_ID = str_split(Ensembl_ID, '\\.', simplify = T) %>% .[,1]) 

gene_anno <- AnnoProbe::annoGene( data_exp$Ensembl_ID, ID_type = 'ENSEMBL') %>%
  select(1,3) %>% set_names(., c('SYMBOL','Ensembl_ID')) 

data_exp <- merge(data_exp, gene_anno, by = 'Ensembl_ID') %>% 
  select(-Ensembl_ID) %>% select(SYMBOL,everything()) %>% 
  deDuplicateGene()
dat = as.matrix(2^data_exp - 1) #逆转log
exp = apply(dat, 2, as.integer) # # 用apply转换为整数矩阵，行名消失
rownames(exp) = rownames(dat)
exp= log(edgeR::cpm(exp)+1) #获取cpm表达矩阵
data_exp=exp
#获取组别信息，将样本分为正常组与肿瘤组
Group <- factor(ifelse(as.numeric(substr(colnames(data_exp),14,15)) < 10,'tumor','normal'))  
Group = factor(Group,levels = c("normal","tumor"))
table(Group) 
## Group
## normal  tumor 
##  72 535
save(data_exp,Group,data_meta,proj,file = 'input.Rdata') 

p1绘制72例正常肾脏组织与522例ccRCC的COL18A1-AS1表达图

# 1.清空环境，加载包，加载数据
rm(list=ls())
library(ggplot2)
library(ggsignif)
library(ggpubr)
load(file = 'input.Rdata') 
  
# 2.整理数据
exp <- data_exp
dat=as.data.frame(exp["COL18A1-AS1",])
dat$group=Group
colnames(dat)=c("COL18A1-AS1","group")

# 3.画图
p1 <- ggboxplot(dat,x="group",y="COL18A1-AS1",palette = c("#00AFBB", "#E7B800"),color = "group")+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(title ="TCGA-KIRC",x="",y="expression level")+
  theme_bw()
p1  

p2绘制71对配对正常与肿瘤样本的COL18A1-AS1基因表达

# 1获取正常样本与肿瘤样本对应的样本ID及矩阵
data_plot=dat
data_pair1 <- data_plot[dat$group == 'normal', ]
data_pair1_id <- data_plot[dat$group == 'normal',] %>% rownames(.) %>% str_sub(., 1, 12) 
data_pair2 <- data_plot[dat$group == 'tumor', ]
data_pair2_id <- data_plot[dat$group == 'tumor', ] %>% rownames(.) %>% str_sub(., 1, 12) 

# 2.取配对样本则意味着配对的正常组与肿瘤组前12位相同，即可以取交集
inters <- intersect(data_pair1_id,data_pair2_id)

# 3.再根据正常组与肿瘤组特性取出各自的样本，组成新的绘图矩阵
N_id <- paste0(inters,'-01A')
p_id =  paste0(inters,'-11A')
datN <- data_plot[rownames(data_plot)%in%N_id,] #一个71
datP <- data_plot[rownames(data_plot)%in%p_id,] #一个72

# 4.去掉NA样本，获取真正配对样本，将配对样本组成新矩阵
nodup <- datN[!str_sub(rownames(datP), 1, 12)%in%str_sub(rownames(datN), 1, 12),]#找不配对的
datP <- datP[str_sub(rownames(datP), 1, 12)%in%str_sub(rownames(datN), 1, 12),]
datPair <- rbind(datN,datP)

# 5.设置比较组别，绘制配对样本COL18A1-AS1基因表达图
compaired <- list(c("normal", "tumor")) 
ggpaired(datPair, x = "group", y = "COL18A1-AS1",
   color = "group", 
   line.color = "black", 
   line.size = 0.4,
   palette = "jco")+ 
  labs(x = NULL, y = "Relative COL18A1-AS1\nexpression level", title="TCGA-KIRC",tag = "B")+
  theme(plot.title = element_text(hjust = 0.5, size = 20),
  plot.tag = element_text(size = 25),
  axis.title.y = element_text(size = 20, 
   color = "black"),
  legend.text = element_text(size = 7),
  legend.position = "none",
  aspect.ratio=.6)+
  stat_compare_means(comparisons=compaired, method = "wilcox.test",
   label ="p.signif")

# 6.保存后续ROC绘图所需的Rdata
save(data_plot,dat,file="input2.Rdata")

p3基于COL18A1-AS1基因表达量分组后生存分析图

# 1.清空环境，加载包，加载数据
rm(list=ls())
library(survival)
library(survminer)
load("input.Rdata")
load("input2.Rdata")

# 2.数据整理
data_sur1 <- data_plot[dat$group == 'tumor', ] #只要肿瘤组
data_sur1 <- as.data.frame(data_sur1)
## 根据COL18A1-AS1的表达与中位数的表达进行生存分析分组构建
data_sur1 <- mutate(data_sur1,
  gene = ifelse(data_sur1$`COL18A1-AS1`<median(data_sur1$`COL18A1-AS1`),"low","high"))

# 3.过滤不符合要求的样本
## 01以样本为主，去除肿瘤重复，取样本ID不一致的肿瘤样本
k = !duplicated(str_sub(rownames(data_sur1),1,12)) 
table(k)
## k
## FALSE  TRUE 
##  5   530
data_sur1<-data_sur1[k,]
data_sur1$sample<-rownames(data_sur1)
## 02取有生存信息的肿瘤样本
k1<-data_meta$sample %in%data_sur1$sample 
data_meta2<-data_meta[k1,]
mg <- merge(data_sur1,data_meta2,by="sample") #拼接成包含临床信息与表达信息的矩阵
## 03以生存时间与生存结局为条件过滤
max(mg$OS.time) #断定其以日为单位，日的话过滤条件取30
## [1] NA
m1 = mg$OS.time >= 30;table(m1)
## m1
## FALSE  TRUE 
## 13   513
m2 = !(is.na(mg$OS.time)|is.na(mg$OS));table(m2)
## m2
## FALSE  TRUE 
##  4   526
mg = mg[m1&m2,]

# 4绘制生存曲线
fit <- survfit(Surv(OS.time, OS) ~ gene, data =mg) #生存建模
ggsurvplot(fit,
  pval = TRUE, conf.int = TRUE,
  risk.table = TRUE, 
  risk.table.col = "strata", 
  linetype = "strata",
  surv.median.line = "hv", 
  ggtheme = theme_bw(), 
  palette = c("#E7B800", "#2E9FDF"),title = "TCGA-KIRC-OS")

p4基于DFS的生存曲线绘制

# 获取含DFS信息的矩阵（从Cell整理的泛癌矩阵寻找）
pancancer=fread("Survival_SupplementalTable_S1_20171025_xena_sp")# 这个不建议运行
kirc=pancancer[pancancer$sample%in%substr(rownames(data_sur1),1,15),]
kirc$sample=paste0(kirc$sample,"A")
m3 = kirc$DSS.time >= 30;table(m3)
## m3
## FALSE  TRUE 
## 15   515
m4 = !(is.na(kirc$DSS.time)|is.na(kirc$DSS));table(m4)
## m4
## FALSE  TRUE 
## 11   519
metag = mg[m3&m4,]
inters <- intersect(metag$sample,kirc$sample)
metag <- metag[metag$sample%in%inters,]
kircg <- kirc[kirc$sample%in%inters,]
kircg <- as.data.frame(kircg)
DFS <- kircg[,c("sample","DSS","DSS.time")]
metaDFS <- merge(metag,DFS)  #构成包含DSS的矩阵
metaDFS=na.omit(metaDFS)
# 准备绘制DSS的生存曲线并保存
sfitDFS=survfit(Surv(DSS.time, DSS)~gene, data=metaDFS)
ggsurvplot(sfitDFS,pval =TRUE, data = metaDFS, risk.table = TRUE)

ggsave(paste0(proj,"_DFS_vo.png"),width = 15,height = 10) 

p5基于GRADE的箱线图

# 重点就是提取与整理GRADE信息，构建绘图数据，进行绘图

# 1.输入数据与GRADE初整理
data_plot <- data_plot %>% as.data.frame(.) 
data_plot$sample<-rownames(data_plot)
data_meta$grade[grepl('11A',data_meta$sample)]='normal'
data_meta3 <- merge(data_plot,data_meta,by="sample")
head(trimws(data_meta3$`COL18A1-AS1`)) # 删除空格
## [1] "0.472770378491576"  "2.15853432121316"   "0.129533994460106" 
## [4] "0.282717185208964"  "0.130205363207768"  "0.0425249765737968"
data_meta3$`COL18A1-AS1` =  as.numeric(trimws(data_meta3$`COL18A1-AS1`))

#  2.整理GRADE分期信息
kp = data_meta3$grade %in% c('normal','G1','G2','G3','G4')
data_meta3 = data_meta3[kp,]
data_meta3$grade = factor(data_meta3$grade,
  levels = c('normal','G1','G2','G3','G4'))

#  3.绘制基于GRADE的箱线图
box1 <- ggboxplot(data_meta3, x = 'grade', y = 'COL18A1-AS1',
   color = "grade",title ="TCGA-KIRC-Grade")+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(x="",y="expression level")+
  theme_bw()
box1  

p6基于T分期的箱线图

#  1.输入数据与T分期整理
data_meta3$T[grepl('11A',data_meta3$sample)]='normal'
table(data_meta3$T) 
## 
## normal  T1  T2  T3 
##  72 271  67 178
data_meta3=na.omit(data_meta3)

#  2.绘制箱线图
box2 <- ggboxplot(data_meta3, x = 'T', y = 'COL18A1-AS1',
   color = "T",title ="TCGA-KIRC-T",)+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(x="",y="expression level")+
  theme_bw()
box2

p7基于N分期的箱线图

## 1.按照N分期整理
data_meta3$N[grepl('11A',data_meta3$sample)]='normal'
table(data_meta$N) 
## 
##  N0  N1  NX 
## 446  30 509

## 2.绘制基于N分期的箱线图
box3 <- ggboxplot(data_meta3, x = 'N', y = 'COL18A1-AS1',
   color = "N",title ="TCGA-KIRC-N",)+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(x="",y="expression level")+
  theme_bw() 
box3

p8基于M分期的箱线图

# 1.将正常组对应的M分期里的行名修改
data_meta3$M[grepl('11A',data_meta3$sample)]='normal'
table(data_meta$M)  
## 
##   M0  M1  MX 
##   2 794 155  34

# 2.绘制基于M分期的箱线图
box4 <- ggboxplot(data_meta3, x = 'M', y = 'COL18A1-AS1',
   color = "M",title ="TCGA-KIRC-M",)+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(x="",y="expression level")+
  theme_bw() 
box4

p9基于stage的箱线图

# 按照stage分期整理
data_meta3$stage[grepl('11A',data_meta3$sample)]='normal'
table(data_meta3$stage)
## 
## normal not reported   stage i  stage ii stage iii  stage iv 
##  71   3 261  55 122  71
data_meta3 = filter(data_meta3, data_meta3$stage !='not reported') 

# 绘制stage分期箱线图
box5<- ggboxplot(data_meta3, x = 'stage', y = 'COL18A1-AS1',
  color = "stage",title ="TCGA-KIRC-stage",)+
  stat_compare_means(method = "t.test",label = "p.signif")+
  labs(x="",y="expression level")+
  theme_bw() 
box5

p10基于指定基因的分类模型评估(ROC曲线)

# 准备ROC曲线绘制所需的数据（xdata，ydata）
rm(list = ls())
load("input2.Rdata")
library(pROC) 
library(ggplot2) 
xdata <-  data_plot$group
ydata <- data_plot$`COL18A1-AS1`
m <- pROC::roc( xdata ,  ydata )
auc(m)
## Area under the curve: 0.9153

# 绘制ROC图
g <- ggroc(list(m ),legacy.axes = T,size = 1) 
gg=g + theme_minimal() +
  scale_color_manual(values = c("#2fa1dd"))+
  geom_segment(aes(x = 0, xend = 1, y = 0, yend = 1), 
   colour = "grey", linetype = "dashed")+
  annotate("text",x = .75, y = .25,
  label = paste("AUC of m = ",format(round(as.numeric(auc(m)),4),nsmall = 2)),color = "#2fa1dd")
gg 

p11.reportROC绘制ROC曲线 

library(reportROC)
reportROC(gold=data_plot$group,
 predictor=1-data_plot$`COL18A1-AS1`,
 important="se",
 plot=TRUE)

##  Cutoff   AUC AUC.SE AUC.low AUC.up     P   ACC ACC.low ACC.up   SEN SEN.low
##   0.386 0.915  0.020   0.876  0.955 0.000 0.883   0.883  0.883 0.882   0.855
##  SEN.up   SPE SPE.low SPE.up   PLR PLR.low PLR.up   NLR NLR.low NLR.up   PPV
##   0.910 0.889   0.816  0.961 7.940   4.128 15.272 0.132   0.104  0.169 0.983
##  PPV.low PPV.up   NPV NPV.low NPV.up   PPA PPA.low PPA.up   NPA NPA.low NPA.up
##    0.972  0.995 0.504   0.417  0.591 0.882   0.855  0.910 0.889   0.816  0.961
##    TPA TPA.low TPA.up KAPPA KAPPA.low KAPPA.up
##  0.883   0.857  0.909 0.580     0.494    0.665

如果你也感兴趣这样的数据挖掘技能，目前只能先预定明年的马拉松授课名额：提前锁定年后马拉松授课名额，赠送价值1600的单细胞数据分析一次

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

•  数据挖掘（GEO,TCGA,单细胞）2022年暑期班（收官之作），快速了解一些生物信息学应用图表
•  生信入门课-2022年暑期班（收官之作），你的生物信息学第一课

写在文末

我在《生信技能树》，《生信菜鸟团》，《单细胞天地》的大量推文教程里面共享的代码都是复制粘贴即可使用的， 有任何疑问欢迎留言讨论，也可以发邮件给我，详细描述你遇到的困难的前因后果给我，我的邮箱地址是 [e-mail]

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。