前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >适合用于FPKM数据求差异基因的ballgown算法

适合用于FPKM数据求差异基因的ballgown算法

作者头像
用户1359560
发布2021-12-06 18:12:14
1.2K0
发布2021-12-06 18:12:14
举报
文章被收录于专栏:生信小驿站

对于FPKM表达数据时,Edger,limma,和deseq等算法并不合适。而ballgown是针对于FPKM数据开发的差异基因算法,可以尝试。 示例数据如下所示:

代码语言:javascript
复制
# --------------------------------------------------------
# 
# 
# 
# --------------------------------------------------------

setwd("D:\\SCIwork\\F5\\DEG")


library(tidyr)
library('ballgown')
load("mRNA_exprSet.Rda")

mRNA_exprSet[1:4,1:4]

# --------------------------------------------------------
# 
# 
# 
# --------------------------------------------------------

mRNA_exprSet <- mRNA_exprSet %>%
  tidyr::separate(gene_id, c("gene_name",
                             "gene_id",
                             "gene_biotype"), 
                  sep = " \\| ")

mRNA_exprSet <- mRNA_exprSet[,-(2:3)]
index <- duplicated(mRNA_exprSet$gene_name)
mRNA.data <- mRNA_exprSet[!index,]
dim(mRNA.data)


# --------------------------------------------------------
# 
# 
# 
# --------------------------------------------------------

#包含基因名的第一列转为行名
BLCA_fpkm_data = mRNA.data
rownames(BLCA_fpkm_data) = BLCA_fpkm_data[,1]
BLCA_fpkm_data  =  BLCA_fpkm_data[c(-1)]

#生成分组文件
load("mRNA_exprSet.Rda")
metadata <- data.frame(names(mRNA_exprSet)[-1])
for (i in 1:length(metadata[,1])) {
  num <- as.numeric(substring(metadata[i,1],14,15))
  if (num %in% seq(1,9)) {metadata[i,2] <- "T"}
  if (num %in% seq(10,29)) {metadata[i,2] <- "N"}
}

names(metadata) <- c("TCGA_id","group")

metadata$group <- as.factor(metadata$group)


# --------------------------------------------------------
# 
# 
# 
# --------------------------------------------------------


result_diff = stattest(gowntable = BLCA_fpkm_data ,
                       pData = metadata , 
                       covariate = "group" , 
                       getFC = TRUE , 
                       log =TRUE,
                       meas='FPKM',
                       feature="gene")


result_diff$LogFC <- log2(result_diff$fc)

result_diff$LogFC_abs <- abs(result_diff$LogFC)

write.csv (result_diff, "mRNA_BLCA_fpkm_diff.csv", row.names = F)

# --------------------------------------------------------
# 
# 
# 
# --------------------------------------------------------

foldChange =0.5

padj=0.05

diffSig=result_diff[which(result_diff$pval< padj & result_diff$LogFC_abs > foldChange),]

dim(diffSig)

write.csv(diffSig, file="diffSig_mRNA_BLCA.csv")
本文参与 腾讯云自媒体同步曝光计划,分享自作者个人站点/博客。
原始发表:2021/7/30 上,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 作者个人站点/博客 前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档