把差异分析换一个单位
转录组差异分析大家应该是都不陌生了,无论是表达量芯片还是转录组测序,拿到了矩阵后下游分析无非就是选择不同的统计学R包,以及让人眼前一亮的可视化!但是常规的依据变化倍数和统计学p值来划分上下调基因列表和对这样的基因集进行go和kegg数据库注释大概率都是被做过了的,这里分享最近看到的一个还算是另辟蹊径的做法!
文献是2021发表在Nature Communications杂志的文章:《Individualized interactomes for network-based precision medicine in hypertrophic cardiomyopathy with implications for other clinical pathophenotypes》,链接是 https://pubmed.ncbi.nlm.nih.gov/33558530/
这个研究的数据在:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE160997
差异分析的实验设计策略是:18 HCM patients and compare them with those in 5 healthy controls.
结果是:We identified 2238 significantly differentially expressed genes between the control and HCM cohorts (1649 are protein-coding genes);
到这里仍然是平平无奇的转录组差异分析而已,但是认真看原文,就可以发现蛛丝马迹,它把差异分析换了一个单位,不再是纠结于具体的某个基因的高低变化,而是看任意的两两之间的基因组合的相关性是否在两个分组里面被破坏了!
先内部计算相关性再差异分析
这个RNA-Seq数据定量得到 44,285 genomic features,在5个正常样品表达量矩阵里面,首先删除0值后得到26,946 genes,再删除低表达基因后,剩下16,457 genes ,这些基因两两组合成为 135,408,196 correlations。
非常这样, 就是首先针对control分组的表达量矩阵,去除低表达量基因,然后两两之间组合成为相关性基因对!
![](http://n4.ikafan.com/assetsj/blank.gif)
然后考虑18个病人组成的另外一个分组,每次挑选一个病人到5个正常样品表达量矩阵里面去重新计算各个基因之间的相关性,就可以得到 patient-specific protein–protein interaction networks ,看看具体哪些基因因为增加了一个病人而不是正常人导致前面的基因配对相关性被很多程度的扰动了!
有意思的是作者给出来了代码:https://github.com/bwh784/HCM/blob/main/myPCC.R
library(WGCNA)
datExpr<-read.table("HCM_Filtered3CountData.txt")
datExpr_hcm<-datExpr[,1:18]
datExpr_ctrl<-datExpr[,19:23]
hcm<-cbind(datExpr_hcm[,1],datExpr_ctrl)
B<-transposeBigData(datExpr_ctrl)
ctrl_A<-adjacency(B,
selectCols = NULL,
type = "unsigned",
power = 1,
corFnc = "cor", corOptions = list(use = "p", method='pearson'),
weights = NULL)
hcm_B<-transposeBigData(hcm)
hcm_A<-adjacency(hcm_B,
selectCols = NULL,
type = "unsigned",
power = 1,
corFnc = "cor", corOptions = list(use = "p", method='pearson'),
weights = NULL)
for (i in 1:16456)
{
for (j in (i+1):16457)
{
delta_p=hcm_A[i,j]-ctrl_A[i,j]
dominator=(1-ctrl_A[i,j]*ctrl_A[i,j])/4
x=c()
p=pnorm(abs(delta_p),mean=0,sd=dominator,lower.tail = FALSE, log.p = FALSE)
if(!is.na(p)&p<3.69e-10&ctrl_A[i,j]<1)
{
x=cbind(x,i,j,p,delta_p,ctrl_A[i,j])
write.table(x, file="HCM1_significant_deltaP.txt",sep = "\t",append = TRUE, row.names = FALSE,col.names = FALSE)
}
}
}
但是我很讨厌这样的双层for循环的代码,仅仅是个人情绪问题,我比较倾向于使用apply系列函数来替代这样的for循环。
学徒作业
下载前面的文章介绍的数据集:GSE160997,里面的5个正常对照和18 HCM patients的表达量文件:
GSM4887610 anterior septal tissues - ctrl1
GSM4887611 anterior septal tissues - ctrl2
GSM4887612 anterior septal tissues - ctrl3
GSM4887613 anterior septal tissues - ctrl4
GSM4887614 anterior septal tissues - ctrl5
GSM4887615 anterior septal tissues - HCM1
GSM4887616 anterior septal tissues - HCM2
GSM4887617 anterior septal tissues - HCM3
GSM4887618 anterior septal tissues - HCM4
GSM4887619 anterior septal tissues - HCM5
GSM4887620 anterior septal tissues - HCM6
GSM4887621 anterior septal tissues - HCM7
GSM4887622 anterior septal tissues - HCM8
GSM4887623 anterior septal tissues - HCM9
GSM4887624 anterior septal tissues - HCM10
GSM4887625 anterior septal tissues - HCM11
GSM4887626 anterior septal tissues - HCM12
GSM4887627 anterior septal tissues - HCM13
GSM4887628 anterior septal tissues - HCM14
GSM4887629 anterior septal tissues - HCM15
GSM4887630 anterior septal tissues - HCM16
GSM4887631 anterior septal tissues - HCM17
GSM4887632 anterior septal tissues - HCM18
自己看文献,写代码,得到每个病人的特异性的被扰动的基因相关性配对情况。
如果你对差异分析都还不了解,可以先看看我们的b站免费课程:
全网最系统的表达芯片数据处理教程
表达芯片数据处理教程,早在2016年我就系统性整理了发布在生信菜鸟团博客:http://www.bio-info-trainee.com/2087.html
配套教学视频在B站:https://www.bilibili.com/video/av26731585/
代码都在:https://github.com/jmzeng1314/GEO
早期目录如下:
第一讲:GEO,表达芯片与R 第二讲:从GEO下载数据得到表达量矩阵 第三讲:对表达量矩阵用GSEA软件做分析 第四讲:根据分组信息做差异分析 第五讲:对差异基因结果做GO/KEGG超几何分布检验富集分析 第六讲:指定基因分组boxplot指定基因list画热图 第七讲:根据差异基因list获取string数据库的PPI网络数据 第八讲:PPI网络数据用R或者cytoscape画网络图 第九讲:网络图的子网络获取 第十讲:hug genes如何找
公众号推文在: