使用DEseq2循环做多组间差异表达分析

    当有多组RNA-seq数据时，有时需要对多个组合进行差异表达分析，例如当我有CIM0/CIM7/CIM14/CIM28四组时，我需要得到每个组合间的差异表达情况，CIM7:CIM0; CIM14:CIM0; CIM14:CIM7等。使用ANOVA的方式也可以进行多组间比较，但由于ANOVA是指定同一个CK，并且不能得到具体是哪组相对于CK有差异表达，不能精准的解决我的需求，因此选择使用DEseq2循环对不同组进行差异表达分析。

一. R脚本

    目前脚本中DEGs（差异表达基因）筛选标准为log2FoldChange>1或log2FoldChange<-1以及pvalue<0.05, qvalue<0.05，可根据需求更改阈值。

###

library(GenomicFeatures)

library(DESeq2)

library(dplyr)

###load and set output file

file_in <- "counts.txt"

file_design <- "design.txt"

file_compare <- "compare2.txt"

file_deg_num = paste("data//","DE_",file_in,sep="")  ##number of DEGs in all compare

file_final_csv  = paste("data//","DE_",file_in, "_Final_Out.csv",sep="")

file_final_genelist = paste("data//","DEG_geneid_allcomapre.txt",sep="")

################################################################

###########get DEGs in different compare

# counts file

data_in = read.table(file_in, head=TRUE,row.names =1, check.names = FALSE)

mycompare=read.table(file_compare,head=TRUE)

mydesign=read.table(file_design, head = TRUE)

##filter counts

countData = as.data.frame(data_in)

dim(countData)

mycounts_filter <- countData[rowSums(countData) != 0,]

dim(mycounts_filter)

##check group

head(mycompare)

head(mydesign)

total_num=dim(mycompare)[1]  ##get the num of groups

tracking = 0

gene_num_out = c()

pvalue_cut = 0.01

condition_name = c()

for (index_num in c(1:total_num)){

  tracking = tracking + 1

  test_name = as.character(mycompare[index_num,1])

  group1 = as.character(mycompare[index_num, 2])

  group2 = as.character(mycompare[index_num, 3])

  sample1 = mydesign[mydesign$Group == group1,]

  sample2 = mydesign[mydesign$Group == group2,]

  allsample = rbind(sample1, sample2)

  counts_sample = as.character(allsample$counts_id)

  groupreal = as.character(allsample$Group)

  countData = mycounts_filter[, counts_sample]

  colData = as.data.frame(cbind(counts_sample, groupreal))

  names(colData) = c("sample", "condition")

  dds <- DESeqDataSetFromMatrix(countData = countData,

                                colData = colData,

                                design = ~ condition)

  dds <- DESeq(dds)

  resSFtreatment <- results(dds, cooksCutoff =FALSE, contrast=c("condition",group2,group1))

  out_test = as.data.frame(resSFtreatment)

  final_each = cbind( out_test$log2FoldChange,out_test$pvalue,out_test$padj)

  rownames(final_each) = resSFtreatment@rownames

  names = c('(logFC)','(pvalue)','(Qvalue)')

  final_name = paste(test_name,'_',names,sep="")

  colnames(final_each) = final_name

  if (tracking == 1){

    final_table = final_each

  }else{

    final_table = cbind(final_table, final_each)

  }

  gene_sel = out_test[((!is.na(out_test$pvalue))&(!is.na(out_test$log2FoldChange)))&out_test$pvalue < 0.05 & abs(out_test$log2FoldChange) >= 1 & out_test$padj<0.05, ]

  gene_sel<-na.omit(gene_sel)  ##delete rows contain NA

  gene_sel_up = gene_sel [gene_sel$log2FoldChange>0,]

  gene_sel_do = gene_sel [gene_sel$log2FoldChange<0,]

  file_out_up = paste("data//DEGsid//","UP_",test_name,".txt",sep="")

  file_out_do = paste("data//DEGsid//","DOWN_",test_name,".txt",sep="")

  gene_list_up = rownames(gene_sel_up)

  gene_list_do = rownames(gene_sel_do)

  all_ub_down = list(gene_list_up, gene_list_do)

  nameup <- paste(test_name,"_up",sep="")

  namedown <- paste(test_name,"_down",sep="")

  names(all_ub_down) = c(nameup, namedown)

  if (tracking == 1){

    final_genelist = all_ub_down

  }else{

    final_genelist = c(final_genelist, all_ub_down)

  }

  gene_num_up = length(gene_list_up)

  gene_num_do = length(gene_list_do)

  write.table(gene_sel_up, file= file_out_up , row.names = TRUE,col.names = TRUE)

  write.table(gene_sel_do, file= file_out_do , row.names = TRUE,col.names = TRUE)

  condition_name = c(condition_name, paste("UP_",test_name,sep=""),paste("DO_",test_name,sep=""))

  gene_num_out = c(gene_num_out,gene_num_up,gene_num_do)

}

final_DEGs_list<-do.call(cbind, lapply(lapply(final_genelist, unlist), `length<-`, max(lengths(final_genelist))))

final_DEGs_list

out_final2 = cbind(condition_name, gene_num_out)

colnames(out_final2) = c("Tests", "DEG number")

write.table(final_DEGs_list, file=file_final_genelist, row.names = FALSE, sep="\t", na = "") ###all DEGs list

write.table(out_final2, file=file_deg_num, row.names = FALSE, sep="\t")    ##DEG number in different compare

write.csv(final_table, file=file_final_csv, row.names = TRUE,quote=TRUE)  ##allgene in all compare

二. 输入文件

2.1.counts.txt

    DEseq2的输入counts数据需为整数；如果是RSEM结果，使用expected

counts取整。

2.2.design.txt

文件共四列，格式如下：

sample_id: 样本id;

counts_id: counts文件中的样品id;

Group: 样本分组；

rep: 样本重复;

2.3.compare.txt

指定比较组

Test：进行比较的组

denominator: 对照组

numerator: 比较组

注：对照组在前，比较组在后，假如做CIM7相对于CIM0的差异表达基因分析，denominator为CIM0，numerator为CIM7。

三.输出文件

    运行脚本前在当前目录下建立data文件夹，并在data下建立DEGsid文件夹，用于存放输出文件。脚本运行完成后会在data文件夹下输出3个文件：

3.1.DE_counts.txt_Final_Out.csv

输出总表，格式如下：

第一列为基因id，包括所有基因，后几列为所有比较组的logFC,pvalue以及qvaue值。

3.2.DE_counts.txt

所有比较组中差异表达基因的数目。

3.3. DEG_geneid_allcomapre.txt

    所有比较组中差异表达基因的基因id，该文件用于筛选不同组间的重叠基因。“-up”表示比较组内上调的基因，“_down”表示比较组内下调的基因。

3.4.DEGsid

    DEGsid文件夹下生成一系列文件，为每个比较组的DEseq2的分析结果, “UP_”和“DOWN”分别表示上调和下调的基因。

四.找不同比较组间的重叠DEGs

得到这些输出结果后我还想找到不同组间的重叠DEGs:

names(mydegs)获得不同比较组的名称，指定mycom1和mycom2后使用intersect（）找到交集。