发布于 由Wolf_rabbit

转录组差异表达分析（UP、DOWN 和NOT_Change）

利用Rstudio新建一个名为Differential_Expression_analysis的project，利用以下代码进行分析（此步骤将获得每个基因在不同比较组合中的差异表达情况）

运行过程如有报错，请补齐   ｝

#清空当前环境中的所有变量

rm(list = ls())

{

 

  #如果没有安装"BiocManager"，则进行安装，否则跳过

{if (!requireNamespace("BiocManager", quietly = TRUE))

  install.packages("BiocManager")

}

#安装并加载R包"DESeq2"

{if (!requireNamespace("DESeq2", quietly = TRUE))

  BiocManager::install("DESeq2")

}

 

library(DESeq2)

#输入TPM数据

tpm_input <- read.table(file = "G:/R-Demo/counts_to_TPM/all_TPM.txt",

                  header = T,row.names = 1,sep = "\t")

#去除各转录组中tpm值均为0的基因

tpm <- tpm_input[rowSums(tpm_input) != 0,]

# 将TPM数据乘以1000并四舍五入，不保留小数

tpm <- round(tpm)

head(tpm)

 

# 创建条件信息:需要因子factor格式

#主要是转录组中包含重复，需要指定哪些重复是对应一个样本的，后续需要对不同样本进行比较

conditions0 <- c("HF1", "HF2", "HF3", "HGP", "HGX", "HL", "HM1", "HM2",

                        "HM3", "HP1", "HP2", "HP3", "HR", "HS", "HZZ")

 

# 嵌套循环比较两两条件

for (i in 1:(length(conditions0)-1)) {

  for (j in (i+1):length(conditions0)) {

    condition1 <- conditions0[i]

    condition2 <- conditions0[j]

 

    # 打印两个条件的比较结果

    print(paste("Comparing", condition1, "and", condition2))

    

    # 在这里执行比较的操作，根据具体需求进行相应的处理

#将需要比较的对象与数据中的各列一一对应

conditions <- as.factor(rep(c(condition1, condition2),each = 3))

#整理countData：将两者数据分别进行查找录入

{

  tpm_DES1 <- tpm[,grep(condition1, colnames(tpm_input))]

  #新增一列Geneid进行合并

  tpm_DES1$Geneid <- rownames(tpm_DES1)

  tpm_DES2 <- tpm[,grep(condition2, colnames(tpm_input))]

  tpm_DES2$Geneid <- rownames(tpm_DES2)

  #将两个数据框进行合并，以Geneid为参考

  tpm_DES <- merge(tpm_DES1,tpm_DES2, by = "Geneid")

  rownames(tpm_DES) <- tpm_DES[,1]

  #去除Geneid列

  tpm_DES <- tpm_DES[,-1]

  #将在各转录组中均为0的基因去掉

  tpm_DES <- tpm_DES[rowSums(tpm_DES) != 0,]

}

# 创建 DESeq2 的数据对象

#countData是表达量矩阵或者数据框

#colData是使用condition作为差异表达分析的条件变量。

#在这种设计下，DESeq2将计算每个基因在不同条件下的差异表达。

dds <- DESeqDataSetFromMatrix(countData = tpm_DES,

                              colData = data.frame(condition = conditions),

                              design = ~ condition)

 

#进行差异分析：这一步耗时较长，45个转录组的count文件花费大约2小时

dds <- DESeq(dds)

#获取差异表达基因的统计结果。

res <- results(dds, cooksCutoff = F, independentFiltering = F)

 

#查看前几行

head(res)

#将res转成dataframe格式

res_dataframe <- data.frame(res)

head(res_dataframe)

#如果没有安装dyplr，则安装dyplr，否则会直接跳过

{if (!requireNamespace("dplyr", quietly = TRUE))

  install.packages("dplyr")

}

#加载dplyr包

library(dplyr)

#在res_dataframe中新建group列

#根据log2FoldChange和padj条件进行写入“up”或者“down”，将结果赋予res_dataframe2

res_dataframe %>%

  mutate(group =case_when(

    log2FoldChange >= 2 & padj <= 0.05 ~ "UP",

    log2FoldChange <= -2 & padj <= 0.05 ~ "DOWN",

    TRUE ~ "NOT_Change"

)) -> res_dataframe2

#统计group列中不同数值的情况

table(res_dataframe2$group)

#输出路径定义：

File_output_path <- "./Differential_Expression_results/"

#输出文件名定义：

File_name <- paste0(condition1, "_vs_", condition2,".csv")

 

output_file <- paste0(File_output_path, File_name)

#导出.csv文件，quote=F：禁止在文本中添加“”

write.csv(res_dataframe2,file = output_file, quote = F)

 

Loading