练习1：从PDAC-NET.txt表达谱中获得PDAC和NET的差异基因，(显著性水平为0.05 ，fdr 控制在0.01)并对上述差基因画热图。
heatmap必备包
library(gplots)
以下内容都是模仿老师的代码，不过我的文件夹结构和老师不完全相同，小改
读入表达谱
gene_exp <- read.table("./data/PDAC-NET.txt", sep="\t", header=T,row.names=1)
gene <- as.matrix(gene_exp)
对输出文件的创造（相比老师，这里把append参数改为FALSE，为了后续电脑重启继续做，跑一遍的时候不会做出重复内容的文件）
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)
获取每个基因在样本中的均值
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后，得到的差异基因
画热图
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
选做：对火山图（volcano plot）颜色重置：高表达红色表示，低表达绿色表示，无差异表达灰色表示
2022.10.19 ↓成功达成颜色要求 顺便学习了一下ggplot2的使用 画了个辅助线
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 ↓ 画的还不是特别满足要求 可能之后会改吧
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
练习2：仔细研读“Differentially Expressed mRNAs (Pass Volcano Plot).xls”文件，回答下述问题
1. 该芯片是哪家芯片公司生产？写出追踪过程
本次的数据是GSE43795，脑子不动，直接甩给GEO的搜索引擎
GSE43795的平台信息
那么就是Illumina咯。
2. 根据这组数据，自己计算P-value, Fold Change这两列，并和这组数据比较，验证自己的结果是否正确
对.xls文件的处理涉及到按sheet读取和给定范围读取，readxl包是个相对更好的方法。（相比openxlsx而言，而且这个只能读取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')
然后就是计算和拼表的过程 看excel可以知道是T vs N的比对结果
excel里的
才疏学浅，只能选取现成的归一化数据开干啦
gene <- full[,c(9,32:37)] #获取归一化的T和N数据
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无要求
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')
结果 左边三列是自己的 右边三列是公司的
- 黄框圈出来的有很多重复出现的基因
- 只有红框圈出来的结果是匹配的，FDR值可以说是全错hhh 看了表格里用的都是BH方法，不至于差这么多吧！应该还是我的问题。先把作业交了慢慢改。。？
