ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
106 stars 30 forks source link

读取不同格式的单细胞转录组数据及遇到问题的解决办法 #3481

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/mpEQU_aKcaq3cbzh8JsSBw

ixxmu commented 1 year ago

读取不同格式的单细胞转录组数据及遇到问题的解决办法 by 单细胞天地

几种不同格式单细胞转录组数据的读取方式,以及最近在读取单细胞转录组数据的时候遇到的问题。

读取h5格式的单细胞文件

以下面这个数据集为例,h5格式的单细胞文件 GSE175687

GEO Accession viewer (nih.gov)

批量导入数据:

rm(list=ls())
options(stringsAsFactors = F
library(Seurat)
library(dplyr)
library(hdf5r)

library(data.table)
dir='GSE175687_raw/' 
samples=list.files( dir )
samples

sceList = lapply(samples,function(pro){ 
  print(pro) 
  sce =CreateSeuratObject(counts =  Read10X_h5( file.path(dir,pro)) ,
                          project =   gsub('_filtered_feature_bc_matrix.h5','',gsub('^GSM[0-9]*_','',pro) ) ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})
sceList
samples

整合数据

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =   gsub('_filtered_feature_bc_matrix.h5','',gsub('^GSM[0-9]*_','',samples)))

读取10X格式的单细胞文件

文件输入格式

10X格式文件结构:

├── [ 160]  Sample1
│   ├── [2.2M]  barcodes.tsv.gz
│   ├── [221K]  features.tsv.gz
│   └── [4.5M]  matrix.mtx.gz
├── [ 160]  Sample2
│   ├── [2.2M]  barcodes.tsv.gz
│   ├── [221K]  features.tsv.gz
│   └── [4.0M]  matrix.mtx.gz

可以看到,每个10X的样品都是一个独立的文件夹,每个文件夹里面都是 matrix.mtx, genes.tsv (or features.tsv), and barcodes.tsv  这样的3个文件,非常有规律。

library(data.table)
sce =CreateSeuratObject(counts =  Read10X("./GSE135710_raw/") ,
                          #project = gsub('.gz','',gsub('^GSM[0-9]*_','') ) ,
                          min.cells = 5,
                          min.features = 300 )

library(stringr)
head(rownames(sce@meta.data))
phe=str_split(rownames(sce@meta.data),'-',simplify = T)
head(phe)
tail(phe)
table(phe[,2])
sce@meta.data$orig.ident=paste0('p',phe[,2])
table(sce@meta.data$orig.ident)

fread 读取txt.gz格式文件

数据集:GSE183625

GEO Accession viewer (nih.gov)

导入数据

rm(list=ls())
options(stringsAsFactors = F
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)

getwd()   
library(data.table)
dir='GSE183625_RAW/' 
samples=list.files( dir )
samples

library(data.table)
sceList = lapply(samples,function(pro){ 
  print(pro)
  ct=fread(file.path( dir ,pro),data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  sce=CreateSeuratObject(counts =  ct ,
                         project = gsub('.expression_matrix.txt.gz','',gsub('^GSM[0-9]*_','',pro) ) ,
                         min.cells = 5,
                         min.features = 300)
  
  return(sce)
})
sceList
samples

整合数据

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =   gsub('.expression_matrix.txt.gz','',gsub('^GSM[0-9]*_','',samples) )  )

fread读取单细胞测序数据csv格式

数据集:GSE129516

GEO Accession viewer (nih.gov)

导入数据:

rm(list=ls())
options(stringsAsFactors = F
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
 
library(data.table)
dir='GSE129516_RAW/' 
samples=list.files( dir )
samples

sceList = lapply(samples,function(pro){ 
  print(pro) 
  ct=fread(file.path( dir ,pro),data.table = F)
  ct[1:4,1:4]
  rownames(ct)=ct[,1]
  ct=ct[,-1]
  sce =CreateSeuratObject(counts =  ct ,
                          project =   gsub('_filtered_gene_bc_matrices.csv.gz','',gsub('^GSM[0-9]*_','',pro) ) ,
                          min.cells = 5,
                          min.features = 300 )
  return(sce)
})
sceList 
samples

整合数据:

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =   gsub('_filtered_gene_bc_matrices.csv.gz','',gsub('^GSM[0-9]*_','',samples) )  )

最近在读取单细胞转录组数据的时候遇到的问题

同样都是10X数据,在文件读取的时候也是会出现问题的,有的GEO上传上去的数据并不标准,标准的文件如下。

最近遇到的数据集是这样的,需要把第一列和第一行都切割掉。

第一行不切割的话,会在读取的时候报错,文件读不进去。而且由于feature是一列的,所以需要加个参数gene.colum=1。

第一列不切割的话,会在使用marker gene做dotplot的时候报错,因为匹配不到细胞和基因。

去head一下数据的话,发现的确是这样。

把第一列切割后就没问题了。

或者有条件的话可以自己跑一下上游定量流程,得到定量的数据。



往期回顾

单细胞转录组揭示了Semaglutide对肥胖小鼠非心肌细胞的影响

单细胞基因集合打分方法-08-UCell包-测试

细胞互作中的“avoidance ”是什么?

手把手带你复现NC图表之Figure 1

bollito—单细胞转录组的流程化工具






如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料
每天都精彩

长按扫码可关注