顺着上周与大家分享的nanostring芯片原始数据的提取, 我看到了曾老师于20年布置的学徒作业~
链接如下:《Nanostring的表达矩阵分析也是大同小异》
Nanostring的表达矩阵分析也是大同小异 - 知乎 (zhihu.com)
就是要复现上图~ 草草一看应该是提取原始数据,取差异基因然后绘图吧。
文章的最后看到了老师的期望,就决定从原始数据开始分析。
那我们就开始吧
在补充材料里面看到的了GSE编号 GSE134129,虽然老师在知乎上也说了..
下载原始数据
这次的初始文件是txt格式的,之前是RCC格式的(这可能就是问题所在????)
rm(list=ls())
options(warn = -1)
library(nanostringr)
#文件解压
library(R.utils)
#filenames=dir(path = "./GSE134129_RAW/",pattern = ".txt.gz$")
#filenames=paste0("./GSE134129_RAW/",filenames)
#sapply(filenames, gunzip)
#文件读取合并成为表达矩阵
filenames2=dir(path = "./GSE134129_RAW/",pattern = ".txt$")
filepath=paste0("./GSE134129_RAW/",filenames2)
list=list()
for (i in 1:length(filepath)) {
list[[i]]=read.delim(file=filepath[i])
}
table(list[[1]][["ID_REF"]]==list[[19]][["ID_REF"]]) #判断一下,应该是一样的
new.list=lapply(list, function(li){
cbind(li[[5]])
})
data=data.frame(x=vector(length = 750))
for (i in 1:length(new.list)) {
data[,i]=as.vector(new.list[[i]])
}#整合数据
rownames(data)=list[[1]][["ID_REF"]] #修改行名
library(stringr)
colname=str_extract(filenames2,"GSM\\d*")
colnames(data)=colname
table(data[,1]==as.vector(list[[1]][["Control_1"]]))
table(data[,21]==as.vector(list[[21]][["STING.KO_4"]]))
exp=data
#获取分组
group=ifelse(str_detect(filenames2,"Control"),"Control",
ifelse(str_detect(filenames2,"STING_KO"),"KO","STING"))
table(group)
图注对火山图的描述
实验一共分3个组,作者分别两两差异分析,将差异结果放在一张火山图上 下载临床数据确认一下具体的实验操作
#看过文献,发现需要详细确认一下分组情况
suppressMessages (library(GEOquery))
gse_number = "GSE134129"
gset=getGEO(gse_number,destdir = ".",AnnotGPL = F,getGPL = F)
# Found 1 file(s)
# GSE134129_series_matrix.txt.gz
# Using locally cached version: ./GSE134129_series_matrix.txt.gz
gset=gset[[1]]
phenoDat <- pData(gset)
table(phenoDat$characteristics_ch1.4)
# treatment: PBS (Sham) treated treatment: STING agonist treated
# 13 8
phenoDat.new=cbind(phenoDat$title,phenoDat$characteristics_ch1.4)
分组情况如下~
#根据分组信息,获取exp.st,exp.ko,group.st,group.ko,组成list
exp.st=exp[,1:17]
dim(exp.st)
# > dim(exp.st)
# [1] 749 17
group.st=group[1:17]
table(group.st)
# group.st
# Control STING
# 9 8
exp.ko=exp[,-(10:17)]
dim(exp.ko)
# [1] 749 13
group.ko=group[-(10:17)]
table(group.ko)
# group.ko
# Control KO
# 9 4
#组成list
exp.list=list(exp.st=exp.st,exp.ko=exp.ko)
exp.list[[1]][1:4,1:4]
# GSM3937676 GSM3937677 GSM3937678 GSM3937679
# A2m 3.138 3.053 1.585 2.722
# Abca1 8.109 9.009 7.624 7.890
# Abcb1a 5.655 6.138 5.472 5.914
# Abcg1 8.175 8.249 7.748 8.208
group.list=list(group.st=group.st,group.ko=group.ko)
group.list[[1]][8:12]
# [1] "Control" "Control" "STING" "STING" "STING"
最后得到以下两个列表用于差异基因分析exp列表
group列表
##批量差异分析
library(limma)
deg.list=list[[]]
for (i in 1:length(exp.list)) {
exp=exp.list[[i]]
Group=group.list[[i]]
design=model.matrix(~Group)
fit=lmFit(exp,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
deg.list[[i]]=deg
}
差异分析结果列表
现在差异分析结果有了,看一下作者关注的基因的表达量与我是否一致。
看了一下图上标记的基因
# 上调
# Angpt1 Pdgfrb Cdh5 Sell Ccr7 Cd86
#下调
# Cdh5 Mcam Angpt1 Sell Vegfa
稍微比对一下
##差异分析结果校对
up.article=c( "Angpt1","Pdgfrb","Cdh5","Sell","Ccr7","Cd86")
up.list=deg.list[[1]]
up.list[rownames(up.list) %in% up.article,]
# logFC AveExpr t P.Value adj.P.Val B
# Cd86 1.3106 7.482 4.647 0.0001868 0.003638 0.8724
# Sell 1.5947 8.168 4.471 0.0002767 0.003638 0.5006
# Angpt1 1.1483 5.628 3.382 0.0032204 0.010221 -1.8104
# Pdgfrb 1.0758 6.130 2.933 0.0086995 0.018832 -2.7327
# Cdh5 1.0715 6.177 2.793 0.0117928 0.023554 -3.0120
# Ccr7 0.6802 4.225 1.770 0.0931071 0.119875 -4.8377
down.article=c( "Cdh5","Mcam","Angpt1","Sell","Vegfa")
down.list=deg.list[[2]]
down.list[rownames(down.list) %in% down.article,]
# logFC AveExpr t P.Value adj.P.Val B
# Sell -1.3211 7.011 -4.552 0.0004431 0.01043 -0.09605
# Mcam -0.7474 6.517 -3.232 0.0059733 0.05227 -2.66977
# Angpt1 -1.0136 4.776 -3.076 0.0081562 0.06234 -2.97260
# Cdh5 -0.5305 5.510 -2.885 0.0119261 0.08264 -3.33949
# Vegfa -0.6953 6.752 -2.365 0.0328493 0.14828 -4.29987
与文章相比趋势是一致的,但数值上有差异。
不,应该说是有明显差异
从需要复现的图中来看,Cdh5的LogFC值应该在3-4左右,但我的差异结果却只有1多一点。同时这个Ccr7也有点奇怪。
于是再次仔细看了看文章的分析,然后发现图的纵坐标标注是倍数的变化...
原来作者没有差异分析,只是将各组表达量进行了简单的相除....
是我从来没有听说过的操作了...
回到文章细看~
作者在文章的补充材料中只提供了表达量和P值,没有logFC的相关信息,一共750个基因
再细看文章,
发现数据被标准化了,也没有详说用的是哪种标准化方式。
关于目标图文章中为数不多的描述
那就从作者提供的矩阵开始复现吧。