基于count数据的基因差异表达分析万能代码

DoubleHelix

发布于 2022-06-13 12:59:04

3.2K0

发布于 2022-06-13 12:59:04

文章被收录于专栏：生物信息云生物信息云

关于差异分析的文章中【一文就会TCGA数据库基因表达差异分析】其实有推送过，这篇文章目前为止，有近千人付费学习。

从付费上来说，该视频教程的反应是很良好的。也说明差异分析是很普通的分析了。关于差异分析没有什么难的，原理看文章【超详细的DESeq2和edgeR包的基本原理和实战案例】中就有详细的讲解，看了就能学会。

主要是该教程是介绍R教程后，为学习我R语言视频的粉丝们录制的视频，学习差异分析的同时顺便巩固R语言，所以该教程介绍处理TCGA数据集的时候，很多R基础函数的应用都在里面。时间久了，数据处理方面过时了，因为TCGA数据库发生重大更新，所以教程在处理数据部分已经不适用。当然，数据处理部分我也写了数据库更新后，我也分享了数据库更新后相关数据的处理的代码。可以看看下面的文章：

2022-TCGA数据库重大更新后miRNA-Seq数据的下载与整理。

2022-TCGA数据库重大更新后RNASeq的STAR-Counts数据的下载与整理

2022-TCGA数据库重大更新后3行代码提取simple nucleotide variation的数据

你得到的数据，可以根据之前的文章的代码稍微改改，也是可以进行差异表达分析的。当然，这里，我把处理数据完全封装成一个函数，方便大家使用。

getTCGA_RNAseq_data = function(dataPath,json,data_type){

  ###从json文件获取信息
  metadata_json <- rjson::fromJSON(file=json)
  json_info <- do.call(rbind, lapply(1:length(metadata_json), function(i){
    TCGA_Barcode <- metadata_json[[i]][["associated_entities"]][[1]][["entity_submitter_id"]]
    json_File_Info <- data.frame(TCGA_Barcode = TCGA_Barcode)
    rownames(json_File_Info) <- metadata_json[[i]][["file_name"]]
    return(json_File_Info)
  }))
  ##读取转录组数据
  filepath = dir(path = dataPath,
                 pattern = "counts.tsv$",
                 full.names = T,
                 recursive = T)
  if(length(filepath)!=0){
    exp <- lapply(filepath,function(wd){
      tempPath <- unlist(strsplit(wd,"/"))
      filename <- tempPath[length(unlist(strsplit(wd,"/")))]
      message(paste0("微信公众号:生物信息云 提示:正在读入文件:\n",filename))
      oneSampExp <- read.table(wd,comment.char = "#",header = T,sep = "\t")
      oneSampExp = oneSampExp[-c(1:4),]
      if(wd == filepath[1]){
        oneSampExp = oneSampExp[,c("gene_id","gene_name","gene_type",data_type)]
        colnames(oneSampExp) <- c("gene_id","gene_name","gene_type",
                                  json_info[filename,"TCGA_Barcode"])
        rownames(oneSampExp) <- oneSampExp[,"gene_id"]
      }else{
        rnm = oneSampExp[,"gene_id"]
        oneSampExp = data.frame(value = oneSampExp[,data_type])
        colnames(oneSampExp) <- json_info[filename,"TCGA_Barcode"]
        rownames(oneSampExp) <-  rnm
      }
      return(oneSampExp)
    })
    data <- do.call(cbind,exp)
    return(data)
  }else{
    message("Please check the datapath")
  }
}

dataPath：下载后的文件解压后的路径。比如我这里的数据是在："F:/BioInfoNotes/data/TCGA/data"路径下，dataPath就可以设置为：

dataPath = "F:/BioInfoNotes/data/TCGA/data"

json：json就是就是下载的json文件的完整路径(包括文件名）

json = "metadata.cart.2022-04-05.json"#当前工作文件夹

data_type参数和前面文章【2022-TCGA数据库重大更新后RNASeq的STAR-Counts数据的下载与整理】中介绍的一样，是下面中的一种。

"unstranded";
"stranded_first";
"stranded_second";
"tpm_unstranded";
"fpkm_unstranded";
"fpkm_uq_unstranded"

前面说一般我们要"fpkm_unstranded"或"tpm_unstranded"的数据。但差异分析，我们要"unstranded"的数据，就是count的数据。所以我们这里的差异分析获取的是unstranded的数据。

关于差异分析以及不是什么新鲜的的了，前面我的文章，以及其他很多人的文章都有介绍。这里我再写这个教程，就是为了随后的分析简单化。

我前面的的分析都是基于TCGA来讲解，分组是按照肿瘤和正常来进行差异分析的。前面的介绍也说过差异分析的3个包：limma、DESeq2和edgeR包。一般转录组，基于count数据的差异分析，我推荐是DESeq2和edgeR，我自己常用DESeq2，你可以两者使用后取交集。

关于差异分析没有什么难的，原理看文章【超详细的DESeq2和edgeR包的基本原理和实战案例】中就有详细的讲解，你肯定能学的会的。

我这里就是将这3种方法一起封装成一个方法。

countsDEAnalysis <- function (counts, group, comparison, method = "limma", filter = TRUE){
#代码:https://mp.weixin.qq.com/s/e05YajG5nR5HhVC1xxisOA
#
}

counts当然就是我们的表达矩阵，自己测序的也行，只要是count数据的表达矩阵就行，行名是基因名，列名是样本编号。

这里，我顺便处理一下TCGA的样本。

json = "metadata.cart.2022-04-05.json"
dataPath = "data/"
data_type = "unstranded"
data = getTCGA_RNAseq_data(dataPath,json,data_type)
head(data)[,1:3]

和前面文章【2022-TCGA数据库重大更新后RNASeq的STAR-Counts数据的下载与整理】得到的表达矩阵一样。获取的数据如下：前3列我们只要gene_name这一列，而且这一列有重复的。原因在【16-gtf文件信息提取】和【生信中各种ID转换】文章中有介绍。基础知识看文章【常用生物信息 ID的介绍】。

处理一下数据得到我们想要的数据：

library(dplyr)
exp = data[apply(data[,-c(1:3)], 1, sum)!= 0,]#删除在所有样本中不表达的基因

exp = arrange(exp,gene_name) #按照gene_name列排序

exp = exp[!duplicated(exp$gene_name),]#我直接去重了的，只要第一个

rownames(exp) = exp$gene_name
anto =  exp[,c(1:3)] #保存基因类型的信息
exp = exp[,-c(1:3)] #删除多余的列用于后续分析

head(exp)[,1:3]

group就是分组信息的一个向量，和列名一一对应。

这里的案例数据是TCGA的。我们可以使用TCGAbiolinks相关的函数来区分样本信息，

#=====================获取样本信息
library(TCGAbiolinks)
##正常组织样本
SamN <- TCGAquery_SampleTypes(barcode = colnames(exp),typesample = "NT")
##肿瘤组织样本
SamT <- setdiff(colnames(exp),SamN)
exp <- exp[,c(SamT,SamN)]#重新排序一下列的顺序
gr_meta = c(rep("Tumor",length(SamT)),rep("Normal",length(SamN)))#分组信息，和列名一一对应

我们需要gr_meta这么一个向量。

> gr_meta
 [1] "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor" 
[10] "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor" 
[19] "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor" 
[28] "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Tumor"  "Normal"
[37] "Normal" "Normal" "Normal" "Normal" "Normal" "Normal" "Normal" "Normal"

comparison就是一个差异分析中谁和谁相比的一个向量，比如我们这里：

comparison = 'Tumor-Normal'，中间用“-"链接，前面比后面。

method参数就是使用什么方法进行差异分析，DESeq2和edgeR中的一种，当然也可以是limma，只是基于count数据的差异分析，我还是建议使用DESeq2和edgeR。

filter参数就是TRUE或者FALSE，表示是否过滤数据。

然后下面就可以进行差异表达分析了。对于其他数据，你只需要整理成为前面这样的代码就可以了。

##使用DESeq2包进行差异表达分析
library(DESeq2)
library(edgeR)
DEG = countsDEAnalysis(counts = exp,
                       group= gr_meta,
                       comparison='Tumor-Normal',
                       method = "DESeq2",
                       filter = F)

筛选差异基因

DEG <- DEG %>%
  mutate(direction = factor(ifelse(FDR < 0.01 & abs(logFC) > 1.5,#添加direction一列
                                   ifelse(logFC > 1.5, "Up", "Down"),"Ns"),
                            levels=c('Up','Down','Ns')))

> head(DEG)
             symbol     baseMean      logFC     lfcSE       stat       PValue          FDR direction
5_8S_rRNA 5_8S_rRNA 1.988501e-02 -1.2378406 3.6569164 -0.3384930 7.349917e-01           NA        Ns
5S_rRNA     5S_rRNA 6.420780e-01 -0.2522589 0.9205390 -0.2740339 7.840586e-01 8.564467e-01        Ns
7SK             7SK 9.501850e-02 -1.3935564 3.1731217 -0.4391752 6.605346e-01           NA        Ns
A1BG           A1BG 3.624819e+01 -4.4617259 0.5602173 -7.9642765 1.661932e-15 4.701700e-14      Down
A1BG-AS1   A1BG-AS1 2.837107e+01 -1.8325840 0.4094068 -4.4761935 7.598562e-06 4.346123e-05      Down
A1CF           A1CF 5.911187e+03 -3.4164799 0.5522670 -6.1862827 6.159951e-10 7.378770e-09      Down

后面就是可视化了【一文就会TCGA数据库基因表达差异分析】，差异分析后的就是富集分析【可参考文章：clusterProfiler包进行KEGG,GO,GSEA富集分析，FunRich数据库：一个主要用于基因和蛋白质的功能富集以及相互作用网络分析的独立的软件工具】等了。

拓展：根据某基因在肿瘤组织样本中表达高低分组后进行差异表达分析。

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2022-05-29，如有侵权请联系 cloudcommunity@tencent.com 删除

数据库