单细胞转录组学习流程&总结——（一）背景原理+软件安装+上游分析

学习资料

• 单细胞转录组上游分析-上游分析流程及知识点整理
• 单细胞转录组(Single cell RNA)概述
• 单细胞实战
• 简书上的一个博主-刘小泽，有很多关于单细胞的学习资料：https://www.jianshu.com/u/d7b77c171c15
• 单细胞的教程，很详细，有配套的代码： Analysis of single cell RNA-seq data
  
• 相关的综述 Identifying cell populations with scRNASeq
• 这是10×Genomics官方的文档：https://www.10xgenomics.com/support/user-guides/single-cell-gene-expression?menu[productNames]=Single%20Cell%20Gene%20Expression
• 单细胞转录组流程一：安装运行cellranger及结果解读
• 介绍软件以及用法：10X单细胞测序分析软件:Cell ranger
• 终于弄懂了文件的含义：单细胞测序上游分析-从原始数据到cellranger定量
• 单细胞转录组之cellranger探索——1.安装
• 单细胞基本分析流程学习总结

上游分析

软件安装

1. 进入官网
   What is Cell Ranger
   填基本信息
   Download Center
   
   下载的最新版cellranger-8.0.1

wget -O cellranger-8.0.1.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-8.0.1.tar.gz?Expires=1719930166&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA&Signature=Eu~0hCqtqPZW518BWLlZlpwtWhrgClaIgnhn9-wEZuqoqnFCZC374rGxOImxocnjdnyvuVhYk-fGgPW2bmQVumBVWsx8Rgj7zj7YqV2YDddWryND3gw2ft0M3fesb3Hs7jYwaluF-qZvmnckKhjnRK2Cff9~zJbcBr73Ef6B7492Kmmmpfjoshirw4BtSqt6BV-W273vYHGIz9CEe6EPoTn343cVXyxW~TlXT118N1y3bh2RUpFuR4MtrcB1UBlEz9FEn6B5YznnncAr4mub8YREjtUfkbk-NwzhodmxSnGdbTYptyFDadXM2PYekU4H~HI8TzH0TW5I81020rIRow__"

tar -xzvf cellranger-8.0.1.tar.gz

虽然安装完毕，但是现在要运行cellranger的话需要在终端输入cellranger的整个路径。
将cellranger的路径加入$PATH，这样只需要在终端输入cellranger就可以执行命令。
将cellranger加入 $PATH有两种方式，分为临时和永久。

临时加载环境变量：
将路径加入$PATH：

(base) [animal_geneticB@master-cli1-x86-agent1 cellranger-8.0.1]$ export PATH=/workspace/home/pig/animal_geneticB/software/cellranger-8.0.1
(base) [animal_geneticB@master-cli1-x86-agent1 cellranger-8.0.1]$ echo $PATH
/workspace/home/pig/animal_geneticB/software/cellranger-8.0.1

注：
现在终端输入cellranger就可以执行命令了。

操作仅临时有用，如果电脑重启后就实效了，得重新加入路径。

永久加载环境变量：
自动执行这个命令，需要将export PATH=自己的路径:$PATH这个命令加入.bashrc文件。

vim编辑器打开.bashrc文件：

vim .bashrc

在.bashrc文件最后加入export PATH=自己的路径:$PATH
source .bashrc立即生效：

source .bashrc

但是vim .bashrc就出现了错误

(base) [animal_geneticB@master-cli1-x86-agent1 ~]$ vim ~/.bashrc
-bash: vim: command not found

终于知道是怎么回事了，是由于我已经临时加载环境变量，导致那些命令都不能用了，所以后面重新进入Xshell，并运行以下命令：

(base) [animal_geneticB@master-cli1-x86-agent1 cellranger-8.0.1]$ vim ~/.bashrc
export PATH=/workspace/home/pig/animal_geneticB/software/cellranger-8.0.1:$PATH

先退出Conda环境，再运行以下命令：

(base) [animal_geneticB@master-cli1-x86-agent1 ~]$ conda deactivate
[animal_geneticB@master-cli1-x86-agent1 ~]$ source .bashrc
(base) [animal_geneticB@master-cli1-x86-agent1 ~]$ echo $PATH
/workspace/home/pig/animal_geneticB/software/cellranger-8.0.1:PATH/edgeturbo-client:/workspace/home/pig/animal_geneticB/.aspera/connect/bin:/workspace/home/pig/animal_geneticB/.aspera/connect/bin:/workspace/home/pig/animal_geneticB/anaconda3/bin:/workspace/home/pig/animal_geneticB/anaconda3/condabin:/workspace/home/pig/animal_geneticB/.local/bin:/workspace/home/pig/animal_geneticB/bin:/opt/batch/cli/bin:/opt/batch/cli/bin:/opt/batch/cli/bin:/usr/share/Modules/bin:/usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin:/home/usersname/.aspera/connect/bin/

可以看到，cellranger-8.0.1已经被添加到环境变量中了。

测试安装是否成功：

(base) [animal_geneticB@master-cli1-x86-agent1 sn_RNAseq]$ cellranger testrun --id=tiny

参考：

• cellranger安装以及使用

文件重命名

FASTQs file naming convention
To serve as inputs for Cell Ranger ARC, FASTQ files should conform to the naming conventions of bcl2fastq and mkfastq described below.
GEX FASTQs
[Sample Name]_S1_L00[Lane Number]_[Read Type]_001.fastq.gz
Where Read Type is one of:
I1: Dual index i7 read (optional)
I2: Dual index i5 read (optional)
R1: Read 1
R2: Read 2

重命名，脚本如下，命名为rename.sh

for file in *f1.fastq.gz; do
    prefix="${file%_f1.fastq.gz}"
    mv "$file" "${prefix}_S1_L001_R1_001.fastq.gz"
done

for file in *r2.fastq.gz; do
    prefix="${file%_r2.fastq.gz}"
    mv "$file" "${prefix}_S1_L001_R2_001.fastq.gz"
done

运行脚本：

(base) [animal_geneticB@master-cli1-x86-agent1 sn_RNAseq]$ chmod 777 rename.sh 
(base) [animal_geneticB@master-cli1-x86-agent1 sn_RNAseq]$ bash rename.sh 
(base) [animal_geneticB@master-cli1-x86-agent1 sn_RNAseq]$ ls -lh
total 72G
-rw------- 1 animal_geneticB pig 5.8G Jun 28 13:00 CRR784294_S1_L001_R1_001.fastq.gz
-rw------- 1 animal_geneticB pig  13G Jun 28 13:20 CRR784294_S1_L001_R2_001.fastq.gz
-rw------- 1 animal_geneticB pig 5.6G Jun 28 12:29 CRR784295_S1_L001_R1_001.fastq.gz
-rw------- 1 animal_geneticB pig  12G Jun 28 12:50 CRR784295_S1_L001_R2_001.fastq.gz
-rw------- 1 animal_geneticB pig 5.5G Jun 28 11:51 CRR784296_S1_L001_R1_001.fastq.gz
-rw------- 1 animal_geneticB pig  13G Jun 28 12:18 CRR784296_S1_L001_R2_001.fastq.gz
-rw------- 1 animal_geneticB pig 5.3G Jun 28 11:14 CRR784297_S1_L001_R1_001.fastq.gz
-rw------- 1 animal_geneticB pig  13G Jun 28 11:37 CRR784297_S1_L001_R2_001.fastq.gz

由于这里作者上传的只有两边Read的fastq，不包含Index，因此只批量命名了两个Read文件名称，如果项目包含了index，可根据需求结合cellranger官方建议进行改写。

参考：

• Specifying Input FASTQ Files for cellranger-arc count
• 单细胞测序上游分析-从原始数据到cellranger定量

构建索引

直接用之前跑RNA_seq的基因组文件和注释文件
cellranger 操作笔记-2：构建绵羊单细胞转录组参考基因组

• 过滤注释文件

cellranger mkgtf Sus_scrofa.Sscrofa11.1.112.gtf Sus_scrofa.Sscrofa11.1.112.filtered.gtf --attribute=gene_biotype:protein_coding

• 构建参考基因组index

#=================================================
#运行脚本
#=================================================

cellranger mkref \
--genome=Sus_scrofa_index \
--fasta="/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/reference/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa" \
--genes="/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/reference/Sus_scrofa.Sscrofa11.1.112.filtered.gtf"

#有几个需要注意的地方：
#genome后只能跟输出文件夹，不能有路径；
#genome的文件夹不能提前创建

(base) [animal_geneticB@master-cli1-x86-agent1 reference]$ dsub -s run_index.sh 

构建索引大概花了40min，检查文件是否齐全：

(base) [animal_geneticB@master-cli1-x86-agent1 reference]$ tree Sus_scrofa_index
Sus_scrofa_index
├── fasta
│   ├── genome.fa
│   └── genome.fa.fai
├── genes
│   └── genes.gtf.gz
├── reference.json
└── star
    ├── chrLength.txt
    ├── chrNameLength.txt
    ├── chrName.txt
    ├── chrStart.txt
    ├── exonGeTrInfo.tab
    ├── exonInfo.tab
    ├── geneInfo.tab
    ├── Genome
    ├── genomeParameters.txt
    ├── SA
    ├── SAindex
    ├── sjdbInfo.txt
    ├── sjdbList.fromGTF.out.tab
    ├── sjdbList.out.tab
    └── transcriptInfo.tab

4 directories, 19 files

Count

Cell Ranger count使用手册
我是有四个样本，所以建了四个目录：

(base) [animal_geneticB@master-cli1-x86-agent1 raw_data]$ cat err_4003.log 
error: the following required arguments were not provided:
  --create-bam <true|false>

Usage: cellranger count --id <ID> --transcriptome <PATH> --create-bam <true|false> --fastqs <PATH> --sample <PREFIX> --nosecondary

之前出现了报错，是因为没有加--create-bam这个参数，最终脚本如下：

#===========================================================
#运行脚本
#===========================================================
# Define the paths to the FASTQ files
cellranger count \
--id=sn_RNA_count \
--transcriptome=/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/reference/Sus_scrofa_index \
--fastqs=/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/raw_data/CRR784294,/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/raw_data/CRR784295,/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/raw_data/CRR784296,/workspace/home/pig/animal_geneticB/Analysis/py/data/sn_RNAseq/raw_data/CRR784297 \
--sample=CRR784294,CRR784295,CRR784296,CRR784297 \
--nosecondary \
--create-bam=false  # 或者使用 --create-bam=true 如果你需要 BAM 文件 
    # id指定输出文件存放目录名
    # transcriptome指定与CellRanger兼容的参考基因组
    # fastqs指定mkfastq或者自定义的测序文件,如果需要对多个fastq进行合并比对的话，可以输入多个存放路径并用逗号隔开，相对应的在--sample中，也输入多个需要合并的sample名称并用逗号隔开，这样就不用我们额外手动合并了
    # sample要和fastq文件的前缀中的sample保持一致，作为软件识别的标志
    # expect-cells指定复现的细胞数量，这个要和实验设计结合起来
    # nosecondary 只获得表达矩阵，不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做)

脚本如下，命名为run_count.sh。