R语言相关性分析

Python09

R语言相关性分析,第1张

1.  R语言自带函数cor(data, method=" ")可以快速计算出相关系数 ,数据类型:data.frame

 如data.frame为:zz, 绘图如下:

a. single protein:线性回归画法

1. ggplot(zz,aes(x=a, y=HDL))+

   geom_point(alpha=1,colour="#FFA54F")+

   geom_smooth(method = lm,colour="#8B658B")+

   #scale_color_brewer(palette = "Set1")+

   theme_bw()+

   labs(x="Ferritin",y="HDL.C",title="Pearson’s correlation test of ferritin and HDL.C")+

   annotate("text", x = 1000, y = 2.5, label = "r = -0.51",colour="black",size=4)

2. library(ggstatsplot)

 ggscatterstats(data = alldata,

               y = TRANSFUSION.UNIT,

                x = NPTXR,

                centrality.para = "mean",  #"mean" or "median"                         

               margins = "both",                                       

                xfill = "#D8BFD8",

                yfill = "#EEDD82",

                #line.size= ,

                line.color="#8B6969",

               point.color="#2F4F4F",

                marginal.size=4,

               marginal.type = "density", # "histogram", "boxplot", "density", "violin", "densigram")

                title = "Relationship between TRANSFUSION.UNIT and NPTXR")

b. ggcorrplot, 全部蛋白 global correlation map 画法

ggcorrplot(cor(alldata))

2.  summary(lm(y~x),method=" ") %>%.[["coefficients"]]   正规线性回归

     (其实就是:a<-lm(y~x1+x2+...,data)

      plot(summary(lm(y~x),method=" ")) #绘图

3.  ggcor部分数据绘图:  数据类型为data.frame,纵坐标为各指标or各蛋白,行为观测值。

data <- fortify_cor(alldata[,10:11],alldata,cluster.type = "col")

ggcor<-ggcor(data,label_size=0.5) +

  geom_colour()+

  theme(axis.text.x = element_text(colour = "black",size = 4.7),

                                                        axis.text.y=element_text(size=5.5),

                                                        axis.ticks=element_blank())+

  geom_num(aes(num=r),colour="black",size=1.5)

4. corrr包画法

datasets::mtcars %>%

  correlate() %>%

  focus(-cyl, -vs, mirror = TRUE) %>%

  rearrange() %>%

  network_plot(min_cor = .2)

给你一些代码,你慢慢研究:install.packages('ggplot2')library(ggplot2)ggplot(a)+geom_bar(aes(x1,y,fill/col=x1/x2),position='dodge',stat='summary',fun='sum'/'mean')条形图+theme(text = element_text(family='Kai'))ggplot(a)+geom_boxplot(aes(x1,y,col=x1/x2))箱线图ggplot(a)+geom_point(aes(x1,y,col=x1/x2),position=position_jitter(width=0.04))散点图1+geom_point(aes(x1,y,col=x1/x2),stat='summary',fun='sum'/'mean')+散点2+geom_line(aes(x1,y,group=1/x2,col=x1/x2),stat='summary',fun='sum'/'mean')+折线3+geom_errorbar(aes(x=x1,ymin=y-se,ymax=y+se,col=x1/x2),position=position_dodge(0.9),width=0.2)+误差棒4+geom_text(aes(x1,y,label=marker,col=x1/x2),position=position_dodge(0.9)vjust=2或y+2)+显著字母ggplot(a,aes(x1,y,fill/col=x1/x2))+geom_bar(position='dodge',stat='summary',fun='sum'/'mean')+geom_errorbar(aes(ymin=y-se,ymax=y+se),position=position_dodge(0.9),width=0.2)+geom_text(aes(label=marker),position=position_dodge(0.9),vjust=-2)条形图+误差棒+显著字母(坐标写一次即可)ggplot(a,aes(x1,y,col=x1/x2))+geom_point(position=position_jitter(width=0.04),stat='summary',fun='sum'/'mean')+geom_line(aes(group=1/x2),stat='summary',fun='sum'/'mean')+geom_errorbar(aes(ymin=y-se,ymax=y+se),position=position_dodge(0.9),width=0.2)+geom_text(aes(label=marker),position=position_dodge(0.9),vjust=-2)散点图+折线+误差棒+显著字母(坐标写一次即可)+geom_density(aes(y=liqi))密度图(1个数值型)+geom_area(aes(x=tan,y=liqi))区域图(2个数值型)+geom_smooth(aes(x=tan,y=liqi,group/col=chong),formula=y~x,method='lm',se=F)拟合图,分组/线条颜色(2个数值型)+facet_wrap(~riqi,ncol/nrow=2,labeller='label_both/value')分面图,每行或每列分面数,分面标题+xlab('自变量1(单位)')+ylab('因变量(单位)')+scale_fill_discrete(name='自变量2')更改轴和图例名称+coord_cartesian(ylim= c(0,80))限定轴范围(fill=x1/x2,有此即可变色)+scale_fill_manual(values = c('grey70', 'grey50', 'grey30'))改变条形填充颜色(颜色数量=分组数量)(col=x1/x2,有此即可变色)+scale_color_manual(values = c('red', 'orange', 'yellow'))改变颜色(颜色数量=分组数量)

计算β多样性指数需要用到phyloseq包。它的安装方式不同于简单的install.packages(“phyloseq”)

有两种方法可以安装

1.先安装BiocManager

install.packages("BiocManager")

library("BiocManager")

BiocManager::install("phyloseq")

library("phyloseq")

2.source("https://bioconductor.org/biocLite.R")

biocLite("phyloseq")

#安装phyloseq

library("phyloseq")

安装并加载了phyloseq包后,开始读取数据,前面计算α多样性,用到的是read.table……

qiimedata <- import_qiime(otufilename = "feature-table.taxonomy.txt", mapfilename = "mapping_file.txt", treefilename = "tree.rooted.nwk", refseqfilename = "dna-sequences.fasta")

#读取数据,参数都是文件名,注意加后缀

#otufilename指定out表格,mapfilename指定map文件(分组数据)

#treefilename指定有根进化树文件

#refseqfilename指定代表序列文件

otu<-qiimedata@[email protected]

#从qiimedata里面提取otu

sum_of_otus<-colSums(t(otu))

#t_转置,colsums计算列的和,即计算各个otu检测到的总序列数,为了筛掉一些总序列数过低的otu(可能是测序错误)

sum_of_otus

#查看otu总序列数

selected_otu<-names(sum_of_otus)[sum_of_otus>10]

#获取总序列数大于10的otu id

sub_qiimedata <- prune_taxa(selected_otu, qiimedata)

#筛选总序列数大于10的otu的phyloseq数据

weighted_unifrac<-distance(sub_qiimedata,method = 'wunifrac')

#计算样本间加权unifrac

unweighted_unifrac<-distance(sub_qiimedata,method = 'unifrac')

#计算样本间非加权unifrac

bray_curtis <- distance(sub_qiimedata, method='bray')

write.table(as.matrix(bray_curtis),"bray_curtis.txt",sep = '\t',quote = FALSE,col.names = NA)

#保存距离矩阵

#计算样本间Bray-Curtis距离矩阵,method 可选" wunifrac ", " unifrac " ,"jaccard"等

pcoa_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "PCoA")

#基于Bray-Curtis距离矩阵的PCoA排序分析

p<-plot_ordination(sub_qiimedata, pcoa_of_bray_curtis, type="samples", color="Group1",shape = "Group1")

#将PCoA排序分析结果可视化

library("ggplot2")

p<-p+ scale_colour_manual(values=c("#DC143C","#808000","#00CED1")) + geom_point(size=2) +ggtitle("PCoA of Bray-Curtis distance")+theme(text = element_text(size = 15))

#修改图形大小,ggtitle加标题,stat_ellipse加椭圆

#用scale_colour_manual(values=c())自定义颜色,可查颜色的16进制对照表

p

nmds_of_bray_curtis<-ordinate(physeq=sub_qiimedata,distance = 'bray',method = "NMDS")

#基于Bray-Curtis距离矩阵的NMDS排序分析

p1<-plot_ordination(qiimedata, nmds_of_bray_curtis, type="samples", color="Group1")

#将NMDS排序分析结果可视化

# color=“Group1”指定不同分组的点染不同颜色

p1

p1<-p1+ geom_point(size=3) +ggtitle("NMDS of Bray-Curtis distance") + stat_ellipse()+theme(text = element_text(size = 15))

#对图片进行适当修饰, stat_ellipse()加椭圆, ggtitle()加标题

ggsave(plot = p1,“nmds_of_bary_curtis.pdf",dpi = 300,width

PCoA中的两个点距离,接近β多样性指数

PCA(Principal Components Analysis)即主成分分析,也称主分量分析或主成分回归分析法,首先利用线性变换,将数据变换到一个新的坐标系统中然后再利用降维的思想,使得任何数据投影的第一大方差在第一个坐标(称为第一主成分)上,第二大方差在第二个坐标(第二主成分)上。这种降维的思想首先减少数据集的维数,同时还保持数据集的对方差贡献最大的特征,最终使数据直观呈现在二维坐标系。

PCoA(Principal Co-ordinates Analysis)分析即主坐标分析,可呈现研究数据相似性或差异性的可视化坐标,是一种非约束性的数据降维分析方法,可用来研究样本群落组成的相似性或相异性。它与PCA类似,通过一系列的特征值和特征向量进行排序后,选择主要排在前几位的特征值,找到距离矩阵中最主要的坐标,结果是数据矩阵的一个旋转,它没有改变样本点之间的相互位置关系,只是改变了坐标系统。两者的区别为PCA是基于样本的相似系数矩阵(如欧式距离)来寻找主成分,而PCoA是基于距离矩阵(欧式距离以外的其他距离)来寻找主坐标。

NMDS图中两个点的距离的排序,接近β多样性指数的排序