讀取表達矩陣

setwd("E:/8.差異表達基因/")
a <- read.table("RNAmatrix.txt",header = T)

去重復

b <- a[a$gene_name,]

去掉第一行第一列

rownames(c) <- c[,1]

c <- c[,-1]

選擇所需數(shù)據(jù)

a <- c[,-()]

c <- c[,10:15]

dat1<-as.data.frame(c)
dim(dat1)
dat1[1:4,1:4]

很多表達量為0的樣本，直接選擇在某個基因如果在3個樣品中的表達量為零，則直接舍去。

apply(dat1,1,function(x){sum(floor(x)==0)>3})
dat1<-dat1[!apply(dat1,1,function(x){sum(floor(x)==0)>3}),]
dim(dat1)
head(dat1)
write.csv(dat1,"dat1.csv")

boxplot(dat1)

差異很大取log歸一化

dat3 <- log(dat1)#下游分析的結(jié)果有缺失值，故選擇log(dat2 + 1)

dat4 <- log(dat1 + 1)

boxplot(dat3)

boxplot(dat4)

write.csv(dat3,"dat3.csv")

write.csv(dat4,"dat4.csv")

差異基因分析

library(limma)

group <- c(rep("normal",50),rep("cancer",374)) 
head(group)
View(group)

group <- factor(group)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
design


contrast.matrix <- makeContrasts(normal - cancer,  
                                 levels=design)

contrast.matrix

fit <- lmFit(dat4,design)
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2)

allDiff1=topTable(fit2,adjust='fdr',coef=1,number=Inf) 
save(dat4,group,allDiff1,file = "RNAmatrix_result.txt")

write.csv(allDiff1, file = "normal_vs_cancer.csv" )