lab2-DEG_heatmap-st

練習(xí)1：從PDAC-NET.txt表達(dá)譜中獲得PDAC和NET的差異基因，(顯著性水平為0.05 ，fdr 控制在0.01)并對(duì)上述差基因畫(huà)熱圖。

heatmap必備包

library(gplots)

以下內(nèi)容都是模仿老師的代碼，不過(guò)我的文件夾結(jié)構(gòu)和老師不完全相同，小改
讀入表達(dá)譜

gene_exp <- read.table("./data/PDAC-NET.txt", sep="\t", header=T,row.names=1)
gene <- as.matrix(gene_exp)

對(duì)輸出文件的創(chuàng)造（相比老師，這里把a(bǔ)ppend參數(shù)改為FALSE，為了后續(xù)電腦重啟繼續(xù)做，跑一遍的時(shí)候不會(huì)做出重復(fù)內(nèi)容的文件）

file.create("./result/PDAC-NET-fc-ttest-all.txt")
write.table("GeneName\tPDAC\tNET\tfold_change\tpvalue","./result/PDAC-NET-fc-ttest-all.txt",row.names=FALSE,col.names =FALSE, sep="\t",quote=FALSE,append=F)

獲取每個(gè)基因在樣本中的均值

mean_gene=rbind(apply(gene[,1:6],1,mean),apply(gene[,7:12],1,mean))
mean_gene=t(mean_gene)

P值

for (i in  1:length(rownames(gene))){# t-test pvalue
  gene_test=t.test(gene[i,1:6],gene[i,7:12])#t test for PDAC vs NET
  mean_ratio=2^( mean_gene[i,1]- mean_gene[i,2])   #  mean fold change
  gt=c(rownames(gene)[i],mean_gene[i,],as.vector(mean_ratio),gene_test[[3]])
  write.table(t(gt),"./result/PDAC-NET-fc-ttest-all.txt",row.names=FALSE,col.names =FALSE, sep="\t",quote=FALSE,append=TRUE)   
}

FDR矯正

result_test<-read.table("./result/PDAC-NET-fc-ttest-all.txt",sep="\t",header=TRUE,row.names=1)
p_value<-result_test[[4]] #此命令需要在as.matrix前面
result_test=as.matrix(result_test)
fdr<-p.adjust(p_value,method="BH",length(p_value))
result_fdr=cbind(result_test,fdr) #增加fdr列
test_fdr=result_fdr[result_fdr[,5]<(0.01)&result_fdr[,4]<(0.05),] 
#fdr控制在0.01，pvalue控制在0.05
col_names=colnames(result_test)
colnames( test_fdr)=c(col_names,"fdr")
write.table(test_fdr,"./result/PDAC-NET-test-fdr.txt",row.names=TRUE,col.names=NA, 
sep="\t",quote=FALSE)
#fdr控制在0.01后，得到的差異基因

畫(huà)熱圖

gene_exp <- read.table("./data/PDAC-NET.txt", sep="\t", header=T,row.names=1)
gene <- as.matrix(gene_exp)
diff_gene <- read.table("./result/PDAC-NET-test-fdr.txt",sep="\t",header=TRUE, row.names=1)
diff_gene=as.matrix(diff_gene) 
diff_order=diff_gene[order(diff_gene[,2]),]
diff_order=diff_order[1:120,]
hp_gene=gene[rownames(gene) %in% rownames(diff_order),] 
hp_gene1=cbind(hp_gene[,1:6],hp_gene[,7:12])
png("./result/PDAC-NET_heatmp.png",width = 800, height = 1000)
hmcols <- colorRampPalette(c("green","black","red"))(20)
heatmap(hp_gene1, col=hmcols,cexRow=1)
dev.off()
PDAC-NET_heatmp.png

選做：對(duì)火山圖（volcano plot）顏色重置：高表達(dá)紅色表示，低表達(dá)綠色表示，無(wú)差異表達(dá)灰色表示

2022.10.19 ↓成功達(dá)成顏色要求 順便學(xué)習(xí)了一下ggplot2的使用 畫(huà)了個(gè)輔助線

library(ggplot2)
library(gplots)
gene = read.table("./result/PDAC-NET-fc-ttest-all.txt", sep="\t", header=TRUE, row.names=1)
gene=as.data.frame(gene)
gene_list=gene[,c("fold_change","pvalue")]
gene_list[,"fold_change"]=log2(gene_list[,"fold_change"])
colnames(gene_list)=c("logFC","P.Value")
fcThreshold <- 2
pvalThreshold <- 0.001
gene_list$GeneSymbol = rownames(gene_list)
gene_list$expression = ifelse(gene_list$logFC >= fcThreshold & gene_list$P.Value < pvalThreshold,"Up-regulated",
                             ifelse(gene_list$logFC <=-fcThreshold & gene_list$P.Value < pvalThreshold,"Down-regulated","Unchanged"))
gene_list$log10pval = -log10(gene_list$P.Value)
gbasic=ggplot(data=gene_list,aes(x=logFC, y=log10pval))+
geom_point(size = 2,aes(color = expression))+
scale_colour_manual(values=c("Unchanged"="grey","Up-regulated"="red","Down-regulated"="green"))
scale_x_continuous(limits = c(-2.5, 2.5))+
geom_hline(yintercept=-log10(0.001),linetype=4)+
geom_vline(xintercept=c(-fcThreshold,fcThreshold),linetype=4)

png("./result/PDAC-NET-volcano.png",width = 800, height = 800)
print(gbasic)
dev.off()
PDAC-NET-volcano.png

2022.10.18 ↓ 畫(huà)的還不是特別滿足要求 可能之后會(huì)改吧

library(ggplot2)
library(ggrepel)
gene = read.table("./result/PDAC-NET-fc-ttest-all.txt", sep="\t", header=TRUE, row.names=1)
gene=as.data.frame(gene)
gene_list=gene[,c("fold_change","pvalue")]
gene_list[,"fold_change"]=log2(gene_list[,"fold_change"])
colnames(gene_list)=c("logFC","P.Value")
gene_list$threshold = as.factor(abs(gene_list$logFC) > 2 & gene_list$P.Value < 0.001)
gene_list$GeneSym=rownames(gene_list)

png("./result/PDAC-NET-volcano.png",width = 800, height = 800)
ggplot(data=gene_list, aes(x=logFC, y=-log10(P.Value), colour=threshold)) + geom_point(alpha=0.1, size=1.75)+
  scale_color_manual(values=c("grey", "red","green"))  + xlim(c(-8, 8)) + ylim(c(0, 12)) +xlab("log2 fold change") + ylab("-log10 p-value")
dev.off()
PDAC-NET-volcano.png

R語(yǔ)言繪制火山圖r
在ggplot2的aeS()函數(shù)中使用顏色

練習(xí)2：仔細(xì)研讀“Differentially Expressed mRNAs (Pass Volcano Plot).xls”文件，回答下述問(wèn)題

1. 該芯片是哪家芯片公司生產(chǎn)？寫(xiě)出追蹤過(guò)程

本次的數(shù)據(jù)是GSE43795，腦子不動(dòng)，直接甩給GEO的搜索引擎
GSE43795的平臺(tái)信息

那么就是Illumina咯。

2. 根據(jù)這組數(shù)據(jù)，自己計(jì)算P-value, Fold Change這兩列，并和這組數(shù)據(jù)比較，驗(yàn)證自己的結(jié)果是否正確

對(duì).xls文件的處理涉及到按sheet讀取和給定范圍讀取，readxl包是個(gè)相對(duì)更好的方法。（相比openxlsx而言，而且這個(gè)只能讀取xlsx）

library(readxl)
library(tidyverse)
path <- "./data/Differentially Expressed mRNAs (Pass Volcano Plot).xls"
up <- read_excel(path,sheet=1,range="A18:I2036",col_names = TRUE)
down <- read_excel(path,sheet=2,range="A18:I993",col_names = TRUE)
full <- merge(up, down, all=TRUE, sort=TRUE)
ref <- full[,c(4,2,3,9)]
final <- merge(information,ref,by.x = 'X', by.y = 'GeneSymbol')

R package：readxl

然后就是計(jì)算和拼表的過(guò)程 看excel可以知道是T vs N的比對(duì)結(jié)果
excel里的

才疏學(xué)淺，只能選取現(xiàn)成的歸一化數(shù)據(jù)開(kāi)干啦

gene <- full[,c(9,32:37)] #獲取歸一化的T和N數(shù)據(jù)
means <- rbind(apply(gene[,2:4],1,mean),apply(gene[,5:7],1,mean))
means <- t(means)
gene_infor <- data.frame(matrix(ncol = 5, nrow = 0))

for (i in  1:length(rownames(gene))){ # t-test pvalue
  gene_test=t.test(gene[i,2:4],gene[i,5:7]) # t test for PDAC vs NET
  mean_ratio=2^( means[i,1]- means[i,2])   #  mean fold change
  gt=c(gene[i,1],means[i,],as.vector(mean_ratio),gene_test[[3]])
  gene_infor<-rbind(gene_infor,gt)
}
cnames <- c("gene_ID","Tmean","Nmean","tfoldchange","pval")
colnames(gene_infor)<-cnames

p_value<-gene_infor[,5] 
result_test=as.matrix(gene_infor)
fdr<-p.adjust(p_value,method="BH",length(p_value))
result_fdr=cbind(result_test,fdr) #增加fdr列
test_fdr=result_fdr[result_fdr[,4]>(2)&result_fdr[,5]<(0.05),] 
#foldchange要2以上 p要0.05 FDR無(wú)要求
col_names=colnames(gene_infor)
colnames(test_fdr)=c(col_names,"fdr")
myresults <- test_fdr[,c(1,4,5,6)]

final <- merge(myresults,ref,by.x = 'gene_ID', by.y = 'GeneSymbol')

不過(guò)結(jié)果還是有點(diǎn)毛病的
結(jié)果 左邊三列是自己的 右邊三列是公司的
  1. 黃框圈出來(lái)的有很多重復(fù)出現(xiàn)的基因
  2. 只有紅框圈出來(lái)的結(jié)果是匹配的，F(xiàn)DR值可以說(shuō)是全錯(cuò)hhh 看了表格里用的都是BH方法，不至于差這么多吧！應(yīng)該還是我的問(wèn)題。先把作業(yè)交了慢慢改。。？
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容