转录组数据提取技巧汇总

最新推荐文章于 2024-07-20 20:04:43 发布

凑齐六个字吧

最新推荐文章于 2024-07-20 20:04:43 发布

阅读量594

点赞数 11

分类专栏：技能树学徒作业文章标签：数据库服务器

本文链接：https://blog.csdn.net/zfyyzhys/article/details/139219162

版权

技能树学徒作业专栏收录该内容

6 篇文章 0 订阅

订阅专栏

数据整理是数据分析的前提，而数据的正确提取则是数据整理的前提。

这次曾老师交代了一个作业，希望对两个GEO转录组的数据(GSE216943, GSE236049)进行差异分析。

从GEO数据库中可以发现，这两个数据集都需要下载研究者整理好的原始数据进行整合。

GSE216943数据集：

GSE236049数据集：

针对GSE216943数据集，我仅下载了count数据，然后循环读取(如下是其中一个count文件的样式，里边的列名是id和样本名称)

然而针对GSE236049数据集，我发现研究者把文件整理的很棒，但是对于生信菜鸟的我有点不知道怎么合理的读取并处理信息(如下是其中一个文件，里边有基因名称，FPKM,TPM等值)。

没办法，只好去向曾老师求助了hhhh 。曾老师立刻发给我了一份其他高手处理的流程6666666。最后这位高手的代码也解决了困扰我的问题。

为此，我也顺便把多种提取转录组数据的方法整理了一下，以供大家参考。

转录组数据提取技巧汇总

提取数据前，我们都需要注意获得的数据是count还是tpm数据(首选count数据），只有这两种格式的可以进行差异分析。其他类型的数据如log2cpm，fpkm，rpkm在进行差异表达分析时，需要重新修正为count或者tpm数据。

情况1.可以通过get_count_txt函数获得count数据

library(tinyarray)
get_count_txt("GSE242333")
exp = read.table("GSE242333_raw_counts_GRCh38.p13_NCBI.tsv.gz",check.names = F,row.names = 1,header = T)

复制网址就可以自动下载数据，之后读取即可。

情况2. 需自行下载数据，数据已经被研究者整合好了

#下载后直接读取
#GSE121949
a = read.delim("GSE121949_RNAseq_Raw_and_RPKM.txt")
#GSE235017
dat = data.table::fread("GSE235017_ref_trans_full_table.txt",check.names = F, header = T)

情况3. 原始数据可直接循环读取合并

file_directory = "~/desktop/training/GSE216943_RAW"
fs=list.files(file_directory)
dat = do.call(cbind,lapply(fs, function(x){
  read.table(file.path(file_directory,x),
             header = T,sep = '\t',row.names = 1,quote = "")
}))
dat[1:4,1:4]

GSE216943数据集中的其中一个样本文件

情况4. Affymetrix的CEL文件读取合并及处理

#BiocManager::install("affy")
library(affy)
dir_cels='~/Desktop/GSE29330/GSE29330_RAW'
affy_data = ReadAffy(celfile.path=dir_cels)
eSet = rma(affy_data) #rma函数去normalization
##要进行过滤
calls = mas5calls(affy_data) # get PMA calls
calls = exprs(calls)
absent = rowSums(calls == 'A') # how may samples are each gene 'absent' in all samples
absent = which(absent == ncol(calls)) # which genes are 'absent' in all samples
eSet = eSet[-absent,] # filters out the genes 'absent' in all samples

情况5. 数据需要自行下载，并且文件内容中没有样本信息

#代码来自 艾迪
# 读取文件信息
exp.dirs<-list.files("data/GSE236049_RAW/",full.names = T) #读取完整文件及子文件路径信息
exp.names<-list.files("data/GSE236049_RAW/") #读取子文件的名称信息(区别在于full.name是否为T)
#提取样本名称信息
#需要提取的名称信息应根据实际情况而定
DirInfor <- stringr::str_split(exp.names,"_",n = Inf, simplify = T)%>%
            as.data.frame()
#清洗raw文件并合并
expList <- mapply(function(sampleDir, sampleName) {
    df1 <- data.table::fread(sampleDir, data.table = F) %>%  
    dplyr::select("Gene Name", "FPKM", "TPM") %>% #选择需要的列
    dplyr::rename(gene = "Gene Name")%>% #重命名列
    dplyr::filter(!is.na(gene))%>% #过滤掉基因名为空的行
    dplyr::filter(!gene == "-")%>%  #过滤掉基因名为 "-" 的行
    dplyr::filter(!duplicated(gene)) #移除基因名重复的行
    names(df1)[2:3] <- paste0(sampleName, "_", names(df1)[2:3]) #修改列名
    return(df1) }, #返回处理后的列表
    sampleDir = exp.dirs, # 输入参数
    sampleName = DirInfor$V1, # 输入参数
    SIMPLIFY = F) #结果不简化为矩阵或数组，而是保持为列表

#获取FPKM/TPM数据
exps <- Reduce(function(x, y) merge(x, y, all = TRUE), expList) 
#获取FPKM数据
fpkm <- exps%>%
        tibble::as_tibble()%>%
        tibble::column_to_rownames("gene")%>% #将tibble的gene列转换为行名，并从数据中移除这一列
        dplyr::select(contains("FPKM"))%>% #选择包含FPKM的列
        dplyr::rename_all(~ gsub("_FPKM", "", .))

#获取TPM数据
tpm <- exps%>%
       tibble::as_tibble()%>%
       tibble::column_to_rownames("gene")%>%
       dplyr::select(contains("TPM"))%>%
       dplyr::rename_all(~ gsub("_TPM", "", .))

#数据检查
tpm%>%colSums()%>%table()

类似于GSE236049数据集提供的原始数据样式

情况6. 如果样本量大，原始文件中存在不想提取的数据

#代码来自 艾迪
count.dirs<-list.files("data/GSE216943_RAW/","counts",full.names = T) #在目录下搜索包含count的文件
fpkm.dirs<-list.files("data/GSE216943_RAW/","FPKM",full.names = T)
#读取count数据
count <- lapply(count.dirs, function(x) df<-data.table::fread(x))%>%
         as.data.frame()%>%
         tibble::column_to_rownames("id")%>%
         dplyr::select(contains("NC")|contains("LPS"))
#读取fpkm数据
fpkm<-lapply(fpkm.dirs, function(x) df<-data.table::fread(x))%>%
  as.data.frame()%>%
  tibble::column_to_rownames("id")%>%
  dplyr::select(contains("NC")|contains("LPS"))
#fpkm转换成TPM
fpkmToTpm <- function(fpkm){(fpkm/sum(fpkm))*10^6} #定义转换函数
tpm <- apply(fpkm,2, fpkmToTpm) #2的含义是针对每个样本
#数据检查
tpm%>%colSums()%>%table()

类似于GSE216943数据集提供的原始数据样式，假设数据集中的样本量很大的时候，我们手动勾选下载会比较麻烦。