使用Seurat的v5来读取多个不是10x标准文件的单细胞项目

生信技能树

发布于 2023-12-28 15:16:58

5891

发布于 2023-12-28 15:16:58

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

前面我们在 初试Seurat的V5版本 的推文里面演示了10x单细胞样品的标准3文件的读取，而且在使用Seurat的v5来读取多个10x的单细胞转录组矩阵 的推文里面演示了多个10x单细胞样品的标准3文件的读取。

但是留下来了一个悬念，就是如果我们的单细胞转录组并不是10x的标准3文件，而是tsv或者csv或者txt等文本文件表达量矩阵信息，就有点麻烦了。接下来我们以2020的文章：《Single-Cell Transcriptome Analysis Reveals Dynamic Cell Populations and Differential Gene Expression Patterns in Control and Aneurysmal Human Aortic Tissue》举例说明，它的数据集是 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE155468

可以看到，作者这个时候给出来的是：

GSM4704931_Con4.txt.gz 9.2 Mb
GSM4704932_Con6.txt.gz 3.0 Mb
GSM4704933_Con9.txt.gz 10.0 Mb
GSM4704934_TAA1.txt.gz 7.7 Mb
GSM4704935_TAA2.txt.gz 5.8 Mb
GSM4704936_TAA3.txt.gz 7.2 Mb
GSM4704937_TAA4.txt.gz 12.5 Mb
GSM4704938_TAA5.txt.gz 11.7 Mb
GSM4704939_TAA6.txt.gz 8.1 Mb
GSM4704940_TAA7.txt.gz 18.7 Mb
GSM4704941_TAA8.txt.gz 6.4 Mb

是11个单细胞转录组样品，8 patients with ATAA (4 women and 4 men) and 3 controls (2 women and 1 man). 每个样品都是一个独立的txt文本文件蕴藏着其表达量矩阵信息。

值得注意的是这个2020的数据集还被2023的文章引用了，感兴趣的可以去看看，Genome-wide association study of thoracic aortic aneurysm and dissection in the Million Veteran Program. Nat Genet 2023 Jul;55(7):1106-1115. PMID: 37308786

前面提到了，如果是没有样品的txt独立读取后，再merge的时候成为的Seurat对象里面的各个样品的表达量矩阵的分开的，就会导致所有的后面的步骤都失败。而它每个样品并不是10x单细胞样品的标准3文件，所以没办法使用前面的策略。

第一种方法是把每个样品的矩阵对齐

每个样品的txt仍然是独立的读取，代码如下所示：

dir='GSE155468_RAW/' 
samples=list.files( dir ,pattern = 'gz')
samples 
library(data.table)
ctList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)
  ct=fread(file.path( dir ,pro),data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  colnames(ct) = paste(gsub('.txt.gz','',pro),
                       colnames(ct) ,sep = '_')
  ct=ct[,-1] 
  return(ct)
})

上面的代码返回了 ctList 这个list，它里面有每个单细胞样品的表达量矩阵，但是每个样品的基因数量和细胞数量都是不一样的哦。然后提前把矩阵合并之前需要首先把基因数量对齐，合并后才构建对象：

lapply(ctList, dim)
tmp =table(unlist(lapply(ctList, rownames)))
cg = names(tmp)[tmp==length(samples)]
bigct = do.call(cbind,
                lapply(ctList,function(ct){ 
                  ct = ct[cg,] 
                  return(ct)
                }))
sce.all=CreateSeuratObject(counts =  bigct, 
                       min.cells = 5,
                       min.features = 300)
sce.all
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all@meta.data$orig.ident)

可以看到，我这个时候做了一个处理，就是每个样品的基因数量都对齐了，如果某个基因在某个样品里面特有其实它不会被考虑，因为我考虑的是绝大部分基因。

因为多个样品合并成为了一个超级大的表达量矩阵，就是 bigct 这个变量，所以后面直接针对它来使用CreateSeuratObject函数去构建Seurat对象，就是完美的下游分析的输入数据啦。

第二种方法是把矩阵还原成为10x的3文件

前面我们指出来了，它每个样品并不是10x单细胞样品的标准3文件，每个样品都是一个独立的txt文本文件蕴藏着其表达量矩阵信息，所以没办法使用前面的策略。但是，我们其实可以根据这个txt文件去把它还原成为10x的3文件，早在2020-03-16其实我就写个一个简单的笔记：表达矩阵逆转为10X的标准输出3个文件，但是那个时候的代码略微有点麻烦，我们其实可以把它写成一个函数，接下来让我们演示一下吧。

每个样品的txt仍然是独立的读取，代码如下所示：

dir='GSE155468_RAW/' 
samples=list.files( dir ,pattern = 'gz')
samples 
library(data.table)
ctList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)
  ct=fread(file.path( dir ,pro),data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  colnames(ct) = paste(gsub('.txt.gz','',pro),
                       colnames(ct) ,sep = '_')
  ct=ct[,-1] 
  return(ct)
})

上面的代码返回了 ctList 这个list，它里面有每个单细胞样品的表达量矩阵，但是每个样品的基因数量和细胞数量都是不一样的哦。接下来我们构造一个自定义函数，把表达量矩阵转为10x的3个文件，如下所示：

to10x <- function(ct)
  {
  write.table(data.frame(rownames(ct),rownames(ct)),file = 'features.tsv',
              quote = F,sep = '\t',
              col.names = F,row.names = F)
  write.table(colnames(ct),file = 'barcodes.tsv',quote = F,
              col.names = F,row.names = F)
  file="matrix.mtx"
  sink(file)
  cat("%%MatrixMarket matrix coordinate integer general\n")
  cat("%\n")
  cat(paste(nrow(ct),ncol(ct),sum(ct>0),"\n")) 
  sink()
  tmp=ct[1:5,1:4]
  tmp
  tmp=do.call(rbind,lapply(1:ncol(ct),function(i){
    return(data.frame(row=1:nrow(ct),
                      col=i,
                      exp=ct[,i]))
  }) )
  tmp=tmp[tmp$exp>0,]
  head(tmp)
  write.table(tmp,file = 'matrix.mtx',quote = F,
              col.names = F,row.names = F,append = T )
}

比较简单，接下来就针对前面的表达量列表去循环使用这个函数即可，如下所示：

 
lapply(samples,function(pro){ 
  # pro=samples[1] 
  pro=gsub('.txt.gz','',pro)
  print(pro)
  ct = ctList[[1]]
  dir.create(pro)
  setwd(pro)
  to10x(ct)
  setwd('../')
  })

说实话，函数运行效率确实有点低，不过没关系，反正是练习的代码，我们大概是还是会选择前面的矩阵合并的方式，并不需要把表达量矩阵转为10x的3个文件。成功后可以看到如下所示的文件夹结构：

│   ├── [ 160]  GSM4704935_TAA2
│   │   ├── [115K]  barcodes.tsv
│   │   ├── [291K]  features.tsv
│   │   └── [ 95M]  matrix.mtx
│   ├── [ 160]  GSM4704936_TAA3
│   │   ├── [115K]  barcodes.tsv
│   │   ├── [291K]  features.tsv
│   │   └── [ 95M]  matrix.mtx

值得注意的是每个样品这个时候里面的3文件其实是并没有压缩，所以很耗费空间哦。而且因为这个时候我给出来的名字是features.tsv所以如果想使用Seurat的Read10X读取，就需要把每个样品文件夹里面的3文件gz压缩一下哦！然后把每个样品的文件夹归纳整理到 outputs 文件夹里面，就可以使用如下所示的代码啦。

library(Seurat)
tmp = list.dirs('outputs')[-1]
tmp
ct = Read10X(tmp) 
sce.all=CreateSeuratObject(counts = ct  , 
                           min.cells = 5,
                           min.features = 300,)

如下所示的文件夹架构哦：