前面我們在 初試Seurat的V5版本 的推文里面演示了10x單細胞樣品的標準3文件的讀取，而且在使用Seurat的v5來讀取多個10x的單細胞轉(zhuǎn)錄組矩陣 的推文里面演示了多個10x單細胞樣品的標準3文件的讀取。

但是留下來了一個懸念， 就是如果我們的單細胞轉(zhuǎn)錄組并不是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》舉例說明，它的數(shù)據(jù)集是 https://www.ncbi.nlm./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個單細胞轉(zhuǎn)錄組樣品，8 patients with ATAA (4 women and 4 men) and 3 controls (2 women and 1 man). 每個樣品都是一個獨立的txt文本文件蘊藏著其表達量矩陣信息。

值得注意的是這個2020的數(shù)據(jù)集還被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對象里面的各個樣品的表達量矩陣的分開的，就會導(dǎo)致所有的后面的步驟都失敗。而它每個樣品并不是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，它里面有每個單細胞樣品的表達量矩陣，但是每個樣品的基因數(shù)量和細胞數(shù)量都是不一樣的哦。然后提前把矩陣合并之前需要首先把基因數(shù)量對齊，合并后才構(gòu)建對象：

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) 

可以看到，我這個時候做了一個處理，就是每個樣品的基因數(shù)量都對齊了，如果某個基因在某個樣品里面特有其實它不會被考慮，因為我考慮的是絕大部分基因。

因為多個樣品合并成為了一個超級大的表達量矩陣，就是 bigct 這個變量，所以后面直接針對它來使用CreateSeuratObject函數(shù)去構(gòu)建Seurat對象，就是完美的下游分析的輸入數(shù)據(jù)啦。

第二種方法是把矩陣還原成為10x的3文件

前面我們指出來了，它每個樣品并不是10x單細胞樣品的標準3文件，每個樣品都是一個獨立的txt文本文件蘊藏著其表達量矩陣信息，所以沒辦法使用前面的策略。但是，我們其實可以根據(jù)這個txt文件去把它還原成為10x的3文件，早在2020-03-16其實我就寫個一個簡單的筆記：表達矩陣逆轉(zhuǎn)為10X的標準輸出3個文件，但是那個時候的代碼略微有點麻煩，我們其實可以把它寫成一個函數(shù)，接下來讓我們演示一下吧。

每個樣品的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，它里面有每個單細胞樣品的表達量矩陣，但是每個樣品的基因數(shù)量和細胞數(shù)量都是不一樣的哦。接下來我們構(gòu)造一個自定義函數(shù)，把表達量矩陣轉(zhuǎn)為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 )
}

比較簡單，接下來就針對前面的表達量列表去循環(huán)使用這個函數(shù)即可，如下所示：

 
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('../')
  })

說實話，函數(shù)運行效率確實有點低，不過沒關(guān)系，反正是練習(xí)的代碼，我們大概是還是會選擇前面的矩陣合并的方式，并不需要把表達量矩陣轉(zhuǎn)為10x的3個文件。成功后可以看到如下所示的文件夾結(jié)構(gòu)：

│   ├── [ 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,) 

如下所示的文件夾架構(gòu)哦：

同樣的，只需有了sce.all對象，后面的降維聚類分群就是我們之前的代碼即可啦。明天中午直播一下這個全部的流程哈！