安装这些包

1
2
3
BiocManager::install(c("openxlsx",'org.Bt.eg.db','GO.db','topGO',
                     'GSEABase','clusterProfiler','fgsea',
                     'tidyr','apeglm','DESeq2','pacman'))
1
2
3
4
rm(list = ls())
pacman::p_load("openxlsx",'org.Bt.eg.db','GO.db','topGO',
               'GSEABase','clusterProfiler','fgsea',
               'tidyr','apeglm','DESeq2')#迅速加载

插入表格

image-20211215193210422

1
2
3
4
5
6
7
8
9
10
gene_list = read.xlsx("circRNA_ALL_Tpm.xlsx",sheet = 1)
colnames(gene_list)
rownames(gene_list) = gene_list$X1
gene_list = gene_list[ ,-1]
head(gene_list)
#rownames(gene_list) <- gene_list$XXX  #这里$出第一列ID，如果错误输出的表格会不完整
length(colnames(gene_list))
gene <- as.matrix(gene_list[,1:21])
head(gene)

image-20211215193438778

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
names(gene)
# fake metadata 分组分列
coldata <- data.frame(
  condition = factor(c(
    rep("Heart", 3),
    rep("Liver", 3),
    rep("Spleen", 3),
    rep("Lung", 3),
    rep("Kidney", 3),
    rep("Muscle", 3),
    rep("Fat", 3)
    )))

trait <- c('Heart','Liver','Spleen','Lung','Kidney','Muscle')#所有的组

coldata$condition <- relevel(coldata$condition, ref = "Fat")#所有的组与Fat相比，如果有需要，替换Fat成为其他即可

image-20211208155201461

计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
str(coldata)

dds <- DESeqDataSetFromMatrix(
  countData = round(gene),#这里使用了round()，取整数不然计算后的TPM值不为整数会报错
  colData = coldata,
  design= ~ condition)

dds <- DESeq(dds, betaPrior = FALSE)

result_extract = as.data.frame(resultsNames(dds)[-1])

levels(coldata$condition)

setLabels <- result_extract

循环对比输出！

1
2
3
4
5
6
7
8
9
10
11
12
13
14
res_ori <- list()

for (set in 1:length(setLabels$`resultsNames(dds)[-1]`)) 
{
  res_ori[[set]] <- as.data.frame(results(dds, name=setLabels$`resultsNames(dds)[-1]`[set]));
  #res_ori_2[[set]] <- data.frame(ENSEMBL=row.names(res_ori_1[[set]]),
                                 #res_ori_1[[set]]);
  #res_ori[[set]] <- dplyr::left_join(res_ori_2[[set]],annot,by=NULL);
  res_ori[[set]]$Change <- ifelse(res_ori[[set]]$padj < 0.05 & 
                                    abs(res_ori[[set]]$log2FoldChange) >= 1,
                                  ifelse(res_ori[[set]]$log2FoldChange > 1,'Up','Down'),
                                  'Stable');
  write.csv(res_ori[[set]],paste0(setLabels$`resultsNames(dds)[-1]`[set],"_FC_2_Padj0.05",".csv"))
}