Tag Archives: 统计

四 14

用R语言做逻辑回归分析

Posted on by

回归的本质是建立一个模型用来预测，而逻辑回归的独特性在于，预测的结果是只能有两种，true or false

在R里面做逻辑回归也很简单，只需要构造好数据集，然后用glm函数(广义线性模型（generalized linear model）)建模即可,预测用predict函数。

我这里简单讲一个例子，来自于加州大学洛杉矶分校的课程

Continue reading →

十二 21

用重抽样+主成分方法来做富集分析

Posted on by

之前我们用超几何分布检验的方法做了富集分析，使用的是GSE63067.diffexp.NASH-normal.txt的logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因来检验kegg的pathway的富集情况

结果是这样的

image001

我们接下来用另外一种方法来做富集分析，顺便检验一下，是不是超几何分布统计检验的富集分析方法就是最好的呢？

这种方法是-重抽样+主成分分析

大概的原理是，比如对上图中的，04380这条pathway来说，总共有128个基因，那么我从原来的表达矩阵里面随机抽取128个基因的表达矩阵做主成分分析，并且抽取一千次，每次主成分分析都可以得到第一主成分的贡献度值。那么，当我并不是随机抽取的时候，我就抽04380这条pathway的128个基因，也做主成分分析，并且计算得到第一主成分分析的重要性值。我们看看这个值，跟随机抽1000次得到的值差别大不大。

这时候就需要用到表达矩阵啦！

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

exprSet=read.table("GSE63067_series_matrix.txt.gz",comment.char = "!",stringsAsFactors=F,header=T)

rownames(exprSet)=exprSet[,1]

exprSet=exprSet[,-1]

我们根据ncbi里面对GSE63067的介绍可以知道，对应NASH和normal的样本的ID号，就可以提取我们需要的表达矩阵

把前面两属于Steatosis的样本去掉即可，exprSet=exprSet[,-c(1:2)]

然后再把芯片探针的id转换成entrez id

exprSet=exprSet[,-c(1:2)]

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

probeset <- rownames(exprSet)

rowMeans <- rowMeans(exprSet)

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

match_row=aggregate(rowMeans,by=list(EGID),max)

colnames(match_row)=c("EGID","rowMeans")

dat=data.frame(EGID,rowMeans,probeset)

tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))

relevantProbesets=as.character(tmp_prob$probeset)

length(relevantProbesets) #hgu133plus2.db  20156

exprSet=exprSet[relevantProbesets,]

EGID_name=as.numeric(lookUp(relevantProbesets, platformDB, "ENTREZID"))

rownames(exprSet)=as.character(EGID_name)

d=exprSet

最后得到表达矩阵表格

image002

我们首先得到1000次随机挑选128个基因的表达矩阵的主成分分析，第一主成分贡献度值。

gene128=sapply(1:1000,function(y) {

dat=t(d[sample(row.names(d), 128, replace=TRUE), ]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

}

)

很快就能得到结果，可以看到数据如下

>  summary(gene128)

Min. 1st Qu.  Median    Mean 3rd Qu.    Max.

19.1    25.8    28.8    29.8    32.5    59.7

image003

 

那么接下来我们挑选这个04380这条pathway特有128个基因来算第一主成分贡献度值

path_04380_gene=intersect(rownames(d),as.character(Path2GeneID[['04380']]))

dat=t(d[path_04380_gene,]);

round(100*summary(fast.prcomp(dat))$importance[2,1],2)

image004

得到的值是38.83，然后看看我们的这个38.83在之前随机得到的1000个数里面是否正常，就按照正态分布检验来算

1-pnorm((38.83-mean(gene128))/sd(gene128))

[1] 0.0625

可以看到已经非常显著的不正常了，可以说明这条通路被富集了。

至少说明超几何分布检验的方法得到的富集分析结果跟我们这次的重抽样+主成分分析结果是一致的，当然，也有不一致的，不然就不用发明一种新的方法了。

如果写一个循环同样可以检验所有的通路，但是这样就不需要事先准备好差异基因啦！！！这是这个分析方法的特点！

十二 21

主成分分析略讲

Posted on by

主成分分析是为了简化变量的个数。
我这里不涉及到任何高级统计知识来简单讲解一下主成分分析，首先我们用下面的代码随机创造一个矩阵：

options(digits = 2)
x=c(rnorm(5),rnorm(5)+4)
y=3*c(rnorm(5),rnorm(5)+4)
dat=rbind(x,y,a=0.1*x,b=0.2*x,c=0.3*x,o=0.1*y,p=0.2*y,q=0.3*y)
colnames(dat)=paste('s',1:10,sep="")
dat
library(gmodels)
pca=fast.prcomp(t(dat))
pca
summary(pca)$importance
biplot(pca, cex=c(1.3, 1.2));

Continue reading →

十二 15

用超几何分布检验做富集分析

Posted on by

我们可以直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

我们直接读取用limma包做好的差异分析结果

setwd("D:\\my_tutorial\\补\\用limma包对芯片数据做差异分析")

DEG=read.table("GSE63067.diffexp.NASH-normal.txt",stringsAsFactors = F)

View(DEG)

image001

我们挑选logFC的绝对值大于0.5，并且P-value小雨0.05的基因作为差异基因，并且转换成entrezID

probeset=rownames(DEG[abs(DEG[,1])>0.5 & DEG[,4]<0.05,])

library(hgu133plus2.db)

library(annotate)

platformDB="hgu133plus2.db";

EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

length(unique(EGID))

#[1] 775

diff_gene_list <- unique(EGID)

这样我们的到来775个差异基因的一个list

首先我们直接使用R的bioconductor里面的一个包，GOstats里面的函数来做超几何分布检验，看看每条pathway是否会富集

library(GOstats)

library(org.Hs.eg.db)

#then do kegg pathway enrichment !

hyperG.params = new("KEGGHyperGParams", geneIds=diff_gene_list, universeGeneIds=NULL, annotation="org.Hs.eg.db",

categoryName="KEGG", pvalueCutoff=1, testDirection = "over")

KEGG.hyperG.results = hyperGTest(hyperG.params);

htmlReport(KEGG.hyperG.results, file="kegg.enrichment.html", summary.args=list("htmlLinks"=TRUE))

结果如下：

image002

但是这样我们就忽略了其中的原理，我们不知道这些数据是如何算出来的，只是由别人写好的包得到了结果罢了。

事实上，这个包的这个hyperGTest函数无法就是包装了一个超几何分布检验而已。

如果我们了解了其中的统计学原理，我们完全可以写成一个自建的函数来实现同样的功能。

超几何分布很简单，球分成黑白两色，数量已知，那么你随机抽有限个球，应该抽多少白球的问题！
公式就是 exp_count=n*M/N
然后你实际上抽了多少白球，就可以计算一个概率值！
换算成通路的富集概念就是，总共有多少基因，你的通路有多少基因，你的通路被抽中了多少基因（在差异基因里面属于你的通路的基因），这样的数据就足够算出上面表格里面所有的数据啦！
tmp=toTable(org.Hs.egPATH)
GeneID2Path=tapply(tmp[,2],as.factor(tmp[,1]),function(x) x)
Path2GeneID=tapply(tmp[,1],as.factor(tmp[,2]),function(x) x)
#phyper(k-1,M, N-M, n, lower.tail=F)
#n*M/N
diff_gene_has_path=intersect(diff_gene_list,names(GeneID2Path))
n=length(diff_gene_has_path) #321 # 这里算出你总共抽取了多少个球
N=length(GeneID2Path) #5870  ##这里算出你总共有多少个球（这里是错的，有多少个球取决于背景基因！一般是两万个）
options(digits = 4)
for (i in names(Path2GeneID)){
 M=length(Path2GeneID[[i]])  ##这个算出你的所有的球里面，白球有多少个
 exp_count=n*M/N  ###这里算出你抽取的球里面应该多多少个是白色
 k=0         ##这个k是你实际上抽取了多少个白球
 for (j in diff_gene_has_path){
 if (i %in% GeneID2Path[[j]]) k=k+1
 }
 OddsRatio=k/exp_count
 p=phyper(k-1,M, N-M, n, lower.tail=F)  ##根据你实际上抽取的白球个数，就能算出富集概率啦！
 print(paste(i,p,OddsRatio,exp_count,k,M,sep="    "))
}
随便检查一下，就知道结果是一模一样的！