网页工具能解决一切问题？

生信技能树

发布于 2022-07-26 10:35:29

2810

发布于 2022-07-26 10:35:29

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

收到学员提问，不知道是质疑网页工具还是质疑我们授课的数据分析代码：

上下调出现冲突

其实我压根不关系学员代码写的是否正确，或者说网页工具是否严谨，既然是在质疑表达量高低，那么很明显自己绘制箱线图肉眼确定一下即可，很明显是一个在肿瘤里面下调的基因啊。

所以，学员使用我们授课的数据分析代码针对这个数据集，可以看到确实是在肿瘤里面下调的基因：

在肿瘤里面下调的基因

但是网页工具里面，可以看到变化倍数是一模一样的，仅仅是上下调基因反过来了：

仅仅是上下调基因反过来了

最后的学员他自己绘制箱线图肉眼确定一下，就需要写代码了，很明显确实是在肿瘤里面下调的基因。这个地方，网页工具处理的不严谨。当然了，这里并不是说要把网页工具贬低的一无是处，只不过是如果你自己会写代码，你的选择更多一点。

所以我安排这个月的实习生处理了这个GSE78092数据集

首先，自己去GSE78092数据集下载一些文件，然后差异分析的代码很简单：

#01读取数据
library(data.table) #传说fread读数据是read.table的一百倍
probeM <- fread(file = "GSE78092_series_matrix.txt",
                header = T,sep="\t",check.names = F)
ID <- fread(file = "GPL21485.txt",
            header = T,sep = "\t",check.names = F)

#查看GSE关于样本分组信息可知前三个样本组是疾病GC组，后三个是正常Normal组
Group = c(rep("GC",times=3),rep("Normal",times=3))
Group = factor(Group,levels = c("Normal","GC"))
head(Group)  ##可以看见组别是疾病到正常，与探针表达矩阵一致

#ID与探针矩阵数目2902观测一致,验证内容是否一致
table(ID$ID%in%probeM$ID_REF)  #2902TRUE，说明一一对应

#获取cicRNA表达矩阵
cicrcM <- merge(ID,probeM,by.x="ID", by.y="ID_REF")

# 查看ASCRP000979这个基因的表达情况
cicrcM[cicrcM$ID=="ASCRP000979",]

library(limma)
design=model.matrix(~Group)
fit=lmFit(data,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)

感兴趣的小伙伴，也可以自己去检验一下哦！

既然同为芯片数据，那circRNA的处理与分析流程应该与基因表达谱芯片测序完全一样，区别只不过是一个通过注释平台转换为基因名，一个转换为circRNA名。今天我们关注的不是流程，而是有个小伙伴反应他用我们的芯片差异分析流程与网页差异分析得到的结果相反，那让我们看看是谁出错了。

在开始具体分析之前，我们得明确如何确定差异分析上下调是否发生了颠倒，究竟是网页工具有问题还是我们的脚本有问题。该小伙伴提供了网页版差异分析的结果与自身使用我们流程的结果。我们从该小伙伴提供的结果中随机挑选出那个除上下调有差异，变化倍数几乎完全一样的基因ASCRP000979，通过查看其原始表达量与我们差异分析的脚本进行明确。

说句题外话，之所以选择ASCRP000979是因为其是差异分析结果的第一个基因，并发生显著下调。好，那我们接下来开始三步分析。