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

生信技能树

发布于 2021-08-25 11:48:33

9671

发布于 2021-08-25 11:48:33

最近在这里 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的重新编号。

下次再分享。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-08-17，如有侵权请联系 cloudcommunity@tencent.com 删除

typescript