生长在一个生物实验室,除了没人教,还有个特点就是杂活多,师兄师姐测完转录组就会找我,-师妹帮我分析下吧~ 这样也好,锻炼自己的机会多多,以生信自学路上若干前辈为目标,期待见证自我蜕变的那天。
1 读取数据,去掉缺失值
setwd("/Users/baiyunfan/desktop")
data<-read.csv("data.csv",header = T,sep = ",")
data$Symbol<-as.character(data$Symbol)
这批数据中,O1,O2,O3,O4为一组,代表老年人四个时间点的细胞,Y1, Y2, Y3, Y4为一组,代表年轻人四个时间点的细胞
2 发现Symbol中有好多缺失值,去掉。同时只保留FPKM几列和symbol
data.na<-data[!is.na(data$Symbol),c(3:11)]
3 方差分析数据整理
- 关于选方差分析还是T检验,在之前的统计学文集中已经讲过了,请需要的同学自行查阅
##由于我想筛选不同时间点具有显著差异的基因,所以1,2,3,4为一组。整理成做方差分析需要的格式
library(tidyr)
data_gather<-gather(data.na,sample,expr,-Symbol)
sample_O<-gsub("Y1_FPKM","O1_FPKM",data_gather$sample)
sample_O<-gsub("Y2_FPKM","O2_FPKM",sample_O)
sample_O<-gsub("Y3_FPKM","O3_FPKM",sample_O)
sample_O<-gsub("Y4_FPKM","O4_FPKM",sample_O)
data_gather$sample<-sample_O
4 方差分析
pvalues<-sapply(data2$Symbol,function(i){
res<-aov(expr~sample,data=subset(data2,Symbol==i))
summary(res)[[1]]$"Pr(>F)"[1]
})
data_p<-cbind(data_gather,pvalues)
##这么算会有很多重复的P值,去掉
data_p<-data_p[!duplicated(data_p$Symbol),]
##挑选显著的基因
data_sig<-data_p[which(data_p$pvalues<0.05),]
##按顺序排列
data_sig_order<-data_sig[order(data_sig$pvalues),]
5 去掉表达量低的基因
###将显著基因及其表达值整理出来
sig_gene<-data_sig_order$Symbol
sig_gene_expr<-data1[data1$Symbol==sig_gene[1],]
for(i in 2:length(sig_gene)){
m<-data1[data1$Symbol==sig_gene[i],]
sig_gene_expr<-rbind(test,m)
}
##去掉表达量<1的基因
library(edgeR)
rownames(sig_gene_expr)<-sig_gene_expr$Symbol
sig_gene_expr<-sig_gene_expr[,-9]
sig_gene_droplow<-sig_gene_expr[rowMeans(sig_gene_expr)>1,]