rm(list = ls())
setwd("/data3/ff/project/XCL/pro")
##########
library(data.table)
all<-fread("pro.txt",data.table = F)
colnames(all)
###########1,8,15,22
HFDVAT.CDVAT<-all[,c(1:7)]
HFD_S.CD_S<-all[,c(8:14)]
HFDKPC.CDKPC<-all[,c(15:21)]
S.HFDKPC.S.CDKPC<-all[,c(22:28)]
################################测试
library(impute)
library(limma)
library(dplyr)
##################################################################对于非向量的循环，推荐构建list，然后用lapply
#####
##################################过滤一些全为0的行
HFD_S.CD_S
class(HFD_S.CD_S[1,2])
#过滤至少有一个样本不为0
keep<-rowSums(HFD_S.CD_S> 0) > 1
HFD_S.CD_S <- HFD_S.CD_S[keep,]
data<-list(HFDVAT.CDVAT=HFDVAT.CDVAT,
HFD_S.CD_S=HFD_S.CD_S,
HFDKPC.CDKPC=HFDKPC.CDKPC,
S.HFDKPC.S.CDKPC=S.HFDKPC.S.CDKPC)
#data<-HFDVAT.CDVAT
#data<- HFD_S.CD_S
data<-HFDKPC.CDKPC
data<- S.HFDKPC.S.CDKPC
##################################################################函数体knn.diff
knn.diff <- function(data){
print(data)
data <- na.omit(data)
#将所有0值替换为NA
data[data == 0] <- NA
####去重基因名
table(duplicated(data$`Gene name`))
#对重复基因名取平均表达量,获得表达矩阵
if(sum(duplicated(data$`Gene name`))>0) #判断是否有重复基因
data<-avereps(data,ID=data$`Gene name`) %>% as.data.frame
##更改行名
row.names(data)<-data[,1]
gene<-data[,1]#取出基因名
data<-data[,-1]
#########补全
data=as.matrix(data)
#数据补缺
mat=impute.knn(data)
data=mat$data %>% as.data.frame
#转为数值型
data=as.data.frame(lapply(data,as.numeric))
row.names(data)<-gene
#对数据进行矫正
data=normalizeBetweenArrays(as.matrix(data))
#class(data[1,1])
######################################################差异
pFilter=0.05 #p值临界值
logFCfilter=0 #logFC临界值
conNum=3 #sample1组样品数目
treatNum=3 #sample2组样品数目
#读取输入文件
outTab=data.frame()
group=c(rep(1,conNum),rep(2,treatNum))
data=as.matrix(data)
###################################################去掉数据是恆量
#编写函数，返回一行有几个不同值
uniq.num <- function(x){
y<-length(unique(x[1:3]))+length(unique(x[4:6]))
return(y)
}
#length () 函数用于获取或设置向量 (列表)或其他对象的长度
#按行运行函数，计算每行的不同值
uniq.num.df <- data.frame(num = apply(data, 1, uniq.num))
#过滤掉uniq.num为1或2的,留下＞2的
uniq.num.df.2 <- filter(uniq.num.df, num>2)
#过滤原数据
data <- data[rownames(uniq.num.df.2), ]
#差异分析
for(i in row.names(data)){
geneName=i
rt=rbind(expression=data[i,],group=group)
rt=as.matrix(t(rt))
tTest<-t.test(expression ~ group, data=rt)
pvalue=tTest$p.value
conGeneMeans=mean(data[i,1:conNum])
treatGeneMeans=mean(data[i,(conNum+1):ncol(data)])
logFC=log2(treatGeneMeans/conGeneMeans)
conMed=median(data[i,1:conNum])
treatMed=median(data[i,(conNum+1):ncol(data)])
diffMed=treatMed-conMed
if( ((logFC>0) & (diffMed>0)) | ((logFC<0) & (diffMed<0)) ){
outTab=rbind(outTab,cbind(gene=i,conMean=conGeneMeans,treatMean=treatGeneMeans,logFC=logFC,pValue=pvalue))
}
}
#table(duplicated(outTab$gene))
#对重复基因名取平均表达量,一模一样
#if(sum(duplicated(outTab$gene))>0) #判断是否有重复基因
#outTab<-avereps(outTab,ID=outTab$gene) %>% as.data.frame
#改行名
row.names(outTab)<-outTab$gene
outTab<-outTab[,-1]
#将表达值和结果合并
results<-merge(data,outTab,by="row.names",all.x=TRUE)
return(results)
}
HFDVAT.CDVAT.results<-results
HFD_S.CD_S.results<-results
HFDKPC.CDKPC.results<-results
S.HFDKPC.S.CDKPC.results<-results
write.csv(HFDVAT.CDVAT.results,"HFDVAT.CDVAT.results.csv")
write.csv(HFD_S.CD_S.results,"HFD_S.CD_S.results.csv")
write.csv(HFDKPC.CDKPC.results,"HFDKPC.CDKPC.results.csv")
write.csv(S.HFDKPC.S.CDKPC.results,"S.HFDKPC.S.CDKPC.results.csv")
###########写函数过滤
###########放入自变量
list<- list(HFDVAT.CDVAT.results,HFD_S.CD_S.results,HFDKPC.CDKPC.results,S.HFDKPC.S.CDKPC.results)
#############写过滤函数
filler <- function(x){
x<-na.omit(x)
x<-x[(as.numeric(as.vector(x$logFC))>0.1 & as.numeric(as.vector(x$pValue))<0.05),]
y<-x$Row.names
return(y)
}
#####应用函数
results.list<- lapply(list,filler)
####取交集
a<-intersect(results.list[[1]],intersect(results.list[[2]],results.list[[3]]))
b<-intersect(results.list[[1]],intersect(results.list[[2]],results.list[[4]]))
网友评论