在学术不端的数据取舍上面反复横跳

生信技能树

发布于 2024-01-11 15:24:36

27100

代码可运行

文章被收录于专栏：生信技能树

运行总次数：0

代码可运行

前面我在《单细胞天地》分享了自己整理的常规单细胞转录组数据分析文章的3个规律，首先都是第一层次降维聚类分群。然后第一种文章是再次对每个亚群继续细分走三件套（亚群注释，拟时序，转录因子）这样就可以凑三五个图表。第二种文章是针对其中一个亚群探索三五个大图。最后一种文章就是结合各种公共数据，包括同样的实验设计的bulk转录组，预后信息，分子分型等分类信息。

然后马上这些策略就被应用到了单细胞转录组数据挖掘层面，因为反正也不需要自己产出数据了，过去三五年间单细胞的火热带动了海量的各种实验设计的公开的表达量矩阵。比如这个文献：《Lipid-related protein NECTIN2 is an important marker in the progression of carotid atherosclerosis: An intersection of clinical and basic studies》就是看了看两个分组的具体的基因的差异，在普通bulk转录组和单细胞转录组两个数据集里面，如下所示：

普通bulk转录组和单细胞转录组两个数据集

首先是单细胞转录组数据集，因为是两个分组所以作者进行了两次独立的降维聚类分群，然后看了看自己的目标基因nectin cell adhesion molecule 2 (NECTIN2) 确实是有特别的表现形式：

NECTIN2 was found to be highly expressed in ECs in the AC and PA regions and in macrophages located in the AC region

也就是说这个NECTIN2基因确实是在AC region是独特的高表达的，在单细胞水平。而且去普通bulk转录组验证，但是作者使用的是GEO2R这样的网页工具（小声比比，这里并没有说网页工具不好的意思哈），也顺利的说明了NECTIN2的高表达这个现象啦：

NECTIN2的高表达

我就随手把这个图表复现安排给了这一期的马拉松授课学员，大家群策群力，居然发现了里面的弯弯绕绕，确实是有猫腻啊！

简单的芯片表达量差异分析

我们这里不使用GEO2R这样的网页工具，直接开始敲代码啦！

rm(list = ls()) 
options(stringsAsFactors = F) 
GSE_number<-"GSE28829"

library(GEOquery)
library(stringr) 
getGEO(filename = paste0(GSE_number,"_series_matrix.txt.gz"),destdir = "./",getGPL = F)->geoObject
#提取临床信息
pData(geoObject)->pd
exprs(geoObject)->exp
geoObject@annotation->GPL
boxplot(exp)
save(exp,pd,GPL,file = 'pre.Rdata')

#检查下pd的每行和exp的每列样本是否对应
table(rownames(pd)==colnames(exp))

#提取简化group_list的分组信息
group_list<-pd$`phenotype:ch1`
table(group_list)
group_list[str_detect(group_list,"early")]<-"early"
group_list[str_detect(group_list,"advanced")]<-"advanced"

可以看到，这个GSE28829数据集一共29个样本，值得注意的是文献筛选了14个样本却没有说为什么。

我们先拿全部29个样本的两个分组走差异表达分析

load(file = 'pre.Rdata')
#若load以下代码则使用作者手动过滤后的14个样本做
#load(file ='filted.Rdata' )
#箱线图看着蛮齐的，做不做下面的影响不大
#limma::normalizeBetweenArrays(exp)->dat
exp->dat
#提取，简化分组信息
group_list<-pd$`phenotype:ch1`
table(group_list)
group_list[str_detect(group_list,"early")]<-"early"
group_list[str_detect(group_list,"advanced")]<-"advanced"
table(group_list)

探针id转化，每个芯片平台都有自己的方法，我们这里使用AnnoProbe包即可：

library(AnnoProbe)
ids=idmap(GPL,'soft')
head(ids)
ids=ids[ids$symbol != '',]

ids[str_detect(ids$symbol,"NECTIN2"),]
# [1] ID     symbol
# <0 行> (或0-长度的row.names)
ids[str_detect(ids$symbol,"Nectin2"),]
# [1] ID     symbol
# <0 行> (或0-长度的row.names)

但是文献中的基因名称NECTIN2，转化完id之后找不到。所以我们去NCBI找找他的别名，其中PVRL2能找到对应有探针：

所以接下来就使用这个PVRL2即可。

#有些探针对应不到基因，去掉symbol为空的
dat=dat[rownames(dat) %in% ids$ID,]
ids=ids[match(rownames(dat),ids$ID),]
head(ids) 
colnames(ids)=c('probe_id','symbol')  
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
dat[1:4,1:4] 
#这时候ids与dat行顺序完全一致
identical(ids$probe_id,rownames(dat))

因为有多个探针对应同一个基因的情况，所以要去重。方法一是所有symbol重复的行中，取中位数最大的作为基因表达量。方法二是取同symbol名行的平均数作为基因表达量。（比较麻烦吧，这里只对NECTIN2（既PVRL2）取了平均）

ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids

ids[str_detect(ids$symbol,"PVRL2"),]#查看PVRL2基因对应的探针及表达量
ids[str_detect(ids$symbol,"PVRL2"),]
#          probe_id symbol   median
# 12598   203149_at  PVRL2 7.386974
# 34676   225418_at  PVRL2 7.111999
# 41333   232078_at  PVRL2 4.986980
# 41334 232079_s_at  PVRL2 4.693793

rbind(exp["203149_at",],exp["225418_at",],exp["232078_at",],exp["232079_s_at",])->PVRL2exp
PVRL2exp
apply(PVRL2exp,2,mean)->PVRL2exp
PVRL2exp 

ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
# dat["PVRL2",]<-PVRL2exp #这一行是使用PVRL2的四个的探针的平均表达量替换中位数最高的探针的表达量
# 以防表达量在中位数最高的探针的组间没有差异，而在表达量较低的剩下三个探针的组间差异显著
# 防止冤假错案
save(dat,group_list,PVRL2exp,file = 'step1-output.Rdata')

差异表达分析：

rm(list = ls())
options(stringsAsFactors = F)
library(clusterProfiler)
library(stringr)
library(limma)

load(file = 'step1-output.Rdata')
table(group_list)
# 每次都要检查数据
dat[1:4,1:4]
##################跑下面这一步的话就是使用四行的平均表达量作为PVRL2的表达量
dat["PVRL2",]<-PVRL2exp

#NECTIN2
pro='advanced-vs-early' 
vs=str_replace(pro,"-vs-","-")
dir.create(pro)
setwd(pro) 
gse_number = pro
path="./"


library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)
exprSet=dat
rownames(design)=colnames(exprSet)
head(design)
contrast.matrix<-makeContrasts(vs,
                               levels = design)
contrast.matrix ##这个矩阵声明，我们要把实验组跟对照组进行差异分析比较
deg = function(exprSet,design,contrast.matrix){
  #step1
  fit <- lmFit(exprSet,design)
  #step2
  fit2 <- contrasts.fit(fit, contrast.matrix) 
  fit2 <- eBayes(fit2)  ## default no trend !!!
  #step3
  tempOutput = topTable(fit2, coef=1, n=Inf)
  nrDEG = na.omit(tempOutput) 
  head(nrDEG)
  return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)

## 火山图 
  nrDEG=deg
  head(nrDEG)
  attach(nrDEG)
  plot(logFC,-log10(P.Value))
  df=nrDEG
  cut_logFC <- with(df,mean(abs(df$logFC)) + 2*sd(abs(df$logFC)) )
  if(cut_logFC > 1){logFC_t = 1}
  if(cut_logFC < 1){logFC_t = cut_logFC}
  pvalue_t = 0.05
  df$v= -log10(P.Value) #df新增加一列'v',值为-log10(P.Value)
  df$g=ifelse(df$P.Value > pvalue_t,'stable', #if 判断：如果这一基因的P.Value>0.01，则为stable基因
              ifelse( df$logFC > logFC_t,'up', #接上句else 否则：接下来开始判断那些P.Value<0.01的基因，再if 判断：如果logFC >1.5,则为up（上调）基因
                      ifelse( df$logFC < -logFC_t,'down','stable') )#接上句else 否则：接下来开始判断那些logFC <1.5 的基因，再if 判断：如果logFC <1.5，则为down（下调）基因，否则为stable基因
  )
  table(df$g)
  df$name=rownames(df)
  head(df)
  
  this_tile <- paste0('Threshold of logFC is ',round(logFC_t,3),
                      '\nThe number of up gene is ',nrow(df[df$g == 'up',]) ,
                      '\nThe number of down gene is ',nrow(df[df$g == 'down',])
  )


target_gene="PVRL2"
for_label <-  df %>% 
  dplyr::filter(df$name %in% target_gene) 
head(for_label)
p1 <- ggplot(data = df, 
             aes(x = logFC, 
                 y = -log10(P.Value))) +
  geom_point(alpha=0.6, size=1.5, 
             aes(color=g)) +
  geom_point(size = 3, shape = 1, data = for_label) +
  ggrepel::geom_label_repel(aes(label = name),data = for_label,
                            max.overlaps = getOption("ggrepel.max.overlaps", default = 20),
                            color="black") +
  ylab("-log10(Pvalue)")+
  scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+
  geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) +
  #geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
  #geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
  #coord_flip()+
  xlim(-5, 5)+
  theme_bw()+
  ggtitle(this_tile )+
  theme(panel.grid = element_blank(),
        plot.title = element_text(size=8,hjust = 0.5),
        legend.title = element_blank(),
        legend.text = element_text(size=8))
p1
ggsave(filename = "volcanoplot.png",p1)

bardf<-data.frame(colnames(dat),dat["PVRL2",],as.factor(group_list))
colnames(bardf)<-c("GSE_number","expression","group")
p2<-ggplot(bardf, aes(x = GSE_number,y = expression))+ 
  theme_classic(base_line_size = 1)+
  geom_bar(aes(fill=group),stat = 'identity')+
  theme(axis.text.x = element_text(size = 13,  # 修改X轴上字体大小，
                             color = "black", # 颜色
                             face = "bold", #  face取值：plain普通，bold加粗，italic斜体，bold.italic斜体加粗
                             vjust = 0.5, # 位置
                             hjust = 0.5,
                             angle = 90))+#角度
  coord_cartesian(ylim = c(5.2, 6.8))
ggsave(filename = "barplot.png",p2)
setwd('../')
getwd()

使用表达量中位数最高的行作为该基因的表达量时，NECTIN2（PVRL2）就离上调很不沾边，在火山图也可以看到它远达不到统计学显著的上调基因的分类标准。

挑选样品

前面我们提到了，这个GSE28829数据集一共29个样本，但是文献筛选了14个样本。。。

（这里我们看不出来筛选的依据！原始数据中的临床信息以及文献中都没有说明）

不过为了复现图表，我们还是不得不屈服于文献的挑选策略：

filt<-c('GSM714070','GSM714072','GSM714073','GSM714077','GSM714078','GSM714079','GSM714080','GSM714086','GSM714088','GSM714090'
,'GSM714094','GSM714095','GSM714097','GSM714098')
pd[rownames(pd)%in%filt,]->pd
exp[,colnames(exp)%in%filt]->exp
table(rownames(pd)==colnames(exp))
save(exp,pd,GPL,file = 'filted.Rdata')

有了这个小的表达量矩阵，然后可以做同样的分析。从上文load(file='filted.Rdata')开始再做一遍看看，复现得到的图片如下。

这里用的是平均数去重的方法，可以发现NECTIN2(PVRL2)已经很接近统计学显著（两个阈值：变化倍数和矫正的p值）上调了。这个时候只需要把logFC的阈值减小就可以当作是上调了。看上图我甚至觉得文献中作者压根没设logFC（若设置阈值纵轴肯定会有两条明显的线区分上下调和稳定表达，他的图里面只有横轴的pvalue为0.05时的-log10Pvalue的阈值）

然后文献中和我们这个时候复现的D图看着是不是实验和对照差异挺明显的，其实是因为坐标轴的问题。如果我们调整y轴让它从零开始，如下所示：

调整y轴让它从零开始

同一个数据出的两张图，却给读者展示两种不同的结果。现在你还能看到明显的差异吗？

那么问题来了，用什么理由可以合理的过滤数据样品呢？

其中这个数据集小任务是给上个月的马拉松授课学员布置的，当时就给大家讲清楚了。而今天之所以想到这个话题了，是因为这一期又是有学员遇到了同样的问题，在分析自己感兴趣的数据集的时候发现计算得到的全部的基因的变化倍数情况小的可怜，如果是调整为0.5的阈值还勉强有一些差异分析结果，如果是调整为1那就约等于全军覆没啦！所以我就让学员给了三张图我们看看，如下所示的pca分析图：

pca分析图

我在生信技能树的教程：《你确定你的差异基因找对了吗？》提到过，必须要对你的转录水平的全局表达矩阵做好质量控制，最好是看到标准3张图：

左边的热图，说明我们实验的两个分组，normal和npc的很多基因表达量是有明显差异的
中间的PCA图，说明我们的normal和npc两个分组非常明显的差异
右边的层次聚类也是如此，说明我们的normal和npc两个分组非常明显的差异

如果分组在3张图里面体现不出来，实际上后续差异分析是有风险的。这个时候需要根据你自己不合格的3张图，仔细探索哪些样本是离群点，自行查询中间过程可能的问题所在，或者检查是否有其它混杂因素，都是会影响我们的差异分析结果的生物学解释。

所以，如果我们看到了这样的三张图有问题的实验设计项目的公共数据集，理论上我们需要根据很多实际情况进行样品的筛选。但是又容易在在学术不端的数据取舍上面反复横跳。。。

学徒作业

如上所示的数据集：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE25906

是 genome-wide gene expression profiles of a large (n=60) set of human placentas ，有很明显的两个分组：23 preeclamptic, and 37 control placentas

使用的是 Illumina human-6 v2.0 expression beadchip ，所以可以下载那个22.3 Mb的文件：GSE25906_normalized.txt.gz 走我们的简单的芯片表达量差异分析，但是如何取舍样品，就需要大家熟读文献，理解这个数据集啦：

Transcriptional profiling of human placentas from pregnancies complicated by preeclampsia reveals disregulation of sialic acid acetylesterase and immune signalling pathways. Placenta 2011 Feb;32(2):175-82. PMID: 21183218

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2024-01-07，如有侵权请联系 cloudcommunity@tencent.com 删除

芯片