「Workshop」第四十期 常用的差异分析方法

实际应用中，当数据量足够大时，对样本正态分布要求不再严格。只要数据分布不是严重偏态，一般来说单样本t检验都是适用的。

#生成数据
x <- c(130,120,130,110,130,135,129,124,130,134)
#t检验
t.test(x,mu = 130)

 One Sample t-test

data:  x
t = -1.1884, df = 9, p-value =
0.2651
alternative hypothesis: true mean is not equal to 130
95 percent confidence interval:
 121.8702 132.5298
sample estimates:
mean of x 
    127.2 
#结果显示，P=0.2651＞0.05。在统计学上说明样本均数与总体均数没有差别。

A <- c(85,84,95,73,77,65,85,93,90,91)
B <- c(87,96,77,80,79,96,93,82,84,86)
#方差齐性检验
#合并数据
y <- c(A,B)
#数据分组标签
group=as.factor(c(rep(1,10),rep(2,10)))
#载入car包
library(car)
#方差齐性检验
leveneTest(y=y,group = group)
#结果显示，P=0.5505＞0.05。说明方差齐性。
Levene's Test for Homogeneity of Variance (center = median)
      Df F value Pr(>F)
group  1  0.3703 0.5505
      18  
#独立样本t检验
t.test(A,B,paired = FALSE)
#结果显示P=0.5639。说明两者没有区别。
 Welch Two Sample t-test

data:  A and B
t = -0.589, df = 16.463, p-value = 0.5639
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -10.100024   5.700024
sample estimates:
mean of x mean of y 
     83.8      86.0 

#生成数据
drug1 <- c(4.4,5,5.8,4.6,4.9,4.8,6,5.9,4.3,5.1)
drug2 <- c(6.2,5.2,5.5,5,4.4,5.4,5,6.4,5.8,6.2)
#配对样本t检验
t.test(drug1,drug2,paired = TRUE)
#结果显示，P=0.1575＞0.05,不能说两者存在显著差别。
 Paired t-test

data:  drug1 and drug2
t = -1.5417, df = 9, p-value = 0.1575
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -1.0609306  0.2009306
sample estimates:
mean of the differences 
                  -0.43 

fit <- aov(yield ~ N, data=npk)
fit <- aov(yield ~ N + block, data=npk)

x<-c(210,312,170,85,223)
chisq.test(x)

x<-matrix(c(46,18,6,8),ncol=2,nrow=2)
chisq.test(x)
chisq.test(x)$expected ###查看理论值
fisher.test(x)  ##进行fisher检验

library(stats)
data("mtcars")
wilcox.test(mpg~am,data = mtcars)

          1.limma包做差异分析要求数据满足正态分布或近似正态分布，如基因芯片、TPM格式的高通量测序数据。
          2.通常认为Count数据不符合正态分布而服从泊松分布。对于count数据来说，用limma包做差异分析，误差较大
          3.DESeq2、和 EdgeR都是基于count，然后两个都是NB（negative binomial）但是在估计dispersion parameter的方法上面不一样。
          4.limma，edgeR，DESeq2三大包基本是做转录组差异分析的金标准，大多数转录组的文章都是用这三个R包进行差异分析。
          5.edgeR差异分析速度快，得到的基因数目比较多，假阳性高（实际不差异，结果差异）。DESeq2差异分析速度慢，得到的基因数目比较少，假阴性高（实际差异，结果不差异）。
          6.需要注意的是制作分组信息的因子向量是，因子水平的前后顺序，在R的很多模型中，默认将因子向量的第一个水平看作对照组

library("edgeR")
expr = read.csv("mRNA_exprSet.csv",sep = ',',header=T)  
head(expr)
expr = avereps(expr[,-1],ID = expr$X) # 对重复基因名取平均表达量，然后将基因名作为行名
expr = expr[rowMeans(expr)>1,] # 去除低表达的基因
library(stringr)
tumor <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) < 10]
normal <- colnames(expr)[as.integer(substr(colnames(expr),14,15)) >= 10]

tumor_sample <- expr[,tumor]
normal_sample <- expr[,normal]

exprSet_by_group <- cbind(tumor_sample,normal_sample)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))

data = exprSet_by_group
group_list = factor(group_list)
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)
DGElist <- DGEList( counts = data, group = group_list)
## Counts per Million or Reads per Kilobase per Million
keep_gene <- rowSums( cpm(DGElist) > 1 ) >= 2 ## 自定义
table(keep_gene)
DGElist <- DGElist[ keep_gene, , keep.lib.sizes = FALSE ]

DGElist <- calcNormFactors( DGElist )
DGElist <- estimateGLMCommonDisp(DGElist, design)##为所有基因都计算同样的dispersion
DGElist <- estimateGLMTrendedDisp(DGElist, design)##根据不同基因的均值--方差之间的关系来计算dispersion，并拟合一个trended model
DGElist <- estimateGLMTagwiseDisp(DGElist, design)##计算每个基因的dispersion，并利用经验性贝叶斯方法（empirical bayes）压缩至Trend Didspersion。个人认为这一项相当于GLM中每个基因的beta值
#####################负二项式广义对数线性模型
fit <- glmFit(DGElist, design)
results <- glmLRT(fit, contrast = c(-1, 1)) 
####################类似然负二项式广义对数线性模型
fit <- glmQLFit(dge, design, robust = TRUE)        #拟合模型
lrt <- glmQLFTest(fit) 

nrDEG_edgeR <- topTags(results, n = nrow(DGElist))
nrDEG_edgeR <- as.data.frame(nrDEG_edgeR)
head(nrDEG_edgeR)

padj = 0.01 # 自定义
foldChange= 2 # 自定义
nrDEG_edgeR_signif  = nrDEG_edgeR[(nrDEG_edgeR$FDR < padj & 
                       (nrDEG_edgeR$logFC>foldChange | nrDEG_edgeR$logFC<(-foldChange))),]
nrDEG_edgeR_signif = nrDEG_edgeR_signif[order(nrDEG_edgeR_signif$logFC),]

library(edgeR)
##跟DESeq2一样，导入数据，预处理（用了cpm函数）
exprSet<- read.table(file = "test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep("Contral",1),rep("Treat",1)))

##设置分组信息，并做TMM标准化
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.4  #设置BCV值
et <- exactTest(exprSet, dispersion=bcv^2)
write.csv(topTags(et, n = nrow(exprSet$counts)), 'result.csv', quote = FALSE)   #输出主要结果