#install.packages("corrplot")

#本代碼需要準(zhǔn)備的輸入數(shù)據(jù)為如下：

#1.基因列表文件，無(wú)表頭的，可以是FPKM或TPM或微陣列數(shù)據(jù)

#2.表達(dá)矩陣文件，第一行為樣本名，第一列為基因名

library(limma)

library(corrplot)

setwd("C:\\MetaCode\\1.Bioinformatics.Correlation.Analysis")?? #設(shè)置工作目錄，注意是雙斜杠，不可以出現(xiàn)空格、中文。

rt=read.table("表達(dá)矩陣.txt",sep="\t",header=T,check.names=F)#讀取表達(dá)矩陣文件，被命名為rt，它是有表頭的，按\t分隔，不對(duì)它的名字進(jìn)行檢查，如果是TCGA的數(shù)據(jù)，沒(méi)有“check.names=F”，那么“-”會(huì)被轉(zhuǎn)為“.”

rt=as.matrix(rt)#將rt轉(zhuǎn)為矩陣

rownames(rt)=rt[,1]#rt的第一列是行名

exp=rt[,2:ncol(rt)]#rt第二列開(kāi)始到最后一列是表達(dá)量

dimnames=list(rownames(exp),colnames(exp))#在構(gòu)建矩陣的時(shí)候確定矩陣的行名和列名

data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)#在構(gòu)建矩陣的時(shí)候確定矩陣的行名和列名，此時(shí)我們的矩陣文件被命名為data文件

data=avereps(data)#如果存在基因名重復(fù)的話，對(duì)重復(fù)的基因取平均值。

#data=data[rowMeans(data)>0,]#R 語(yǔ)言中的 rowMeans() 函數(shù)用于找出 DataFrame 、矩陣或數(shù)組的每一行的均值，也就是說(shuō)，通過(guò)這行代碼可以去除表達(dá)均值為負(fù)數(shù)的基因。在現(xiàn)實(shí)應(yīng)用中，我們往往去除低表達(dá)的基因，因?yàn)榈捅磉_(dá)的基因可能是干擾、噪音，特別是在轉(zhuǎn)為T(mén)PM格式且與微陣列數(shù)據(jù)進(jìn)行combat算法的合并的時(shí)候，因?yàn)镚EO的數(shù)據(jù)往往偏大。

data=log2(data+1)#進(jìn)行l(wèi)og2標(biāo)準(zhǔn)化，可以將數(shù)據(jù)拉到一定的范圍，為了避免出現(xiàn)表達(dá)量為0導(dǎo)致報(bào)錯(cuò)的情況，一般統(tǒng)一加個(gè)1。數(shù)據(jù)的標(biāo)準(zhǔn)化(normalization)是將數(shù)據(jù)按比例縮放，使之落入一個(gè)小的特定區(qū)間。在某些比較和評(píng)價(jià)的指標(biāo)處理中經(jīng)常會(huì)用到，去除數(shù)據(jù)的單位限制，將其轉(zhuǎn)化為無(wú)量綱的純數(shù)值，便于不同單位或量級(jí)的指標(biāo)能夠進(jìn)行比較和加權(quán)。標(biāo)準(zhǔn)一般化不會(huì)影響你進(jìn)行相關(guān)性和差異分析。

gene=read.table("gene.txt", header=F, check.names=F, sep="\t")#獲取基因表達(dá)量，命名為gene文件

sameGene=intersect(as.vector(gene[,1]),rownames(data))#as.vector函數(shù)可以將矩陣進(jìn)行向量化，此代碼的含義是，將gene文件第一列與data文件的行名進(jìn)行取交集，交集基因命名為sameGene

geneExp=data[sameGene,]#對(duì)data文件的每一行進(jìn)行比對(duì)，以提取data文件中含有sameGene的行，即我們目標(biāo)基因的表達(dá)矩陣

out=rbind(ID=colnames(geneExp),geneExp)#利用函數(shù)cbind()和rbind()?把向量或矩陣拼成一個(gè)新的矩陣：cbind()把矩陣橫向合并成一個(gè)大矩陣（列方式），而rbind()是縱向合并（行方式）。在這里就是把目標(biāo)基因的表達(dá)矩陣和geneExp的列名進(jìn)行合并，并定義為out文件。

write.table(out,file="geneExp.txt",sep="\t",quote=F,col.names=F)#輸出out文件以便后續(xù)分析，它的名字是"geneExp.txt"，是按sep="\t",不需要使用quote進(jìn)行指定引號(hào)

rt=read.table("geneExp.txt",sep="\t",header=T,row.names=1,check.names=F)#讀取目標(biāo)基因的表達(dá)矩陣文件，被命名為rt，它是有表頭的，按\t分隔，不對(duì)它的名字進(jìn)行檢查。

rt=t(rt)#對(duì)rt文件進(jìn)行轉(zhuǎn)置，類似于excel的行列倒置

res1 <- cor.mtest(rt, conf.level = 0.95)#著性檢驗(yàn)，為每個(gè)值生成 p 值和置信區(qū)間 輸入要素對(duì)

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()