一文打通单细胞上游：从软件部署到上游分析

虽然我分享了：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，但是我发现我一直是从ENA数据库使用aspera进行高速下载，但是它有两个弊端，首先是ena数据库的不稳定，其次是部分数据集的单细胞样品在ena上面并不是R1,R2,I1的3个fastq文件形式。

所以补上这个“生信随笔”的笔记，如果你需要从SRA数据库下载sra文件 ：

下面是公众号“生信随笔”的分享

导语

本文介绍在Linux环境下10X单细胞上游分析环境部署，主要从软件安装和上游分析示例讲起。上游分析示例部分包括两个内容：

• 常规的SRA文件下载，转fastq格式，然后走cellranger流程；

• 此外，有些时候作者提供的是bam文件，我们可以先把bam转为fastq，再走cellranger流程。

一. 部署10x上游分析环境

(1) conda创建分析环境及常用软件安装

#conda install mamba
conda create -n 10x
conda activate 10x
mamba install -y -c bioconda aspera-cli bwa samtools bedtools sambamba sra-tools bowtie2 samblaster fasterq-dump

(2) CellRanger及参考基因组下载

2021年12月的最新版本是cellranger软件6.2.1，hg38及mm10参考基因组。

这里cellranger的下载链接会失效，需要去官网获取：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest

注：服务器下载慢的话可以在windows环境下载，然后上传压缩文件至服务器

##1.1 cellranger软件6.2.1
wget -O cellranger-6.1.2.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-6.1.2.tar.gz?Expires=1636331862&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci02LjEuMi50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MzYzMzE4NjJ9fX1dfQ__&Signature=LVCM9NLNJQ4BgDFmIdLG7jYI~hoT4QB4nNIT9U6krn~gVRHQwLP25YDI~msH6cq6fpd8OC4MTNlLNUrqVFTlBK0~qq9Cq~MA-ofUBxYt74yECf4vgGd5uRTUKsiKy2af~6~kicBaj5ltCoIjmguKVQPof4M2E81IdkMUHsVBp6NoCfQCI68sRtIMyYTx7~fpvY1MZv6PipPUXMctB85usluXc~vMQT6f-W4q-gYKG31y8wqS~4Idh4l1oQpMU4T-I16E-EO04~aLQP9lJSMuswZsUxmOP8~wkR93ULmW16UmCH6YYAXAsVZipc6xe1k12yPZ--ioAT9x4JwZW8EIUQ__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"

md5sum cellranger-6.1.2.tar
#310d4453acacf0eec52e76aded14024c cellranger-6.1.2.tar
tar -xzvf cellranger-6.1.2.tar

##1.2 mm10
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz

md5sum refdata-gex-mm10-2020-A.tar.gz
#886eeddde8731ffb58552d0bb81f533d refdata-gex-mm10-2020-A.tar.gz
tar -xzvf refdata-gex-mm10-2020-A.tar.gz

##1.3 hg38
wget https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-GRCh38-2020-A.tar.gz

md5sum refdata-gex-GRCh38-2020-A.tar.gz
#dfd654de39bff23917471e7fcc7a00cd refdata-gex-GRCh38-2020-A.tar.gz
tar -xzvf refdata-gex-GRCh38-2020-A.tar.gz

##1.4 环境变量配置：PATH的路径依据自己实际的情况而定
ln -s ~/cellrange_soft/cellranger-6.1.2/bin/cellranger ~/anaconda3/envs/10x/bin/

cellranger

二. 常规的单细胞SRR文件格式上游分析

创建分析所用的文件夹备用：

mkdir 1.sra 2.raw_fastq 3.cellranger

(1) 数据下载

示例数据GSE117988：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117988

Metadata是表型数据，Accession List提供了SRA ID：

cd 1.sra

## 写入需要下载的文件名
cat >download_file
SRR7722937
SRR7722938
SRR7722939
SRR7722940
SRR7722941
SRR7722942

# 批量下载
cat download_file |while read id;do (prefetch $id &);done  # 后台下载

(2) fasterq-dump SRA转fastq

常规的SRA转fastq文件，用的是fastq-dump软件，速度非常慢，4-5个小时才能处理完一个样本。

这里用新办法fasterq-dump，2分钟完成一个样本：

cd 2.raw_fastq
ln -s ../1.sra/SRR* ./

### 方法1：多个文件批量做
cat >fastq.sh
ls SRR* | while read id;do ( nohup fasterq-dump -O ./ --split-files -e 6 ./$id  --include-technical & );done
## 运行脚本
bash fastq.sh

### 方法2：for循环一个一个做
for i in `ls SRR*`
do
i=$i
echo "fasterq-dump -O ./ --split-files -e 40 ./$i --include-technical"
done >fastq.sh
## 运行脚本
bash fastq.sh

### 大概耗时2分钟完成
ls -lh

如果觉得占空间，可压缩一下：

### 压缩一下
ls SRR*fastq | while read id; do gzip $id; done

拿到了这些每个样品的多个fastq文件，后面的流程就跟我分享的：10X单细胞转录组原始测序数据的Cell Ranger流程（仅需800元）以及一个10x单细胞转录组项目从fastq到细胞亚群，是一回事啦。

(3) Cellranger count流程

对10X的fq文件运行CellRangerd的counts流程，先做一个测试：

首先需要对fastq.gz文件改名字，SampleName_S1_L001_R1_001.fastq.gz：

mv SRR7722937_1.fastq.gz SRR7722937_S1_L001_I1_001.fastq.gz
mv SRR7722937_2.fastq.gz SRR7722937_S1_L001_R1_001.fastq.gz
mv SRR7722937_3.fastq.gz SRR7722937_S1_L001_R2_001.fastq.gZ

再运行cellranger：

cd 3.cellranger

ref=/home/data/vip10t17/software_install/10x_refernce/refdata-gex-GRCh38-2020-A
id=SRR7722937
cellranger count --id=$id \
--transcriptome=$ref \
--fastqs=/home/data/vip10t17/GEO_data/10x_test/2.raw_fastq \
--sample=$id \
--nosecondary \
--localmem=30

会一个的话，就可以写shell脚本批量完成：

##第一步 改fastq文件名字：
cd 2.raw_fastq

cat ../1.sar/download_file | while read i ;do (mv ${i}_1*.gz 
${i}_S1_L001_I1_001.fastq.gz;mv ${i}_2*.gz ${i}_S1_L001_R1_001.fastq.gz;mv 
${i}_3*.gz ${i}_S1_L001_R2_001.fastq.gz);done

##第二步 批量运行cellranger：
ref=/home/data/vip10t17/software_install/10x_refernce/refdata-gex-GRCh38-2020-A
ls *.fastq.gz | cut -d "_" -f 1 | uniq | while read id;
do
nohup cellranger count --id=$id \
--transcriptome=$ref \
--fastqs=/home/data/vip10t17/GEO_data/10x_test/2.raw_fastq \
--sample=$id \
--nosecondary \
--localcores=10 \ #设置核心数
--localmem=30 &
done

##第三步可选 打包压缩文件方便传输
tar -zcvf Output.tar.gz SRR7722937 SRR7722938 SRR7722939 SRR7722940 SRR7722941 SRR7722942

最主要的几个参数：

• --id 指定输出文件夹的名字

• --transcriptome 指定参考基因组的路径

• --sample 指定需要处理的fastq文件的前缀

• --expect-cell 指定预期的细胞数目，默认参数是3000个

• --localcores 指定计算的核心数

• --mempercore 指定内存大小 GB

• --nosecondary 不需要进行降维聚类（因为后期会用R可视化）

• --chemistry，默认是自动识别chemistry，但是有些时候识别不出chemistry的时候，需要加入参数特别标明

使用cellranger count --help可查看更多参数，Cellranger count结果解读见文末。

三. 如果提供的是BAM数据

注意，bam转fastq，再走cellranger的流程非常耗费计算机资源和时间，八个样本预计会占2T的硬盘空间。练习的话，笔者建议下载1-2个样本即可。

(1) 下载BAM数据：

示例数据来自：https://www.ebi.ac.uk/ena/browser/view/PRJNA727404?show=reads

cd ~/Projects_gao/newhbvhc
##3.1构建下载list
cat >download_file
/vol1/run/SRR144/SRR14424777/740_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424778/725_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424779/713_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424780/119_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424781/114_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424782/106_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424783/104_possorted_genome_bam.bam
/vol1/run/SRR144/SRR14424784/095_possorted_genome_bam.bam

##3.2批量下载
nohup ascp -v -QT -l 300m -P33001 -k1 -i /home/data/ssy40/anaconda3/envs/10x/etc/asperaweb_id_dsa.openssh --mode recv --host fasp.sra.ebi.ac.uk --user era-fasp --file-list download_file ./ & >download.log

ls -lh | cut -d" " -f 5-
 56G 5月   6  2021 095_possorted_genome_bam.bam
 23G 5月   9  2021 104_possorted_genome_bam.bam
 22G 5月  11  2021 106_possorted_genome_bam.bam
 51G 5月  22  2021 114_possorted_genome_bam.bam
 45G 5月   6  2021 119_possorted_genome_bam.bam
 22G 5月   8  2021 713_possorted_genome_bam.bam
 27G 5月   8  2021 725_possorted_genome_bam.bam
 53G 12月 11  2021 740_possorted_genome_bam.bam
449K 12月 10 13:45 wget-log

(2) bam转fastq

参考官网流程：https://support.10xgenomics.com/docs/bamtofastq?src=pr&lss=none&cnm=&cid=NULL

cellranger产生的bam文件里是带有barcode与UMI的，储存在tag标签里：

samtools view 104_possorted_genome_bam.bam | less -SN
samtools view 104_possorted_genome_bam.bam | head -3 | tr "\t" "\n" | cat -n

• CB、CR、CY表示barcode，一般是16个碱基；

• UB、UR、UY表示UMI，一般是10个碱基。

• R一般代表原始测序数据，Y代表质量分数，而B代表校正后的R，可能对应碱基质量分数太低等因素。一般来说R与B都是相同的。

cellranger bamtofastq:

## 3.3构建一个name list文件
cat >name.list
740
725
713
119
114
106
104
095

##3.4创建文件夹
mkdir fastq_file

##3.5准备shell脚本
cat > bamtofastq.sh
cat name.list |while read id
do
cellranger bamtofastq --nthreads 30 --traceback ${id}_possorted_genome_bam.bam ./fastq_file/${id}
done
#运行脚本
nohup bash bamtofastq.sh & 

#必要时可批量kill任务
#ps -ef | grep bamtofastq | awk '{print $2}' | while read id;do kill $id;done  #批量Kill

(3) cellranger count

批量完成：

### 看一下文件夹地址
find ~/Projects_gao/newhbvhc/fastq_file/*/*count*/ | grep [1-1000] | grep -v XX/bam 

##3.7 输入剩余未完成的文件list
cat >other_file.list
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/095/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/104/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/106/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/114/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/119/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/713/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/725/
/home/data/ssy40/Projects_gao/newhbvhc/fastq_file/740/

##3.8 批量shell代码
cat other_file.list |while read id
do
ref=/home/data/ssy40/cellrange_soft/10x_refernce/refdata-gex-GRCh38-2020-A
sample_name=${id:0-36:3}_out_file
echo "cellranger count --id=$sample_name \
--transcriptome=$ref \
--fastqs=$id \
--sample=bamtofastq \
--nosecondary \
--localmem=20 \
--localcores=30"
done>other_file.sh

#运行
nohup bash other_file.sh>log_106.114.119.log 2>&1 &

四. Cellranger count结果解读

例如其中一个样本产生的结果：

最重要的结果在outs/文件夹下：

结果解读（参考生信技能树Jimmy的帖子）：

• web_summary.html：必看，官方说明 summary HTML file ，包括许多QC指标，预估细胞数，比对率等；

• metrics_summary.csv：CSV格式数据摘要，可以不看；

• possorted_genome_bam.bam：比对文件，用于可视化比对的reads和重新创建FASTQ文件，可以不看；

• possorted_genome_bam.bam.bai：索引文件；

• filtered_gene_bc_matrices：是重要的一个目录，下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件，是经过Cell Ranger过滤后构建矩阵所需要的所有文件；

• filtered_feature_bc_matrix.h5：过滤掉的barcode信息HDF5 format，可以不看；

• raw_feature_bc_matrix：原始barcode信息，未过滤的可以用于构建矩阵的文件，可以不看；

• raw_feature_bc_matrix.h5：原始barcode信息HDF5 format，可以不看；

• analysis：数据分析目录，下面又包含聚类clustering（有graph-based & k-means）、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne，因为我们自己会走Seurat流程，所以不用看；

• molecule_info.h5：可用于整合多样本，使用cellranger aggr函数；

• cloupe.cloupe：官方可视化工具Loupe Cell Browser 输入文件，无代码分析的情况下使用，会代码的同学通常用不到。