生信代码分享，有手就行，列表基因内部相关性分析
#install.packages("corrplot")
#本代码需要准备的输入数据为如下：
#1.基因列表文件，无表头的，可以是FPKM或TPM或微阵列数据
#2.表达矩阵文件，第一行为样本名，第一列为基因名
library(limma)
library(corrplot)
setwd("C:\\MetaCode\\1.Bioinformatics.Correlation.Analysis") #设置工作目录，注意是双斜杠，不可以出现空格、中文。
rt=read.table("表达矩阵.txt",sep="\t",header=T,check.names=F)#读取表达矩阵文件，被命名为rt，它是有表头的，按\t分隔，不对它的名字进行检查，如果是TCGA的数据，没有“check.names=F”，那么“-”会被转为“.”
rt=as.matrix(rt)#将rt转为矩阵
rownames(rt)=rt[,1]#rt的第一列是行名
exp=rt[,2:ncol(rt)]#rt第二列开始到最后一列是表达量
dimnames=list(rownames(exp),colnames(exp))#在构建矩阵的时候确定矩阵的行名和列名
data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)#在构建矩阵的时候确定矩阵的行名和列名，此时我们的矩阵文件被命名为data文件
data=avereps(data)#如果存在基因名重复的话，对重复的基因取平均值。
#data=data[rowMeans(data)>0,]#R 语言中的 rowMeans() 函数用于找出 DataFrame 、矩阵或数组的每一行的均值，也就是说，通过这行代码可以去除表达均值为负数的基因。在现实应用中，我们往往去除低表达的基因，因为低表达的基因可能是干扰、噪音，特别是在转为TPM格式且与微阵列数据进行combat算法的合并的时候，因为GEO的数据往往偏大。
data=log2(data+1)#进行log2标准化，可以将数据拉到一定的范围，为了避免出现表达量为0导致报错的情况，一般统一加个1。数据的标准化(normalization)是将数据按比例缩放，使之落入一个小的特定区间。在某些比较和评价的指标处理中经常会用到，去除数据的单位限制，将其转化为无量纲的纯数值，便于不同单位或量级的指标能够进行比较和加权。标准一般化不会影响你进行相关性和差异分析。
gene=read.table("gene.txt", header=F, check.names=F, sep="\t")#获取基因表达量，命名为gene文件
sameGene=intersect(as.vector(gene[,1]),rownames(data))#as.vector函数可以将矩阵进行向量化，此代码的含义是，将gene文件第一列与data文件的行名进行取交集，交集基因命名为sameGene
geneExp=data[sameGene,]#对data文件的每一行进行比对，以提取data文件中含有sameGene的行，即我们目标基因的表达矩阵
out=rbind(ID=colnames(geneExp),geneExp)#利用函数cbind()和rbind() 把向量或矩阵拼成一个新的矩阵：cbind()把矩阵横向合并成一个大矩阵（列方式），而rbind()是纵向合并（行方式）。在这里就是把目标基因的表达矩阵和geneExp的列名进行合并，并定义为out文件。
write.table(out,file="geneExp.txt",sep="\t",quote=F,col.names=F)#输出out文件以便后续分析，它的名字是"geneExp.txt"，是按sep="\t",不需要使用quote进行指定引号
rt=read.table("geneExp.txt",sep="\t",header=T,row.names=1,check.names=F)#读取目标基因的表达矩阵文件，被命名为rt，它是有表头的，按\t分隔，不对它的名字进行检查。
rt=t(rt)#对rt文件进行转置，类似于excel的行列倒置
res1 <- cor.mtest(rt, conf.level = 0.95)#著性检验，为每个值生成 p 值和置信区间 输入要素对
pdf("correlation.pdf",height=16,width=16)#保存图片的文件名称
corrplot(corr=cor(rt),
method = "circle",
order = "hclust",
tl.col="black",
addCoef.col = "black",
p.mat = res1$p,
sig.level = 0.001,
insig = "pch",
number.cex = 1,
type = "upper",
col=colorRampPalette(c("blue", "white", "red"))(50),
)
dev.off()
