首页 > 学术百科

UCSCXena数据下载后的后续分析代码——转载自麦麦大人

UCSCXena数据下载后的后续分析代码——转载⾃麦麦⼤⼈#已有UCSCXena数据，进⾏后续分析

rm(list=ls())

options(stringsAsFactors = F)

集中空调系统##可以直接读没有解压的数据，mRNA-SEQ数据

mRNA_HiSeqV2=read.table(“HiSeqV2”,header = T,sep = ‘t‘)

dim(mRNA_HiSeqV2)

mRNA_HiSeqV2[1:4,1:4]

#查看NA的数据

it(mRNA_HiSeqV2))

save(mRNA_HiSeqV2,file=“mRNA_HiSeqV2.Rdata”)

固态高频

##临床信息

mRNA_clinical=read.table(“CESC_clinicalMatrix” ,header = T,sep = ‘t‘, quote = “”)

dim(mRNA_clinical)

mRNA_clinical[1:4,1:4]

save(mRNA_clinical,file=“mRNA_clinical.Rdata”)

##⽣存相关信息

mRNA_survival=read.table(“” ,header = T,sep = ‘t‘)

dim(mRNA_survival)

mRNA_survival[1:4,1:4]

save(mRNA_survival,file=“mRNA_survival.Rdata”)

rm(list=ls())

options(stringsAsFactors = F)

load(“mRNA_HiSeqV2.Rdata”)

>>#

# 这⾥需要解析TCGA数据库的ID规律，来判断样本归类问题。mp.weixin.qq/s/Ph1O6V5RkxkyrKpVmB5ODA

#01–09是癌症，10–19是正常，20–29是癌旁

expr <- mRNA_HiSeqV2

expr[1:4,1:4]

rownames(expr)=expr[,1]

expr=expr[,-1]

expr[1:4,1:4]

#通过数字来判断样本类型

ls <- unlist(substr(colnames(expr),14,15))

table(ls)

group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,‘tumor’,‘normal’)

table(group_list)

it(expr)

dim(exprSet)

exprSet[1:4,1:4]

exprSet=exprSet[,group_list==“tumor”]

dim(exprSet)

##xena下载的数据是经过log的count值，所以要还原才可以进⾏差异分析

exprSet <- 2^exprSet-1

dim(exprSet)

exprSet[1:4,1:4]

#并且还要取整数

exprSet <- floor(exprSet)

exprSet[1:4,1:4]

save(exprSet,file=“exprSet_DEG.Rdata”)

#四、分组

拉线绝缘子

#根据感兴趣的基因的表达量进⾏分组，之后进⾏差异分析。

load(“exprSet_DEG.Rdata”)

####根据感兴趣的基因⾼低表达做差异分析，这⾥⽤了中位数，也可以⽤平均值

group_list=ifelse(as.numeric(exprSet[“MYC”,])>median(as.numeric(exprSet[“MYC”,])),‘high’,‘low’) table(group_list)

class(group_list)

group_list <- factor(group_list)

class(group_list)

#五、差异分析

#转录组数据最常⽤DESeq2包进⾏差异分析

library(DESeq2)

colData <- data.frame(row.names=colnames(exprSet),

group_list=group_list)

colData

dds <- DESeqDataSetFromMatrix(countData = exprSet,

colData = colData,

design = ~ group_list)

dds <- DESeq(dds)

res <- results(dds,

contrast=c(“group_list”,“high”,“low”))

head(res)

##sort by padj 进⼀步优化，安装标准化的P值排序

resOrdered <- res[order(res$padj),]

head(resOrdered)

DEG =as.data.frame(resOrdered)

DESeq2_DEG = na.omit(DEG)

head(DESeq2_DEG)

nrDEG=DESeq2_DEG#[,c(2,6)]

#colnames(nrDEG)=c(‘log2FoldChange’,‘pvalue’)

head(nrDEG)

#六、筛选差异基因

#对得到的差异基因，按照fold change和P值，制定⾃⼰的标准，进⾏筛选。

##筛选

#get diff_gene 得到差异基因，选P值及logFC⼤于1的

chunqingintdiff_gene <- subset(nrDEG,pvalue<0.01 & (log2FoldChange>=1 | log2FoldChange<=-1)) head(diff_gene)

女夭#给数据标记上调下调

resdata <- diff_gene

resdata$significant <- c(rep(“unchange”,nrow(resdata)))

resdata$significant[resdata$pvalue<=0.01&resdata$log2FoldChange>=1] <- “up”resdata$significant[resdata$pvalue<=0.01&resdata$log2FoldChange<=-1] <- “down”

#查看最新数据的前五⾏

head(resdata)

resdata[1:40,]

diff_gene <- resdata

table(diff_gene$significant)

#输出⽂件为TXT格式

write.table(x=diff_gene,

file=“diff_result_”,

quote=F,

sep = “t“,

row.names = T,

col.names = T)

save(nrDEG,DESeq2_DEG,diff_gene,file=“DESeq2_DEG_result.Rdata”)

load(“DESeq2_DEG_result.Rdata”)

#七、热图

rm(list=ls())

options(stringsAsFactors = F)

load(“DESeq2_DEG_result.Rdata”)

load(“exprSet_DEG.Rdata”)

####根据感兴趣的基因⾼低表达做差异分析，这⾥⽤了中位数，也可以⽤平均值

group_list=ifelse(as.numeric(exprSet[“MYC”,])>median(as.numeric(exprSet[“MYC”,])),‘high’,‘low’) table(group_list)

class(group_list)

group_list <- factor(group_list)

class(group_list)

>#3提取差异最⼤的50个基因画热图,⽤P值判断

## heatmap

library(pheatmap)

choose_gene=head(rownames(diff_gene),50) ## 50 maybe better

choose_matrix=exprSet[choose_gene,]

#把太⼩的值去掉

boxplot(choose_matrix)

n=t(scale(t(choose_matrix)))

boxplot(n)

n[n>2]=2

n[n< –2]=-2

n[1:4,1:4]

boxplot(n)

#分组信息

annotation_col=data.frame(group=group_list)

row.names(annotation_col) <- colnames(exprSet)

annotation_col

#pheatmap(choose_matrix,filename = paste0(n,‘_need_DEG_top50_heatmap.png’)) pheatmap(n,annotation_col=annotation_col,show_rownames=F,show_colnames=F)

>#3提取差异最⼤的50个基因画热图,⽤logFC判断

## heatmap

library(pheatmap)

diff_gene<- diff_gene[order(abs(diff_gene$log2FoldChange),decreasing = T),]

head(diff_gene)

exprSet <- exprSet[,order(group_list)]

exprSet[1:4,1:4]

choose_gene=head(rownames(diff_gene),50) ## 50 maybe better

choose_matrix=exprSet[choose_gene,]

#把太⼩的值去掉

boxplot(choose_matrix)

n=t(scale(t(choose_matrix)))

boxplot(n)

n[n>2]=2

n[n< –2]=-2

n[1:4,1:4]

boxplot(n)

#分组信息

annotation_col=data.frame(group=group_list)

郑道传row.names(annotation_col) <- colnames(exprSet)

annotation_col

#pheatmap(choose_matrix,filename = paste0(n,‘_need_DEG_top50_heatmap.png’)) pheatmap(n,annotation_col=annotation_col,

show_rownames=F,show_colnames=F)

#⼋、⽕⼭图

#做⽕⼭图数据处理

>#⽕⼭图

rm(list = ls()) ## 魔幻操作，⼀键清空~

options(stringsAsFactors = F)

#导⼊数据,⽕⼭图需要的数据前期处理

load(“DESeq2_DEG_result.Rdata”)

head(DESeq2_DEG)

library(ggpubr)

df<-DESeq2_DEG

df<-DESeq2_DEG[-1,] ##MYC的P值太⼩，为了作图好看

attach(df)

##⽕⼭图横坐标是logFC，纵坐标是-log10(P.Value)

plot(log2FoldChange,-log10(pvalue))

df$v=-log10(pvalue) #df新增加⼀列‘v’，值为-log10(P.Value)

本文发布于:2024-09-23 11:12:08，感谢您对本站的认可！

本文链接：https://www.17tex.com/xueshu/570802.html

上一篇：利用UCSCXena做TCGA数据库的生存曲线分析

下一篇：UCSB、UCI、UCSD、UCD 的区别

标签：数据分析差异判断表达样本标准化基因

留言与评论（共有 0 条评论）