firebrowse下載TCGA數(shù)據(jù)
##用R處理數(shù)據(jù)
dat = read.table("mirna.txt",header=T,sep="\t")
View(head(dat))
dim(dat)
dat1 = dat[-1,c(1,seq(3,ncol(dat)-1,3))]
View(head(dat1))
dat2 = dat1[,-1]
rownames(dat2) = dat1[,1]
View(head(dat2))
library(stringr)
dat.tumor = dat2[,str_detect(colnames(dat2),'.01')]
dim(dat.tumor)
dat.normal = dat2[,str_detect(colnames(dat2),'.11A')]
dim(dat.normal)
dat.final = cbind(dat.tumor,dat.normal)
write.table(dat.final,"datExp.txt",col.names = T,row.names = T,quote=F,sep="\t")
dat.exp = read.table("datExp.txt",header=T,sep="\t")
View(head(dat.exp))
dat.exp[dat.exp == 0] = 1
dat.exp = log(dat.exp,2)
View(dat.exp)
#####-----
pheno = NULL
pheno$sample = colnames(dat.exp)
pheno$group = c(rep('tumor',ncol(dat.tumor)),rep('normal',ncol(dat.normal)))
pheno = as.data.frame(pheno)
head(pheno)
library(limma)
Group = factor(pheno$group,levels=c('tumor','normal'))
design = model.matrix(~0+Group)
colnames(design) <- c('tumor','normal')
design
fit <- lmFit(dat.exp, design)
#構(gòu)建比對(duì)模型，比較兩個(gè)條件下的表達(dá)數(shù)據(jù)
contrast.matrix <- makeContrasts(tumor-normal,#1
levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
#貝葉斯檢驗(yàn)
fit2 <- eBayes(fit2)
#找出差異基因檢驗(yàn)結(jié)果并輸出符合條件的結(jié)果
all.deg = NULL
diff = topTable(fit2,adjust.method="fdr",coef=1,p.value=0.05,
lfc=log(2,2),number=50000,sort.by = 'logFC')
write.table(diff,'DEG.txt',col.names=T,row.names=T,quote=F,sep="\t")
##火山圖
diff0 = diff
P.value = diff0$adj.P.Val
FC = diff0$logFC
df <- data.frame(P.value,FC)
df$threshold = as.factor(abs(df$FC) > log(1.5,2) & df$P.value < 0.05)
levels(df$threshold) = c('grey','red')
df$logp = -log10(df$P.value)
#using base
plot(x = df[,2],y = df[,4], pch=16, col=df$threshold,
xlab="Fold change", ylab="-Log10(Pvalue)",
cex=0.5, main="Volcano Plot")
abline(v=c(-log(1.5,2),log(1.5,2)), h=-log10(0.05),col="green")
#using ggplot2
library(ggplot2)
ggplot(data=df,aes(x=FC,y=-log10(P.value)))+
geom_point(color=df$threshold)
