贝叶斯t-test
setwd("E:/TCGA-mRNA/new/merged_tpm/")
data<-read.table("COAD.sort.rpkm.txt",head=T,sep="\t",row.names=1)
source("d:/desktop/bayesreg.R")
library(multtest)
T<-bayesT(data,numC= ,numE= ,ppde=TRUE,betaFit=1,bayes=TRUE,doMulttest=1)
write.table(T,"bayes",quote=F,col.names=TRUE,sep="\t"))
limma
library(edgeR)
library(limma)
lncRNA<-read.table("D:/Desktop/BRCA/RNAseq/filter_lncRNA.txt",header = TRUE,sep = "\t",row.names = 1)
group<-factor(c(rep("C",113),rep("E",1109)))
genelist <- DGEList(counts = lncRNA,group = group)
keep <- rowSums(cpm(genelist)>0.5)>=2
genelist.filted <- genelist[keep,,keep.lib.sizes=FALSE]
x <- calcNormFactors(genelist.filted,method = "TMM")
design<-model.matrix(~0+group)
colnames(design)<-c("C","E")
logCPM <- cpm(x,log = TRUE,prior.count = 3)
fit<- lmFit(logCPM,design)
cont.matrix=makeContrasts("C-E",levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit<-eBayes(fit2,trend = TRUE)
res<-topTable(fit,number = nrow(lncRNA))
write.table(res,"D:/Desktop/BRCA/RNAseq/DEG_lncRNA.txt",quote = FALSE,sep = "\t")
网友评论