不走寻常路的单细胞表达量矩阵读取

生信技能树

发布于 2024-03-06 11:20:39

4110

发布于 2024-03-06 11:20:39

文章被收录于专栏：生信技能树

之前在在单细胞天地教程：表达矩阵逆转为10X的标准输出3个文件，详细介绍过 10X技术的单细胞转录组的3个标准文件，虽然说绝大部分文献提供其数据的时候并不是标准的文件名字，但是3个文件的文件名字还是通常会遵循以下模式：

Feature / Gene-Barcodes Matrix 文件：这个文件的命名通常包含了数据类型（例如基因表达量）和文件格式（例如稀疏矩阵）。一般情况下，这个文件名中可能包含 "matrix"、"gene_bc_matrix" 或类似的关键词。有时也会包含数据集的名称或样本编号。
Barcode 文件：这个文件通常命名为 "barcodes" 或者包含 "barcode" 关键词。里面有每个样品里面的每个细胞的标签信息，这个信息其实是无所谓的。
Feature / Gene ID 文件：这个文件通常命名为 "features"、"genes" 或包含 "gene" 关键词。取决于不同版本的cellranger定量结果。

需要把每个样品都整理成为3个标准文件，文件名字和文件格式如下所示：

3个标准文件

所以很容易批量读取这样的文件，代码如下所示：

dir='GSE201048_RAW/outputs/'
samples=list.files( dir )
samples 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  tmp = Read10X(file.path(dir,pro )) 
  if(length(tmp)==2){
    ct = tmp[[1]] 
  }else{ct = tmp}
  sce =CreateSeuratObject(counts =  ct ,
                          project =  pro  ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})

也有一些文献给的是txt或者csv格式的，如下所示的文件：

  10M  7 30  2019 GSM3984317_NO.1.expression_matrix.txt.gz
  4.6M  7 30  2019 GSM3984318_NO.2.expression_matrix.txt.gz
  5.3M  7 30  2019 GSM3984319_NO.3.expression_matrix.txt.gz
  7.3M  7 30  2019 GSM3984322_NO.7.expression_matrix.txt.gz
   13M  7 30  2019 GSM3984326_NO.12.expression_matrix.txt.gz
  5.0M  7 30  2019 GSM3984330_NO.16.expression_matrix.txt.gz
   14M  4 24  2020 GSM4495152_NO.20.expression_matrix.txt.gz

代码也是很简单的，批量读取，如下所示：

dir='GSE135045_RAW/'
samples=list.files( dir )
samples 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  f= file.path(dir,pro);print(f)
  ct <- data.table::fread( f, data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  sce <- CreateSeuratObject(ct,project = gsub('.expression_matrix.txt.gz','',pro),
                            min.cells = 5,
                            min.features = 500) #后面就可以单细胞处理的标准流程啦
 
  return(sce)
}) 
do.call(rbind,lapply(sceList, dim))
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids = samples  )  
sce.all <- JoinLayers(sce.all)

但是我看到了一个比较狡猾的数据集（GSE133283），它官网给出来了的文件如下所示：

  32M  6 24  2019 GSM3904816_Adult-1_gene_counts.tsv.gz
   32M  6 24  2019 GSM3904817_Adult-2_gene_counts.tsv.gz
   33M  6 24  2019 GSM3904818_Adult-3_gene_counts.tsv.gz
   37M  6 24  2019 GSM3904819_Aged-1_gene_counts.tsv.gz
   37M  6 24  2019 GSM3904820_Aged-2_gene_counts.tsv.gz
   64M  6 24  2019 GSM3904821_EAE-1_gene_counts.tsv.gz
   63M  6 24  2019 GSM3904822_EAE-2_gene_counts.tsv.gz
   30M  6 24  2019 GSM3904823_Juvenile-1_gene_counts.tsv.gz
   29M  6 24  2019 GSM3904824_Juvenile-2_gene_counts.tsv.gz
   44M  6 24  2019 GSM3904825_Juvenile-3_gene_counts.tsv.gz
   40M  6 24  2019 GSM3904826_Juvenile-4_gene_counts.tsv.gz

它看起来是每个样品一个独立的文本文件，但里面并不是行列式的表达量矩阵文件，读入简单肉眼看了看：

>   f= file.path(dir,pro);print(f)
[1] "GSE133283_RAW/GSM3904816_Adult-1_gene_counts.tsv.gz"
ct <- data.table::fread( f, data.table = F)
>   head(ct)
           gene             cell count
1 0610005C13Rik AACCGCGTCCGTTGCT     1
2 0610005C13Rik ATTATCCAGTTGTCGT     1
3 0610005C13Rik CAGTAACAGCAGCCTC     1
4 0610005C13Rik CATCGAAGTAGCGTAG     1
5 0610005C13Rik CCACTACTCACATGCA     1
6 0610005C13Rik CCTAAAGCAAGTAATG     1
>   ct[1:4,1:4]
Error in `[.data.frame`(ct, 1:4, 1:4) : undefined columns selected 
>   dim(ct)
[1] 6182813       3

是稀疏矩阵的简化版，我略微思考了一下，做了一个简单的变幻：

  library(reshape2)
  tmp = dcast(ct,gene~cell)
  tmp[1:4,1:4]

可以看到，很多NA，其实就是单细胞转录组里面的0值，需要替换一下：

>   tmp[1:4,1:4]
           gene AAACCTGAGATGTGTA AAACCTGAGGTACTCT AAACCTGAGTGTTAGA
1 0610005C13Rik               NA               NA               NA
2 0610007N19Rik               NA               NA               NA
3 0610007P14Rik                1                2               NA
4 0610009B14Rik               NA               NA               NA
> dim(tmp)
[1] 21892  3664

这个时候我们可以借助于r编程语言里面的reshape2包的dcast函数进行数据转换，在 R 语言中，reshape2 包提供了 dcast() 函数，用于将数据框从长格式（long format）转换为宽格式（wide format）。长格式数据通常包含多行和少列，每行对应一个观察值，并且包含一个用于标识不同组的变量；而宽格式数据通常包含少行和多列，每行对应一个唯一的标识符，并且包含多个变量。最后的完整的代码是：

dir='GSE133283_RAW'
samples=list.files( dir )
samples 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  f= file.path(dir,pro);print(f)
  ct <- data.table::fread( f, data.table = F)
  head(ct)
  dim(ct)
  #ct[1:4,1:4]
  library(reshape2)
  tmp = dcast(ct,gene~cell)
  tmp[1:4,1:4]
  ct = tmp
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  ct[is.na(ct)]=0
  sce <- CreateSeuratObject(ct,project = gsub('_gene_counts.tsv.gz','',pro),
                            min.cells = 5,
                            min.features = 500) #后面就可以单细胞处理的标准流程啦
  
  return(sce)
})

有了这个seurat的对象，后面就是我们常规的！代码在：（链接: https://pan.baidu.com/s/1pKEnPmWXi-pTab0WZUWzgg?pwd=a7s1）相信大家很容易跟着去复现一次！