疾病前期阶段到疾病发生的临界点预测以实现疾病预警（二）＜代码演示＞

1、数据分析环境

matlab和R语言

2、演示数据

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30550

该数据集包含流感病毒感染过程的时间序列数据，在该数据集中，17 名健康人类志愿者接种了 H3N2 流感病毒，并在 16 个时间点检测了每个个体的宿主外周血样本中的基因表达谱，17 名志愿者（受试者）中有 9 人出现了严重的流感症状，其余 8 人没有。

3、具体步骤

3.1、数据下载

a、先下载代码和参考样本

直接下载：https://github.com/Junxian-Li-0/sPGGM_project/archive/refs/heads/main.zip这个压缩包，里面包含了sPGGM计算代码以及癌症进展的数据。

或者使用git：git clone https://github.com/Junxian-Li-0/sPGGM_project.git

b、PPI网络下载：

到浏览器打开https://cn.string-db.org/cgi/download?sessionId=bxIUoWbWcLle，

下载PPI网络。

c、下载流感病毒感染过程的时间序列数据和预处理

数据下载这边用到r语言的GEOquery包，因为下完以后的数据不仅包含基因表达矩阵，还有整理好的临床信息，分组信息等。

library(GEOquery)
library(GEOquery)
library(clusterProfiler)
library(org.Hs.eg.db)
library(limma)
library(biomaRt)

#########下载流感病毒感染过程的时间序列数据
H3N4<-getGEO("GSE30550")#自动下载数据，如果下载失败可以到数据库下载到本地，在读取进来
expr<-exprs(H3N4[[1]])|>as.data.frame()##提取表达矩阵
fdata<-fData(H3N4[[1]])|>as.data.frame() ##提取平台文件
head(fdata)#里面的Gene_ID是Entrez ID，我们给它转成Symbol
rownames(expr)<-fdata$Gene_ID
genes<-bitr(rownames(expr),fromType = "ENTREZID",toType = "SYMBOL",OrgDb = org.Hs.eg.db)##ID转换
expr<-expr[genes$ENTREZID,]##过滤掉没有匹配上的
rownames(expr)<-genes$SYMBOL##将行名替换成Symbol
expr[1:4,1:4]

#########数据标准化
data_all<-normalizeBetweenArrays(expr)|>as.data.frame()##标准化
##将参考样本和目标样本分开
pdata<-pData(H3N4[[1]])|>as.data.frame()
ref_sample<-pdata$geo_accession[pdata$`time_hpi:ch1`=="Baseline"]
expr<-data_all[,!colnames(data_all)%in%ref_sample]
ref<-data_all[,ref_sample]
write.csv(ref,"ref.csv",quote = F)
write.csv(expr,"H3N4.csv",quote = F)

3.2、sPGGM计算

sPGGM用到matlab

matlab里将工作目录切换到代码所在的文件夹

case_data=readtable('H3N4.csv','ReadRowNames', true);
case_data=table2array(case_data);
case_data=case_data(:,1:32);
ref_data=readtable('ref.csv','ReadRowNames',true);
ref_data=table2array(ref_data);
network_path='network.txt';
local_sPGGM = get_LocalsPGGM(case_data,ref_data,network_path);
count=20 %选择前k个模块
patient_label=(1:32); 
patient_num=ones(32,1);
sPGGM = calc_GlobalsPGGM(local_sPGGM,count,patient_label,patient_num);
%% local_sPGGM: the output matrix from get_LocalsPGGM function
%% count: signaling molecules num
%% patient_label: the output array from get_stage function
%% patient_num: the output array from get_stage function
T = array2table(sPGGM);
writetable(T, 'global_sPGGM.csv');

 3.3、sPGGM可视化

library(GEOquery)
library(ggplot2)
H3N4<-getGEO("GSE30550")
meta_data<-pData(H3N4[[1]])
sPGGM<-read.csv("global_sPGGM.csv")|>t()|>as.data.frame()
sPGGM$sample_id<-rownames(meta_data)[1:32]
sPGGM$time<-meta_data$`time_hpi:ch1`[1:32]
sPGGM$time<-factor(sPGGM$time,levels = sPGGM$time[1:16])
sPGGM$sample<-gsub(",.*","",meta_data$title)[1:32]
ggplot(sPGGM,aes(x=time,y=V1,color=sample))+
  geom_point()+
  geom_line(aes(color=sample,group=sample),size=1)+
  theme_bw()+labs(y="sPGGM")+
  theme(axis.text.x = element_text(angle = 45,hjust=1))

Subject 01是中途出现症状的，Subject 02是没有出现症状的。从结果中可以看出Subject 01的sPGGM得分会出现剧烈波动。

此外还有另外一个基于python写的DNB的代码：https://colab.research.google.com/github/hiroshi-yamashita/dnb-tools/blob/master/tutorial_dnb_timeseries.ipynb，里面具体原理可以看一下它的文献。