这个文件大到R语言已经无能为力

Posted on by

最近在这里 COVID19 相关单细胞文献的数据集，看到了一个迄今为止数据量最大的。题为“COVID-19 immune features revealed by a large-scale single cell transcriptome atlas”的研究性论文，通过对196例新冠肺炎病人284个样本进行单细胞转录组测序，绘制了新冠肺炎病人的免疫图谱。我看了看，作者们还提供了一个网页工具，里面有全部的数据，简单下载了一下，发现是python格式文件，如下所示：

58G 2月 21 13:39 COVID19_ALL.h5ad
 14G 8月 12 17:04 COVID19_ALL.h5ad.tar.gz
 52G 8月 13 11:15 COVID19_ALL.h5seurat

本来是准备对它进行很简单的转换，然后愉快的进入rstudio处理它，如下所示的代码：

Convert("covid19.h5ad", dest = "h5seurat", overwrite = TRUE)
so <- LoadH5Seurat("covid19.h5seurat")

尴尬的是，报错了，搜索了一下报错信息发现了这个 ，https://github.com/satijalab/seurat/issues/4030

没办法，重新看了看文献， 原来是有gse号啊， 重新下载如下所示文件：


9.0M 5月 17 13:09 GSE158055_cell_annotation.csv.gz
5.8M 5月 17 13:09 GSE158055_covid19_barcodes.tsv.gz
 20G 8月 13 16:41 GSE158055_covid19_counts.mtx
7.7G 5月 17 13:48 GSE158055_covid19_counts.mtx.gz
 71K 5月 17 13:09 GSE158055_covid19_features.tsv.gz

 49K 5月 17 13:09 GSE158055_sample_metadata.xlsx

这个时候就很熟悉了，标准的 barcodes.tsv.gz和counts.mtx.gz，和features.tsv.gz，足够走我们的单细胞转录组流程啦！

但是悲剧的是，仍然是报错，如下所示：

library(Seurat)
# https://hbctraining.github.io/scRNA-seq/lessons/readMM_loadData.html
library(Matrix)
library(data.table)

mtx=readMM('GSE158055_covid19_counts.mtx.gz')
mtx
dim(mtx)

cl=fread('GSE158055_covid19_barcodes.tsv.gz',
 header = T,data.table = F ) 
head(cl)
rl=fread('GSE158055_covid19_features.tsv.gz',
 header = T,data.table = F ) 
head(rl)

dim(mtx)
colnames(mtx) <- rl$V1
rownames(mtx) <- cl$V1
cov19=CreateSeuratObject(counts = t(mtx),
 pro='cov19')
head(cov19@meta.data)

再次搜索了报错信息，:https://github.com/satijalab/seurat/issues/4030

I want to load a large-scale single cell data using Read10X. But there is an error message as follows:

Error in scan(file, nmax = 1, what = what, quiet = TRUE, ...) :
scan() expected 'an integer', got '2319108599'

Do you have any idea to solve this problem?

发现居然是文件太大了，超出了R语言本身的限制，真尴尬啊！

没办法，我只能说 根据 GSE158055_cell_annotation.csv.gz 里面的细胞亚群信息，进行拆分，总共是 1,462,702 cells ，首先分成6个大亚群：

  • B cells (MS4A1, CD79A, CD79B),
  • myeloid cells (CST3, LYZ),
  • NK cells (GNLY, NKG7, TYROBP)
  • epithelial cells (KRT18, KRT19),
  • CD4 and CD8 T cells (CD3D, CD3E, CD3G, CD40LG, CD8A, CD8B).

都是耳熟能详的基因啦。

使用shell命令简单统计大类细胞亚群

 zcat GSE158055_cell_annotation.csv.gz |cut -d"," -f4|sort |uniq -c|sort -n -rk1,1
 403700 B
 374454 CD8
 294788 Mono
 260141 CD4
 59062 NK
 21471 Macro
 13684 DC
 13146 Plasma
 13056 Mega
 5652 Epi
 3531 Neu
 17 Mast

如果是 继续细分：

 zcat GSE158055_cell_annotation.csv.gz |cut -d"," -f3|sort |uniq -c


 227948 B_c01-TCL1A
 92913 B_c02-MS4A1-CD27
 68283 B_c03-CD27-AIM2
 14556 B_c04-SOX5-TNFRSF1B
 11330 B_c05-MZB1-XBP1
 1816 B_c06-MKI67

 563 DC_c1-CLEC9A
 10254 DC_c2-CD1C
 186 DC_c3-LAMP3
 2681 DC_c4-LILRA4

 25 Epi-AT2
 507 Epi-Ciliated
 1538 Epi-Secretory
 3582 Epi-Squamous

 11221 Macro_c1-C1QC
 6727 Macro_c2-CCL3L1
 1030 Macro_c3-EREG
 1282 Macro_c4-DNAJB1
 671 Macro_c5-WDR74
 540 Macro_c6-VCAN
 17 Mast
 13056 Mega

 41222 Mono_c1-CD14-CCL3
 84402 Mono_c2-CD14-HLA-DPB1
 136158 Mono_c3-CD14-VCAN
 8721 Mono_c4-CD14-CD16
 24285 Mono_c5-CD16


 1477 Neu_c1-IL1B
 543 Neu_c2-CXCR4(low)
 297 Neu_c3-CST7
 376 Neu_c4-RSAD2
 59 Neu_c5-GSTP1(high)OASL(low)
 779 Neu_c6-FGF23


 57449 NK_c01-FCGR3A
 966 NK_c02-NCAM1
 647 NK_c03-MKI67


 107008 T_CD4_c01-LEF1
 16693 T_CD4_c02-AQP3
 12123 T_CD4_c03-ITGA4
 20463 T_CD4_c04-ANXA2
 20199 T_CD4_c05-FOS
 25174 T_CD4_c06-NR4A2
 16907 T_CD4_c07-AHNAK
 5080 T_CD4_c08-GZMK-FOS_h
 8466 T_CD4_c09-GZMK-FOS_l
 323 T_CD4_c10-IFNG
 12165 T_CD4_c11-GNLY
 13102 T_CD4_c12-FOXP3
 2247 T_CD4_c13-MKI67-CCL5_l
 191 T_CD4_c14-MKI67-CCL5_h


 74146 T_CD8_c01-LEF1
 27431 T_CD8_c02-GPR183
 41662 T_CD8_c03-GZMK
 18724 T_CD8_c04-COTL1
 62454 T_CD8_c05-ZNF683
 14480 T_CD8_c06-TNF
 61358 T_CD8_c07-TYROBP
 30489 T_CD8_c08-IL2RB
 18303 T_CD8_c09-SLC4A10
 6463 T_CD8_c10-MKI67-GZMK
 517 T_CD8_c11-MKI67-FOS
 1181 T_CD8_c12-MKI67-TYROBP
 1097 T_CD8_c13-HAVCR2
 16149 T_gdT_c14-TRDV2

还是蛮复杂的！

首先切割 GSE158055_covid19_barcodes.tsv.gz

因为细胞注释文件里面的1,462,702 cells 都是有唯一的细胞barcodes，而且注释到了不同的细胞亚群，所以可以先拆分出来全部的 细胞barcodes信息，如下所示：

zcat GSE158055_cell_annotation.csv.gz > GSE158055_cell_annotation.csv

awk -F',' '{print NR","$0 >> ("GSE158055_cell_annotation_"$4)}' GSE158055_cell_annotation.csv

awk -F',' '{print $1 >> ("GSE158055_cell_barcodes_"$4)}' GSE158055_cell_annotation.csv

awk -F',' '{print NR >> ("GSE158055_cell_number_"$4)}' GSE158055_cell_annotation.csv

可以看到全部的：

 wc -l GSE158055_cell_annotation_*
 403700 GSE158055_cell_annotation_B
 260141 GSE158055_cell_annotation_CD4
 374454 GSE158055_cell_annotation_CD8
 13684 GSE158055_cell_annotation_DC
 5652 GSE158055_cell_annotation_Epi
 21471 GSE158055_cell_annotation_Macro 
 17 GSE158055_cell_annotation_Mast
 13056 GSE158055_cell_annotation_Mega
 294788 GSE158055_cell_annotation_Mono
 3531 GSE158055_cell_annotation_Neu
 59062 GSE158055_cell_annotation_NK
 13146 GSE158055_cell_annotation_Plasma

同时

$ wc -l GSE158055_cell_barcodes_*
 403700 GSE158055_cell_barcodes_B
 260141 GSE158055_cell_barcodes_CD4
 374454 GSE158055_cell_barcodes_CD8
 13684 GSE158055_cell_barcodes_DC
 5652 GSE158055_cell_barcodes_Epi
 21471 GSE158055_cell_barcodes_Macro 
 17 GSE158055_cell_barcodes_Mast
 13056 GSE158055_cell_barcodes_Mega
 294788 GSE158055_cell_barcodes_Mono
 3531 GSE158055_cell_barcodes_Neu
 59062 GSE158055_cell_barcodes_NK
 13146 GSE158055_cell_barcodes_Plasma

有了不同细胞亚群的各自的细胞barcodes就可以去拆分counts.mtx 啦。

然后根据barcode去mtx里面拆分

有意思的是，下面的代码报错了！

ls GSE158055_cell_number_* |while read id
do
id=GSE158055_cell_number_Mast
cat $id GSE158055_covid19_counts.mtx |perl -alne '{$h{$F[0]}=1;print if exists $h{$F[1]} }' >${id%%.*}.mtx
done

我仔细看了看，因为这个20G的GSE158055_covid19_counts.mtx 文件并不是制表符分割，是普通的空格分割。

而且就算是我切开了，还需要修改 counts.mtx 文件的表头3行信息，以及细胞barcodes的重新编号。

下次再分享。

其实最简单的办法就是信息scanpy

详见《生信菜鸟团》公众号的：scanpy 学习笔记

​

Comments are closed.