社区首页 >专栏 >你以为它是表达量芯片的原始信号值矩阵吗？

你以为它是表达量芯片的原始信号值矩阵吗？

生信菜鸟团

发布于 2023-09-09 09:07:10

26000

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

1-背景

顺着上周与大家分享的nanostring芯片原始数据的提取，我看到了曾老师于20年布置的学徒作业~

链接如下：《Nanostring的表达矩阵分析也是大同小异》

Nanostring的表达矩阵分析也是大同小异 - 知乎 (zhihu.com)

就是要复现上图~ 草草一看应该是提取原始数据，取差异基因然后绘图吧。

文章的最后看到了老师的期望，就决定从原始数据开始分析。

那我们就开始吧

2- 找原始数据

在补充材料里面看到的了GSE编号 GSE134129，虽然老师在知乎上也说了..

下载原始数据

3- 获取总表达矩阵

这次的初始文件是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

4- 获取总分组

#获取分组
group=ifelse(str_detect(filenames2,"Control"),"Control",
              ifelse(str_detect(filenames2,"STING_KO"),"KO","STING"))
table(group)

5- 提取exp&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列表

6- limma 批量差异分析

##批量差异分析
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
}

差异分析结果列表

7- 差异分析结果校对

现在差异分析结果有了，看一下作者关注的基因的表达量与我是否一致。

看了一下图上标记的基因

# 上调
# 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