这个文件大到R语言已经无能为力 / 四六文摘

最近在这里 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 学习笔记

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

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

首先切割 GSE158055_covid19_barcodes.tsv.gz

然后根据barcode去mtx里面拆分

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

相关推荐