library(limma)
mRNA <- read.table("表达矩阵.txt",sep = "\t",header = T,comment.char = "!",encoding = "UTF-8")
group <- read.table("GROUP.txt",header=T,sep = "\t",encoding = "UTF-8")
group_CP <- group$treat
m_design<- model.matrix(~0+factor(group_CP))
colnames(m_design) = levels(factor(group_CP))
rownames(m_design)= group$ID
contrast.matrix<-makeContrasts("P-C",levels=m_design) #注意P-C顺序，实验组要在前面否则影响上下调结果
m_fit <- lmFit(mRNA,m_design)
m_fit <- contrasts.fit(m_fit, contrast.matrix)
m_fit <- eBayes(m_fit)
m_genlist <- topTable(m_fit, coef = 1, n=Inf) #limma结果
#将表达矩阵与差异分析结果合并
ID_REF <- rownames(m_genlist)
m_genlist <- data.frame(ID_REF,m_genlist)
ID_REF <- rownames(mRNA)
mRNA <- data.frame(ID_REF,mRNA)
test <-merge(mRNA,m_genlist,by = "ID_REF")
result <- subset(test,P.Value<0.05)
row.names(result) <- result[,1]
#绘制热图
heatmap <- result[2:(nrow(group)+1)]
annotation <- data.frame(Factor = factor(group$treat)) #标注样本的分组信息
rownames(annotation) <- colnames(heatmap)
library(pheatmap)
filename <- paste("文件名",".pdf",sep="")
pdf(filename)
pheatmap(heatmap,
annotation=annotation,
annotation_legend = TRUE,
main=filename ,
scale = "row",
show_rownames = F,
color = colorRampPalette(c("green","black","red"))(100))
dev.off()
#表达矩阵与GROUP文件如下所示
1 读取，计算均值，箱图观察
2 查看数据分布
2.1 hist直方图
2.2 qqnorm散点图
3 Shapiro-Wilk正态性检验
4 方差齐性检验
意义：方差分析就是在大家误差水平差不多的条件下看控制和对照组是不是有显著差异。那方差其实就是误差水平了。当方差不一致的时候，这个方法就没法分辨出究竟是控制造成的差异还是，内在的波动造成的差异。
参考： https://www.zhihu.com/question/21195390
参考： https://blog.csdn.net/tiaaaaa/article/details/58130363
4.1 F检验
使用条件：数据正态分布，只可以检验两个样本
4.2 bartlett检验
使用条件：正态分布的数据，多个样本
4.3 levene检验
没有条件：数据可不具正态性，可以检验多个总体的方差齐性
SPSS的默认方差齐性检验方法
5 差异检验
5.1 参数检验：T检验
使用条件：两样本来自正太分布总体，方差齐
5.2 非参数检验：Wilcoxon秩和检验（两样本）
参数：
参考： https://www.jianshu.com/p/f30d1fe877ea
5.3 非参数检验：Kruskal-Wallis（KS）秩和检验（多样本）
5.4 Deseq两组reads count差异分析
#PCoA 分析在R语言中进行主要依赖于以下得包，进行这个分析得主要可以应用于形态学数据得相似与差异性分析。library(ade4)
library(ggplot2)
library(RColorBrewer)
library(vegan)
这里我们使用R自带得数据iris
data(iris)
在R语言中通常都会使用这个数据进行案例分析
#iris
data(iris)
iris
data01<-iris[,-5]#数据预处理，去掉最后一列得数据标签
data01
dis01<-vegdist(data01,method = "euclidean")#这里是为了算矩阵距离，方法根据数据选择合适得方法
dis01
pcoa1<- dudi.pco(dis01, scan = FALSE,nf=3)#进行PCoA分析
pcoa1
pcoa1_eig<-pcoa1$eig[1:2]/sum(pcoa1$eig)#算一下前两列对整个数据得解释比例
pcoa1_eig
samplesite1<-data.frame({pcoa1$li})[1:2]#将前两列的数据分析结果放到sample_site1里面
sample_site1
sample_site1$names<-rownames(sample_site1)#设置名称
sample_site1$names
iris$Species
sample_site1$level<-factor(iris$Species,levels = c("setosa","versicolor","virginica"))#设置level的标签
sample_site1$level
names(sample_site1)[1:2]<-c("PCoA1","PCoA2")
p<-ggplot(sample_site1, mapping=aes(PCoA1, PCoA2,color=level))+theme_classic()
p<-p+geom_point()#绘制散点图
p
