R语言DESeq2包差异表达,参考了1
首先需要安装DESeq2包,这个包位于bioconductor中,通过BiocManager来安装
# 安装BiocManager
install.packages("BiocManager")
# 安装DESeq2
BiocManager::install("DESeq2")
安装过程中可能需要选择源,选择一个比较近的源即可。安装失败很可能是网络问题,这时候需要手动指定源。
然后载入包
library(DESeq2)
载入有可能报错缺依赖包,根据报错信息装上去即可。
安装完成后开始进行分析,以TCGA乳腺癌数据为例,假设数据文件位于当前目录
# 载入数据, a是counts矩阵,c是样本标签
# 具体来说,a是一个以Tab分隔的矩阵,行是基因,列是样本;c是一个只有两列的矩阵,一列是样本,另一列是样本标签,样本类别。
a <- read.table("./TCGA-BRCA.txt", header=T, check.name=F, row.name=1)
c <- read.table("./TCGA-BRCA.htseq_counts_samples.txt", header=T, check.name=F, row.name=1)
进行一定的处理构建dds对象
# 列信息,第二个参数表示有一列叫group_list,值是c的label这一列
colData <- data.frame(row.names=rownames(c), group_list=c$label)
# 解析一下参数,countData是counts矩阵,colData是表达矩阵的列信息,这个矩阵的行需要跟counts矩阵的列能对应;design是一个式子,使用到的列名必须是colData里面的
dds <- DESeqDataSetFromMatrix(countData = a, colData = colData, design = ~ group_list)
dds对象大约是这样:
> dds
class: DESeqDataSet
dim: 19731 1222
metadata(1): version
assays(1): counts
rownames(19731): ENSG00000000003 ENSG00000000005 ... ENSG00000281883
ENSG00000281887
rowData names(0):
colnames(1222): TCGA-A8-A07U-01 TCGA-D8-A1JM-01 ... TCGA-BH-A1FD-11
TCGA-A7-A13F-11
colData names(1): group_list
然后进行归一化
dds2 <- DESeq(dds)
dds2对象大约是这样:
> dds2
class: DESeqDataSet
dim: 19731 1222
metadata(1): version
assays(6): counts mu ... replaceCounts replaceCooks
rownames(19731): ENSG00000000003 ENSG00000000005 ... ENSG00000281883
ENSG00000281887
rowData names(23): baseMean baseVar ... maxCooks replace
colnames(1222): TCGA-A8-A07U-01 TCGA-D8-A1JM-01 ... TCGA-BH-A1FD-11
TCGA-A7-A13F-11
colData names(3): group_list sizeFactor replaceable
取结果
res <- results(dds2)
保存结果
write.table(res, "de_result.txt", sep="\t", quote=F)
保存的是所有的结果,后续可以根据需要筛选log2FC,如果需要筛选FDR,需要自己计算,参考2。
R语言中对于一个函数的功能和参数不明白可以用help,如
help(DESeqDataSetFromMatrix)
会打开文档的网页。