专栏首页生信小驿站多分组差异分析解决方案(2)分批次差异基因后取交集

多分组差异分析解决方案(2)分批次差异基因后取交集

主要方法:如果不同分组代表着一定的趋势,例如group1,group2,group3的样本严重程度越来越重。那么就可以求group1和group2的差异基因,group2和group3的差异基因,group1和group3的差异基因,最后把三次得到的上调差异基因和下调差异基因求交集。

第一步读取数据,合并表达矩阵和分组文件

#===========================================================================


#===========================================================================


rm(list = ls(all.names = TRUE))

options(stringsAsFactors = F)

library(Matrix)



setwd('D:\\SCIwork\\F38KRT\\s2')


data <- read.csv('cdata.csv', header = T, row.names = 1)

data <- as.data.frame(t(data))

data[1:4,1:4]


normalized<-function(y) {
  
  x<-y[!is.na(y)]
  
  x<-(x - min(x)) / (max(x) - min(x))
  
  y[!is.na(y)]<-x
  
  return(y)
}

db <-  as.data.frame(apply(data,2,normalized))


data <- db

data$sample <- rownames(data)

data$sample <- chartr(old='.', new='-', x=data$sample)


setwd('D:\\SCIwork\\F38KRT\\s3')

group <- read.csv('group2.csv', header = T)

names(group)[1] <- 'sample'

group$sample <- chartr(old='.', new='-', x=group$sample)

group <- subset(group, select=c("sample", "group"))

group$subtype <- group$group

group$group  <- NULL




dt <- merge(group, data, by='sample')

dt[1:4,1:4]

dt$sample <- NULL

table(dt$subtype)

dt_total <- dt





normalized<-function(y) {
  x<-y[!is.na(y)]
  x<-(x - min(x)) / (max(x) - min(x))
  y[!is.na(y)]<-x
  return(y)}

求group1和group2的差异基因

#===========================================================================


#===========================================================================


table(dt$subtype)

dt <- dt[order(dt$subtype), ]

dt[1:4,1:4]



dt_Con <- subset(dt, dt$subtype == 'Subtype1')

dt_Con[1:4,1:4]

dt_Exp <- subset(dt, dt$subtype == 'Subtype2')

dt_Exp[1:4,1:4]




dt_Con$subtype <- paste0(dt_Con$subtype, rownames(dt_Con))

rownames(dt_Con) <- dt_Con$subtype

dt_Con$subtype <- NULL

dt_Con <- as.data.frame(t(dt_Con))




dt_Exp$subtype <- paste0(dt_Exp$subtype, rownames(dt_Exp))

rownames(dt_Exp) <- dt_Exp$subtype

dt_Exp$subtype <- NULL

dt_Exp <- as.data.frame(t(dt_Exp))


Pvalue<-c(rep(0,nrow(dt_Con)))

log2_FC<-c(rep(0,nrow(dt_Con)))

for(i in 1:nrow(dt_Con)){
  
  y=t.test(as.numeric(dt_Con[i,]),as.numeric(dt_Exp[i,]))
  Pvalue[i] <- y$p.value
  log2_FC[i] <-log2(mean(as.numeric(dt_Exp[i,]))/(mean(as.numeric(dt_Con[i,]))))
  
}


library(dplyr)

library(tidyr)

library(tibble)

# 对p value进行FDR校正
fdr=p.adjust(Pvalue, "BH") 
# 在原文件后面加入log2FC,p value和FDR,共3列;
out<- as.data.frame(cbind(log2_FC,Pvalue,fdr))
out$gene <- rownames(dt_Con)
# out <- out %>%
#   dplyr::filter(log2_FC > 0.5 & Pvalue < 0.05)


setwd('D:\\SCIwork\\F38KRT\\s5')

write.csv(out, file = 'out_S1.csv')

求group2和group3的差异基因

#===========================================================================


#===========================================================================


table(dt$subtype)

dt <- dt[order(dt$subtype), ]

dt[1:4,1:4]



dt_Con <- subset(dt, dt$subtype == 'Subtype2')

dt_Con[1:4,1:4]

dt_Exp <- subset(dt, dt$subtype == 'Subtype3')

dt_Exp[1:4,1:4]




dt_Con$subtype <- paste0(dt_Con$subtype, rownames(dt_Con))

rownames(dt_Con) <- dt_Con$subtype

dt_Con$subtype <- NULL

dt_Con <- as.data.frame(t(dt_Con))




dt_Exp$subtype <- paste0(dt_Exp$subtype, rownames(dt_Exp))

rownames(dt_Exp) <- dt_Exp$subtype

dt_Exp$subtype <- NULL

dt_Exp <- as.data.frame(t(dt_Exp))


Pvalue<-c(rep(0,nrow(dt_Con)))

log2_FC<-c(rep(0,nrow(dt_Con)))

for(i in 1:nrow(dt_Con)){
  
  y=t.test(as.numeric(dt_Con[i,]),as.numeric(dt_Exp[i,]))
  Pvalue[i] <- y$p.value
  log2_FC[i] <-log2(mean(as.numeric(dt_Exp[i,]))/(mean(as.numeric(dt_Con[i,]))))
  
}


library(dplyr)

library(tidyr)

library(tibble)

# 对p value进行FDR校正
fdr=p.adjust(Pvalue, "BH") 
# 在原文件后面加入log2FC,p value和FDR,共3列;
out<- as.data.frame(cbind(log2_FC,Pvalue,fdr))
out$gene <- rownames(dt_Con)
# out <- out %>%
#   dplyr::filter(log2_FC > 0.5 & Pvalue < 0.05)


setwd('D:\\SCIwork\\F38KRT\\s5')

write.csv(out, file = 'out_S2.csv')

求group1和group3的差异基因

#===========================================================================


#===========================================================================


table(dt$subtype)

dt <- dt[order(dt$subtype), ]

dt[1:4,1:4]



dt_Con <- subset(dt, dt$subtype == 'Subtype1')

dt_Con[1:4,1:4]

dt_Exp <- subset(dt, dt$subtype == 'Subtype3')

dt_Exp[1:4,1:4]




dt_Con$subtype <- paste0(dt_Con$subtype, rownames(dt_Con))

rownames(dt_Con) <- dt_Con$subtype

dt_Con$subtype <- NULL

dt_Con <- as.data.frame(t(dt_Con))




dt_Exp$subtype <- paste0(dt_Exp$subtype, rownames(dt_Exp))

rownames(dt_Exp) <- dt_Exp$subtype

dt_Exp$subtype <- NULL

dt_Exp <- as.data.frame(t(dt_Exp))


Pvalue<-c(rep(0,nrow(dt_Con)))

log2_FC<-c(rep(0,nrow(dt_Con)))

for(i in 1:nrow(dt_Con)){
  
  y=t.test(as.numeric(dt_Con[i,]),as.numeric(dt_Exp[i,]))
  Pvalue[i] <- y$p.value
  log2_FC[i] <-log2(mean(as.numeric(dt_Exp[i,]))/(mean(as.numeric(dt_Con[i,]))))
  
}


library(dplyr)

library(tidyr)

library(tibble)

# 对p value进行FDR校正
fdr=p.adjust(Pvalue, "BH") 
# 在原文件后面加入log2FC,p value和FDR,共3列;
out<- as.data.frame(cbind(log2_FC,Pvalue,fdr))
out$gene <- rownames(dt_Con)
# out <- out %>%
#   dplyr::filter(log2_FC > 0.5 & Pvalue < 0.05)


setwd('D:\\SCIwork\\F38KRT\\s5')

write.csv(out, file = 'out_S3.csv')

取交集

#===========================================================================


#===========================================================================


diff1 <- read.csv('out_S1.csv', header = T, row.names = 1)

diff2 <- read.csv('out_S2.csv', header = T, row.names = 1)

diff3 <- read.csv('out_S3.csv', header = T, row.names = 1)



diff1on <- subset(diff1, diff1$log2_FC > 0.2)
diff2on <- subset(diff2, diff2$log2_FC > 0.2)
diff3on <- subset(diff3, diff3$log2_FC > 0.2)

up_gene <- intersect(diff1on$gene, diff2on$gene)
up_gene <- intersect(up_gene, diff3on$gene)




diff1down <- subset(diff1, diff1$log2_FC < -0.2)
diff2down <- subset(diff2, diff2$log2_FC < -0.2)
diff3down <- subset(diff3, diff3$log2_FC < -0.2)

down_gene <- intersect(diff1down$gene, diff2down$gene)
down_gene <- intersect(down_gene, diff3down$gene)

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 多分组差异分析解决方案(1)循环T检验

    用户1359560
  • 单基因生信分析流程(2)一文解决差异分析、基因相关分析问题

    (1)读取基因表达矩阵 (2)根据基因表达量设置样本分组 (3)设置差异倍数、生成差异分析结果 (4)绘制火山图和热图

    用户1359560
  • 转录组分析的正确姿势

    生信宝典
  • 综述:高维单细胞RNA测序数据分析工具(上)

    当你的才华还撑不起你的野心时,请潜下心来,脚踏实地,跟着我们慢慢进步。不知不觉在单细胞转录组领域做知识分析也快两年了,通过文献速递这个栏目很幸运聚集了一些小伙伴...

    生信技能树jimmy
  • 【Scikit-Learn 中文文档】分解成分中的信号(矩阵分解问题) - 无监督学习 - 用户指南 | ApacheCN

    2.5. 分解成分中的信号(矩阵分解问题) 2.5.1. 主成分分析(PCA) 2.5.1.1. 准确的PCA和概率解释(Exact PCA and p...

    片刻
  • 单细胞转录组测序中的批次效应知多少? (上)

    当然,这样做的第一件问题就是数据来自不同样本,实验室和实验的数据并不能“很好地混合”。数据不能很好的联合分析主要是基于细胞注释和使用UMAP 和/或 tSNE ...

    生信技能树
  • 单细胞RNA-seq数据分析最佳实践(中)

    Luecken MD, Theis FJ. Current best practices in single-cell RNA-seq analysis: a ...

    生信技能树jimmy
  • 单细胞转录组测序数据质控

    总的来说,单细胞转录组体现异质性(个体),Bulk转录组体现平均程度(总体)。因此Bulk转录组不能区分一个样本中的不同细胞类型。

    生信交流平台
  • GEO数据库可能遇到的问题 (一)

    昨天介绍完GEO2R之后其实该和大家说一下富集分析相关的东西了(昨日链接:GEO2R差异表达分析软件)。但是,由于GEO数据库里面的数据种类比较多...

    医学数据库百科

扫码关注云+社区

领取腾讯云代金券