我是一個(gè)很菜的鳥藏雏，為了計(jì)算一個(gè)基因與其他基因的相關(guān)分析，在網(wǎng)上找了很多資料作煌！我看不懂腳本诉稍，就只能一個(gè)個(gè)死看，一個(gè)個(gè)嘗試最疆！最后杯巨，柳暗花明，找到這個(gè)朋友推送的方法努酸。雖然我不懂他寫了什么服爷，也不懂每一步到底什么意思，但跟著這個(gè)朋友的分享的腳本跑了出來获诈！非常感謝他Ｈ栽础！（留念第一篇簡書舔涎，希望大家批評指正笼踩！也希望給需要的朋友帶來幫助！Ｍ鱿印）

備注：1嚎于、我的數(shù)據(jù)是自然群體的FPKM矩陣，橫向?yàn)榛虮磉_(dá)量挟冠，縱列為樣本于购，第一行是樣本名，第一列為基因名知染；

2肋僧、以下唯一要改的是將“gene A”修改成自己的目標(biāo)基因，隨后就是計(jì)算其他所有基因與目標(biāo)基因的相關(guān)系數(shù)及p值？氐嫌吠！

https://mp.weixin.qq.com/s?src=11&timestamp=1675943258&ver=4340&signature=L2Ybm0h5bdumwp798BLbTQOi655Ehu0gG7ck5sDUA39ZpyxwVHvtxVneSc9EJYrvTme5-rYJdT-95bsj7lMByZlKLDpipE7DL1GyY5y2udZaDj3-T2ergQsY-bT08G0H&new=1

> rm(list = ls())

> options(stringsAsFactors = F)

> tcga_foxo1 <- read.table(file ="gwas.txt",header = T, sep = "\t", quote = "", check.names = F)

> class(tcga_foxo1)

[1] "data.frame"

> row.names(tcga_foxo1) <- tcga_foxo1[,1]

> exprSet1 <- tcga_foxo1 [,-1]

> test1 <- exprSet1[1:10,1:10]

> test1

? ? ? ? ? ? ? ? ?B1? ? ? ?B10? ? ? ? B11? ? ? ?B12? ? ? B13? ? ? ?B14

LG0101766 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101798 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101802 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101812 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101832 0.3205862 3.6745567 0.67330880 0.5555175 1.100432 0.4851674

LG0101836 0.1826809 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101844 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101875 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

LG0101909 0.1714331 0.0963218 0.07501063 0.0825173 0.000000 0.0000000

LG0101924 0.0000000 0.0000000 0.00000000 0.0000000 0.000000 0.0000000

? ? ? ? ? ? ? ? ?B15? ? ? ? B16? ? ? B17? ? ? ? B18

LG0101766 0.00000000 0.00000000 0.000000 0.00000000

LG0101798 0.00000000 0.00000000 0.000000 0.00000000

LG0101802 0.00000000 0.00000000 0.000000 0.00000000

LG0101812 0.00000000 0.00000000 0.000000 0.00000000

LG0101832 0.19367102 0.18983320 1.517738 0.76362680

LG0101836 0.18393392 0.00000000 0.000000 0.77703675

LG0101844 0.09542112 0.00000000 0.000000 0.00000000

LG0101875 0.00000000 0.00000000 0.000000 0.00000000

LG0101909 0.25891335 0.08459422 0.000000 0.07291937

LG0101924 0.00000000 0.00000000 0.000000 0.00000000

> exprSet2<- t(exprSet1 )

> class(exprSet1 )

[1] "data.frame"

> class(exprSet2)

[1] "matrix" "array"?

> exprSet <- as.data.frame(exprSet2)

> class(exprSet)

[1] "data.frame"

> test2 <- exprSet[1:10,1:10]

> test2

? ? LG0101766 LG0101798 LG0101802 LG0101812 LG0101832 LG0101836? LG0101844

B1? ? ? ? ? 0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.3205862 0.1826809 0.00000000

B10? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 3.6745567 0.0000000 0.00000000

B11? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.6733088 0.0000000 0.00000000

B12? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.5555175 0.0000000 0.00000000

B13? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 1.1004324 0.0000000 0.00000000

B14? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.4851674 0.0000000 0.00000000

B15? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.1936710 0.1839339 0.09542112

B16? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.1898332 0.0000000 0.00000000

B17? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 1.5177376 0.0000000 0.00000000

B18? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0? ? ? ? ?0 0.7636268 0.7770368 0.00000000

? ? LG0101875? LG0101909 LG0101924

B1? ? ? ? ? 0 0.17143306? ? ? ? ?0

B10? ? ? ? ?0 0.09632180? ? ? ? ?0

B11? ? ? ? ?0 0.07501063? ? ? ? ?0

B12? ? ? ? ?0 0.08251730? ? ? ? ?0

B13? ? ? ? ?0 0.00000000? ? ? ? ?0

B14? ? ? ? ?0 0.00000000? ? ? ? ?0

B15? ? ? ? ?0 0.25891335? ? ? ? ?0

B16? ? ? ? ?0 0.08459422? ? ? ? ?0

B17? ? ? ? ?0 0.00000000? ? ? ? ?0

B18? ? ? ? ?0 0.07291937? ? ? ? ?0

> y <-as.numeric(exprSet[,"gene A"])

> colnames <- colnames(exprSet)

> cor_data_df <- data.frame(colnames)

> for (i in 1:length(colnames)){

+ test <- cor.test(as.numeric(exprSet[,i]),y,type="spearman")

+?

+? cor_data_df[i,2] <- test$estimate

+?

+? cor_data_df[i,3] <- test$p.value

+?

+ }

There were 50 or more warnings (use warnings() to see the first 50) #忽略不用管這個(gè)

> names(cor_data_df) <- c("symbol","correlation","pvalue")

> head(cor_data_df)

? ? ?symbol? correlation? ? ?pvalue

1 LG0101766? ? ? ? ? ?NA? ? ? ? ?NA

2 LG0101798? ? ? ? ? ?NA? ? ? ? ?NA

3 LG0101802? 0.108265485 0.06947002

4 LG0101812? 0.006733649 0.91036654

5 LG0101832 -0.144771832 0.01496744

6 LG0101836? 0.004308967 0.94257104

> write.table(cor_data_df, file = "cor_data_df", sep = "\t",row.names = TRUE,col.names = NA)