【原】拆分大规模单细胞数据

前言

去年8月份有下载单细胞转录组和免疫组测试数据集的需求， 首先从10X官网下载，发现数据量小。

恰好发现284个样本的新冠图谱文章发表，所以决定借用此篇文献的数据，拆分为10X官方的三个文件：features.tsv.gz、barcodes.tsv.gz和matrix.mtx.gz。

1、数据下载

文章标题：COVID-19 immune features revealed by a large-scale single cell transcriptome atlas
数据链接：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE158055

# 数据下载
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_sample_metadata.xlsx
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_cell_annotation.csv.gz
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_covid19_barcodes.tsv.gz
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_covid19_features.tsv.gz
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_covid19_counts.mtx.gz
wget https://ftp.ncbi.nlm./geo/series/GSE158nnn/GSE158055/suppl/GSE158055_covid19_BCR_TCR.tar.gz

由于数据量过大，矩阵文件无法读取，遇到和此文一样的问题^[1]：https:///2021/08/20210831071133930b.html

2、以样本为单位拆分转录组数据

2.1、新建文件夹与文件重命名

mkdir -p ~/Single_cells_2021/Small_works_2021/R-package_dataset/filtered_feature_bc_matrix && cd
mv ~/GSE158055_covid19_features.tsv ./
mv ~/GSE158055_covid19_barcodes.tsv ./
mv ~/GSE158055_covid19_counts.mtx ./
rename "GSE158055_covid19_" "" *.gz
rename "counts.mtx" "matrix.mtx" *

2.2、查看大矩阵文件

按样本查看注释文件：cat GSE158055_cell_annotation.csv |sed '1d'| cut -d',' -f2| sort | uniq -c
删除注释文件第一行：sed -i '1d' GSE158055_cell_annotation.csv
# 提取矩阵文件的前两行以备后用
cat matrix.mtx | head -2 > ../header2

2.3、根据注释文件，以样本为单位拆分大矩阵文件

awk -F',' '{print NR","$0 >> ("GSE158055_cell_annotation_"$2)}' GSE158055_cell_annotation.csv
awk -F',' '{print $1 >> ("GSE158055_cell_barcodes_"$2)}' GSE158055_cell_annotation.csv
awk -F',' '{print NR >> ("GSE158055_cell_number_"$2)}' GSE158055_cell_annotation.csv

2.4、准备单样本barcode文件

根据每个样本的细胞序号文件，从大矩阵文件中提取细胞基因表达矩阵信息

for id in `ls GSE158055_cell_number_*`; do awk -F" " 'NR==FNR{a[$1]; next} $2 in a' $id matrix.mtx > Sample.${id##*_}.mtx; done
rm GSE158055_cell_number_*

2.5、准备单样本matrix文件

matrix文件前三行依次为样本基因数、barcode数、UMI总数，三列依次为基因编号、barcode编号、UMI数目。

修改每个样本的细胞基因矩阵文件，shell脚本文件名为produce_final_bc_matrix.sh。

# 1 修改每个样本matrix文件中barcode编号顺序，要求从1开始，止于样本本身细胞数目，因此每一行第二列加一后都减去文件第一行的第二列值。
# 提取样本ID，保存于文件final_choose_samples.txt
cat final_choose_samples.txt | while read id; do
maximum=`cat *${id}.mtx | awk 'NR==1 {print $2; exit}'`
tmp_value=$[$maximum-1] # tmp_value=$(($maximum-1))
cat *${id}.mtx | awk -v var="$tmp_value" '{print $1, $2-var, $3}' > Sample.${id%.*}_matrix.mtx


# 2 添加每个样本matrix文件的前三行，即统计样本barcode数、UMI总数
# 27943 $cell_number $umi_counts，基因数目为固定值
cell_number=`cat Sample.${id%.*}_matrix.mtx | awk 'END {print $2; exit}'` # 细胞数目取修改后的barcode文件最后一行的第二列
umi_counts=`awk -F' ' '{sum+=$3;} END{print sum;}' Sample.${id%.*}_matrix.mtx` # UMI数目取修改后的barcode文件第三列求和
echo 27943 $cell_number $umi_counts > third_line_tmp
cat header2 third_line_tmp > header3_tmp
cat header3_tmp Sample.${id%.*}_matrix.mtx > ${id%.*}_matrix.mtx


# 3 批量为每个文件创建目录。每个样本文件夹下单独存放基因、barcode、矩阵文件及细胞注释文件，并要求重命名和压缩文件
mkdir -p ${id%_*}/filtered_feature_bc_matrix
mv GSE158055_cell_barcodes_${id} ${id}_barcodes.tsv
cat GSE158055_cell_annotation_${id} | cut -d',' -f2- > ${id}_annotation.csv
cp features.tsv.gz ${id%_*}/filtered_feature_bc_matrix
mv ${id}_barcodes.tsv ${id%.*}_matrix.mtx ${id%_*}/filtered_feature_bc_matrix
mv ${id}_annotation.csv ${id%_*}


# 4 去前缀和跨目录压缩文件
# 跨目录去文件前缀
rename "${id}_" "" ${id%_*}/filtered_feature_bc_matrix/*
rename "${id}_" "" ${id%_*}/*
done


# 跨目录压缩文件
find $(pwd) -name "S-*" -type f -exec gzip \{\} \;

# 最后删除中间文件
rm *_tmp
rm Sample.*
rm GSE158055_cell_*

3、以样本为单位拆分TCR数据

3.1 、读取TCR文件

dataset_path <- 'F:/R_Language/R_Practice/R_Packages/testPackage/data/'
tcr_vdj <- read.delim(file.path(dataset_path, 'GSE158055_covid19_tcr_vdjnt_pclone.tsv'))
dim(tcr_vdj)
table(tcr_vdj$sampleID)
tcr_vdj[1:4, 1:4]

3.2、转换文件格式

从每个细胞一行信息转为两行

tcr_vdj2 <- tcr_vdj %>% tidyr::unite(TCRA, TCRA_cgene, TCRA_vgene, TCRA_dgene, TCRA_jgene, TCRA_cdr3aa, TCRA_cdr3nt) %>%
  tidyr::unite(TCRB, TCRB_cgene, TCRB_vgene, TCRB_dgene, TCRB_jgene, TCRB_cdr3aa, TCRB_cdr3nt) %>%
  tidyr::unite(TCRAB, TCRA, TCRB, sep='.') %>%
  tidyr::separate_rows(TCRAB, sep='\\.', convert=T) %>%
  dplyr::mutate(barcode = cellBarcode) %>%
  dplyr::select(barcode, sampleID, TCRAB) %>%
  dplyr::mutate(chain = rep(c('TRA', 'TRB'), times=n()/2)) %>%
  dplyr::mutate(is_cell='TRUE') %>%
  dplyr::mutate(high_confidence='TRUE') %>%
  dplyr::mutate(full_length='TRUE') %>%
  dplyr::mutate(productive='TRUE') %>%
  tidyr::separate(TCRAB, c("c_gene", "v_gene",'d_gene','j_gene', 'cdr3', 'cdr3_nt'), sep='\\_') %>%
  # tibble::add_column(is_cell='TRUE', .after='barcode') %>%
  dplyr::select(barcode, sampleID, is_cell, high_confidence, chain, v_gene, d_gene, j_gene, c_gene, full_length, productive, cdr3, cdr3_nt)

3.3 、按样本拆分文件

# 2-3 Split by sample
# trc_list <- tcr_vdj %>% split(.$sampleID)
tcr_list <- tcr_vdj2 %>% split(f=list(tcr_vdj2$sampleID))
length(names(tcr_list))
head(names(tcr_list))
tcr_list[[1]][1:4, 1:4]

3.4、 按样本保存拆分文件

# 2-4 Save small files
saveFiles <- function(inputObject, inputName){
  dir.create(inputName)
  write.csv(inputObject, file=file.path(inputName, 'filtered_contig_annotations.csv'), row.names=F, quote=F)
}

setwd('F:/R_Language/R_Practice/R_Packages/testPackage/data/data_cleaning_from_literature/tcr2')
purrr::pmap(list(inputObject=tcr_list, inputName=names(tcr_list)), saveFiles)

4、以样本为单位拆分BCR数据

4.1、 读取BCR文件

dataset_path <- 'F:/R_Language/R_Practice/R_Packages/testPackage/data/'
bcr_vdj <- read.delim(file.path(dataset_path, 'GSE158055_covid19_bcr_vdjnt_pclone.tsv'), check.names=T)
dim(bcr_vdj)
table(bcr_vdj$sampleID)
bcr_vdj[1:2, ]
bcr_vdj_sub <- bcr_vdj[1:6, ]
head(bcr_vdj_sub, 2)

4.2 、转换文件格式

从每个细胞一行信息转为两行

bcr_vdj2 <- bcr_vdj %>% tidyr::unite(TCRA, BCRH_cgene, BCRH_vgene, BCRH_dgene, BCRH_jgene, BCRH_cdr3aa, BCRH_cdr3nt) %>%
  tidyr::unite(TCRB, BCRL.K_cgene, BCRL.K_vgene, BCRL.K_dgene, BCRL.K_jgene, BCRL.K_cdr3aa, BCRL.K_cdr3nt) %>%
  tidyr::unite(TCRAB, TCRA, TCRB, sep='.') %>%
  tidyr::separate_rows(TCRAB, sep='\\.', convert=T) %>%
  dplyr::mutate(barcode = cellBarcode) %>%
  dplyr::select(barcode, sampleID, TCRAB) %>%
  dplyr::mutate(chain = rep(c('IGH', 'IGK/L'), times=n()/2)) %>%
  dplyr::mutate(is_cell='TRUE') %>%
  dplyr::mutate(high_confidence='TRUE') %>%
  dplyr::mutate(full_length='TRUE') %>%
  dplyr::mutate(productive='TRUE') %>%
  tidyr::separate(TCRAB, c("c_gene", "v_gene",'d_gene','j_gene', 'cdr3', 'cdr3_nt'), sep='\\_') %>%
  # tibble::add_column(is_cell='TRUE', .after='barcode') %>%
  dplyr::select(barcode, sampleID, is_cell, high_confidence, chain, v_gene, d_gene, j_gene, c_gene, full_length, productive, cdr3, cdr3_nt)

4.3 、按样本拆分文件

bcr_list <- bcr_vdj2 %>% split(f=list(bcr_vdj2$sampleID))
length(names(bcr_list))
head(names(bcr_list))
bcr_list[[1]][1:4, 1:4]

4.4 、按样本保存拆分文件

saveFiles <- function(inputObject, inputName){
  dir.create(inputName)
  write.csv(inputObject, file=file.path(inputName, 'filtered_contig_annotations.csv'), row.names=F, quote=F)
}

setwd('F:/R_Language/R_Practice/R_Packages/testPackage/data/data_cleaning_from_literature/bcr2')
purrr::pmap(list(inputObject=bcr_list, inputName=names(bcr_list)), saveFiles)

---

参考资料

[1]

This file is so big that R language can’t help it awk split a file into multiple files: https:///2021/08/20210831071133930b.html