analysis-r's People
analysis-r's Issues
fig kegg
rm(list = ls())
options(stringsAsFactors = F)
#step1 get genes_expr
library(AnnoProbe)
#suppressPackageStartupMessages(library(GEOquery))
#gset=AnnoProbe::geoChina('GSE6364')
a=read.table('rifpatient_tissue.txt',sep='\t',quote = "",fill = T,
comment.char = "!",header = T) # 提取表达矩阵
rownames(a)=a[,8]
a <- a[,-1]
expMatrix <- a[,1:6]
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpms <- apply(expMatrix,2,fpkmToTpm)
tpms[1:3,]
colSums(tpms)
genes_expr<-tpms
group_list=c(rep('CTL',3),rep('RIF',3))
save(genes_expr,group_list,file = 'step1-genes_expr.Rdata')
load('step1-genes_expr.Rdata')
temps=as.data.frame(genes_expr)
temps$id=row.names(temps)
b=a[,7:10]
temps=merge.data.frame(temps,b,by.x="id",by.y = 'symbol')
head(temps)
rownames(temps)=temps[,1]
temps <- temps[,-1]
degexpr=temps
save(degexpr,group_list,file = 'deg-genes_expr.Rdata')
deg=degexpr[,7:9]
head(deg)
save(deg,file = 'deg.Rdata')
rm(list = ls())
options(stringsAsFactors = F)
load('deg.Rdata')
load('step1-genes_expr.Rdata')
head(deg)
visualization
need_deg=data.frame(symbols=rownames(deg), logFC=deg$logFC, p=deg$Pvalue)
deg_volcano(need_deg,1)
deg_volcano(need_deg,2)
deg_heatmap(deg,genes_expr,group_list)
deg_heatmap(deg,genes_expr,group_list,30)
check_diff_genes('GPX4',genes_expr,group_list)
check_diff_genes('FOXO1',genes_expr,group_list)
check_diff_genes('ACSL4',genes_expr,group_list)
check_diff_genes('MAP3K4',genes_expr,group_list)
check_diff_genes(c('GPX4','ACSL4','LPCAT3'),genes_expr,group_list)
#tonglujiyinbianhua
if(F){
#ferroptosis
library(KEGGREST)
exprSet<-genes_expr
#cg <- KEGGREST::keggGet("hsa04216")[[1]]$GENE
cg <- KEGGREST::keggGet("hsa04350")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,exprSet,group_list)
p1<-check_diff_genes( cg ,exprSet,group_list)
p1
ggsave("RIFCTL TGFb.png",p1,height = 10, width = 5)
cg=cg[1:88]
ferroptosisrelategene=cg
ferrodeg=a[ferroptosisrelategene,]
save(ferrodeg,file = 'ferrodeg.Rddta')
check<-genes_expr[cg,]
check<-na.omit(genes_expr[cg,])
ac=data.frame(g=group_list)
rownames(ac)=colnames(check)
pheatmap(check,show_colnames =F,show_rownames = T) #对dat按照cg取行,所得到的矩阵来画热图
n=t(scale(t(check)))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于scale这个函数应用在不同组数据间存在差异时,需要行名为样本,因此需要用t(dat[cg,])来转换,最后再转换回来
n=t(scale(t(check)))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于scale这个函数应用在不同组数据间存在差异时,需要行名为样本,因此需要用t(dat[cg,])来转换,最后再转换回来
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = T)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #将ac的行名也就分组信息(是‘no TNBC’还是‘TNBC’)给到n的列名,即热图中位于上方的分组信息
pheatmap(n,show_colnames =F,
show_rownames = T,
cluster_cols = F,
annotation_col=ac,
filename = 'Pre ESE ferroptosis deg.png',width = 8, height= 9) #列名注释信息为ac即分组信息
#autophagy
library(KEGGREST)
cg <- KEGGREST::keggGet("hsa04140")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,gene_exCPEMS,group_list)
#apotosis
library(KEGGREST)
cg <- KEGGREST::keggGet("hsa04210")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,gene_exCPEMS,group_list)
细胞衰老
map04218
library(KEGGREST)
cg <- KEGGREST::keggGet("map04933")[[1]]$GENE
cg=as.character(sapply(cg[seq(2,length(cg),by=2)], function(x) strsplit(x,';')[[1]][1]))
check_diff_genes( cg ,genes_expr,group_list)
}
for volcano
if(F){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(Pvalue))
library(ggpubr)
df=nrDEG
df$v= -log10(Pvalue) #df新增加一列'v',值为-log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$Pvalue>0.05,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
ifelse( df$logFC >1,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
ifelse( df$logFC < -1,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1.5 的基因,再if 判断:如果logFC <1.5,则为down(下调)基因,否则为stable基因
)
table(df$g)
df$name=rownames(df)
head(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = c('ATF3'), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano2.png')
}
for heatmap
if(F){
group_list
table(group_list)
x=deg$logFC #deg取logFC这列并将其重新赋值给x
names(x)=rownames(deg) #deg取probe_id这列,并将其作为名字给x
cg=c(names(head(sort(x),100)),#对x进行从小到大排列,取前100及后100,并取其对应的探针名,作为向量赋值给cg
names(tail(sort(x),100)))
library(pheatmap)
dat=genes_expr
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行,所得到的矩阵来画热图
n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于scale这个函数应用在不同组数据间存在差异时,需要行名为样本,因此需要用t(dat[cg,])来转换,最后再转换回来
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #将ac的行名也就分组信息(是‘no TNBC’还是‘TNBC’)给到n的列名,即热图中位于上方的分组信息
pheatmap(n,show_colnames =T,
show_rownames = F,
cluster_cols = T,
annotation_col=ac,filename = 'heatmap_top200_DEG2.png') #列名注释信息为ac即分组信息
}
write.csv(DEG,file = 'DEG.csv')
rm(list = ls())
load('RIFanno_DEG.Rdata')
a=DEG
row.names(a)=DEG$ENTREZID
b=as.numeric(a$logFC)
names(b)=row.names(a)
head(b)
library('pathview')
p <- pathview(gene.data = b,
pathway.id = "04216",
species = "hsa",
out.suffix = "RIFCTL",
#limit = list(gene=as.integer(max(abs(a))))
)
p <- pathview(gene.data = b,
pathway.id = "04110",
species = "hsa",
out.suffix = "RIFCTL2",
limit = 2
#limit = list(gene=as.integer(max(abs(a))))
)
p <- pathview(gene.data = b,
pathway.id = "04350",
species = "hsa",
out.suffix = "RIFCTL1",
limit = 1
#limit = list(gene=as.integer(max(abs(a))))
)
check_diff_genes <- function(gene,genes_expr,group_list ){
if(length(gene)==1){
# gene='NKILA'
if(! gene %in% rownames(genes_expr)){
stop(paste0(gene,' in not in your expression matrix'))
}
df=data.frame(value=as.numeric(genes_expr[gene,]),
group=group_list)
library(ggpubr)
ggboxplot(df, "group", "value",
color = "group", palette =c("#00AFBB", "#E7B800"),
add = "jitter", shape = "group")
}else{
library(pheatmap)
cg=gene
cg=cg[cg %in% rownames(genes_expr) ]
warning(paste0('Only ',length(cg),' in ',length(gene),' genes are in your expression matrix'))
if(length(cg)<1){
stop('None of the gene in your expression matrix')
}
n=t(scale(t(genes_expr[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
ac=data.frame(group_list=group_list)
rownames(ac)=colnames(n)
pheatmap(n,show_colnames =F,show_rownames = T,cluster_cols= F,
annotation_col=ac)
}
}
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.