RNA-seq数据差异表达分析
分析转录组测序数据时，通常使用p值/q值和foldchange值来衡量基因的差异的表达水平。目前，大家普遍都认为转录组数据的read counts（即基因的reads数量）符合泊松分布。几个用于差异表达分析的R包如DESeq2和edgeR等，都是基于负二项分布模型设计的，整体而言结果相差不大。Limma包也可以用来分析RNA-seq数据，但主要用于分析芯片数据，现在用的人不多了。当然如果用泊松分布来做差异表达分析的话，也存在缺点，可能会忽视生物学样本间的个体差异。
这里，我将RNA-seq数据差异表达分析大体分为差异表达基因鉴定和后续分析两个部分。
01
差异表达基因鉴定
首先准备好软件的输入数据：表达矩阵（counts/FPKM/RPKM等），文件名为count_test.txt。
具体格式如下：
1
DESeq2
DESeq2要求的输入数据是raw count，无需对数据进行标准化处理，如FPKM/TPM/RPKM等。分析的代码如下：
- ##加载DESeq2包library(DESeq2)##读取数据datacount <- read.table(file ="count_test.txt", sep ="\t", header = T, row.names =1)##构建dds对象datacount <- round(as.matrix(datacount))##分为control和cold两组groups<-c("control-1", "control-2", "cold-1", "cold-2")Data <- data.frame(condition = as.factor(groups))rownames(Data) <- colnames(datacount)dds <-DESeqDataSetFromMatrix(countData = datacount, colData = Data, design= ~condition)## 去掉所有条件都没有read的基因dds <- dds[ rowSums(counts(dds))>1,]##使用DESeq函数预估离散度dds <-DESeq(dds)res <- results(dds)##设定阈值筛选差异基因res <- res[order(res$padj), ]Gene_diff <- subset(res, padj <0.05 & (log2FoldChange >1 | log2FoldChange < -1))Gene_diff_DESeq2 <- row.names(Gene_diff)resData <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)# 写入文件write.csv(resData, file ="control_vs_cold_DESeq2.csv", row.names = F)
其中可以用多重检验校正后的padj来筛选显著差异表达的基因。鉴定完后利用p value和log2FoldChange值还可以画个火山图来直观显示。
其中输入数据"dif_exp.csv"文件具体格式和画图的代码如下所示：
火山图代码（python）：
import seaborn as snsimport numpy as npimport pandas as pd
difexp1= pd.read_table('dif_exp.csv', sep='\t')difexp = difexp1.dropna()
##将所有基因分为三部分difexp['sig'] = 'normal'difexp.loc[(difexp.fold > 1 ) & (difexp.pvalue < 0.05),'sig'] = 'up'difexp.loc[(difexp.fold < -1 ) & (difexp.pvalue < 0.05),'sig'] = 'down'difexp['log(pvalue)'] = -np.log10(difexp['pvalue'])
ax = sns.scatterplot(x="fold", y="log(pvalue)", hue='sig', hue_order = ('down','normal','up'), palette=("#377EB8","grey","#E41A1C"), data=difexp)ax.set_ylabel('-log10(pvalue)',fontweight='bold')ax.set_xlabel('log2(FoldChange)',fontweight='bold')
2
edgeR
edgeR包也是分析RNA-seq数据最常用的R包，它的input数据也是原始的gene counts。分析的代码如下：
##加载edgeR包library(edgeR)##读取数据exprSet <- read.table(file = "count_test.txt", sep = "\t", header = TRUE, row.names = 1, stringsAsFactors = FALSE)##分为两组groups <- factor(c(rep("control", 2), rep("cold", 2)), levels = c("control", "cold"))exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2, ]##构建DGEList对象exprSet <- DGEList(counts = exprSet, group = groups)exprSet <- calcNormFactors(exprSet)##估计离散度exprSet <- estimateCommonDisp(exprSet)exprSet <- estimateTagwiseDisp(exprSet)##找差异基因，并按阈值进行筛选list <- exactTest(exprSet)Gene_diff <- topTags(list, n=nrow(exprSet))##筛选差异表达基因Gene_diff_edgeR <- subset(Gene_diff$table, FDR < 0.05 & (logFC > 1 | logFC < -1))Gene_diff_edgeR <- row.names(Gene_diff_edgeR)write.csv(Gene_diff$table, file = "control_vs_cold_edgeR.csv")3
cufflinks软件包
基因表达标准化一般会计算CPM（counts per million）值，即read counts数除以总reads数再乘以1,000,000。什么RPKM，FPKM，TPM，都是基于CPM对表达值进行归一化。之前有人发现用cuffdiff计算筛选出的一些差异表达基因其实在样本间差异并不显著，但不知怎么地会计算出一个显著的p value值，这也是现在很多人弃用cuffdiff的一个重大原因吧。
02
后续分析
在筛选得到差异表达基因list后，一般接下来就需要分析这些基因参与了哪些功能，代表着怎样的生物学意义。我们可以做：
01
富集分析
包括GO和KEGG富集分析，可以用R里的clustProfiler包进行，也可以利用已有的一些网站。
02
聚类分析
基于差异表达基因间表达模式的相似性，分为不同类，后续可以对这些不同类的基因list分别进行分析。
03
基因共表达网络分析（WGCNA）
基因共表达网络是基于基因间表达模式的相似性构建的网络。通过构建基因共表达网络，可以深入地研究基因间的相互关系并挖掘关键途径中的关键功能模块或核心基因。