python基因差异分析_TCGA的差异基因分析

在分析TCGA数据库里的RNA-seq数据之前，先了解一下TCGA样本的id名称里的小秘密：参考文章：TCGA的样本id里藏着分组信息

根据文章里的内容，我查看了前一篇文章里我下载的count文件(利用R包TCGAbiolinks进行各种数据下载)，打开是这样的：

列是样品名称，格式举例：TCGA.BA.A4IF.01A.11R.A266.07，这里面包含分组信息，比如说这个样品是肿瘤样品还是正常组织。分组信息是在这个id的第14-15位，01-09是tumor，10-29是normal。比如第一个样品，是01A，说明这个样品是个肿瘤样品。

OK,知道了哪里是分组信息，就可以开始进行操作了。

(一)准备R包

if(!require(stringr))install.packages('stringr')

if(!require(ggplotify))install.packages("ggplotify")

if(!require(patchwork))install.packages("patchwork")

if(!require(cowplot))install.packages("cowplot")

if(!require(DESeq2))BiocManager::install('DESeq2')

if(!require(edgeR))BiocManager::install('edgeR')

if(!require(limma))BiocManager::install('limma')

(二)载入数据

这里的数据是我上一篇文章下载好的count文件，不知道怎么下的请移步：利用R包TCGAbiolinks进行各种数据下载

#加载数据

> rm(list = ls())

> a= read.csv("TCGAbiolinks_HNSC_counts.csv")

> dim(a)

[1] 56512 547

#另外要关注一下，你的这个表达矩阵的第一列是样品还是基因名，如果第一列是基因名，要把第一列设置成为行名，不然只有的差异分析会出错

#将第一列换成行名

> row.names(a)

> a

NOTE：如果你没有check你的矩阵，导致了你的矩阵第一列是基因名，后面DESeq2运行会显示：“Error in DESeqDataSet(se, design = design, ignoreRank) : some values in assay are negative”这样的报错。

(三)将样品分组

> group_list

> group_list

> table(group_list)

group_list

normal tumor

44 502

(四)差异分析

DESeq2方法做差异分析

> library(DESeq2)

> colData

condition=group_list) #condition是你的分组信息

> dds

countData = a,

colData = colData,

design = ~ condition)

> dds

estimating size factors

estimating dispersions

gene-wise dispersion estimates

mean-dispersion relationship

final dispersion estimates

fitting model and testing

-- replacing outliers and refitting for 3594 genes

-- DESeq argument 'minReplicatesForReplace' = 7

-- original counts are preserved in counts(dds)

estimating dispersions

fitting model and testing

# 两两比较

> res

> resOrdered

> DEG

> head(DEG)

baseMean log2FoldChange lfcSE stat pvalue padj

ENSG00000231887 233.5976 -7.414328 0.3208937 -23.10525 4.100424e-118 1.594778e-113

ENSG00000107159 1962.6861 5