尽管10Xgenomics已经占据了单细胞测序的主流市场，Cellranger工具的开发也一站式解决了单细胞数据上游分析的难题。但是Smart-seq2仍然具有独特的优势，并且在GEO数据库中也沉淀了很多优秀的Smart-seq2的单细胞转录组数据集，有时作者并不直接提供表达矩阵，因此，从上游开始获取表达矩阵有时候成为唯一的方法。
但是相比Cellranger，处理Smart-seq2的工具却十分有限，且网上的教程并不详尽，容易引向歧途。

准备

• 下载相关软件

conda search kallisto
conda install kallisto=0.48.0
conda install bustools
conda install sra-tools
conda install axel

• 下载测序文件

prefetch --option-file SRR_Acc_List.txt

• 下载参考基因组

cd ~/ref
axel https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/GRCh38.p13.genome.fa.gz

• 下载参考转录组

axel https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.transcripts.fa.gz

• 下载基因注释信息

axel https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz

不同种属的相关文件已经有人整理并发布Releases · pachterlab/kallisto-transcriptome-indices

序列比对

Kallisto序列比对
Kallisto比对也称为pseudo-alignment（假比对）。即它的比对方法和STAR不同，不是将完整reads比对到参考转录本上，而是将kmer（特定长度的序列）比对到参考转录本上。

1.SRR转为fastq

ls SRR/SRR* | while read id;do \
    (fastq-dump --gzip --split-3 -A `basename $id` -O rawdata $id &);done

2. 构建参考基因组索引

mkdir indices/Kallisto
kallisto index -i indices/Kallisto/transcripts.idx ~/ref/transcripts.fa

只需要提供转录本的fasta格式的序列即可。-k参数指定kmer的长度，-i参数指定输出的索引的名字，注意kallisto建立的索引为一个文件。或者可以通过kallisto-transcriptome-indices进行下载。

2. Bus格式转化

随着kallisto版本的更新，开始支持越来越多的单细胞测序技术如10X genomic(V1-V3)、CELSEQ(V1-V3)、DROPSEQ、SMARTSEQ。你可以用如下的命令查看该软件支持的所有测序技术。需要注意的是，在kallisto-0.44.0版本以前，并不支持此功能。本文章所用的为0.48.0版本。

kallisto bus -l
List of supported single-cell technologies

short name       description
----------       -----------
10xv1            10x version 1 chemistry
10xv2            10x version 2 chemistry
10xv3            10x version 3 chemistry
Bulk             Bulk RNA-seq or Smart-seq2 (multiplexed)
BDWTA            BD Rhapsody WTA
CELSeq           CEL-Seq
CELSeq2          CEL-Seq version 2
DropSeq          DropSeq
inDropsv1        inDrops version 1 chemistry
inDropsv2        inDrops version 2 chemistry
inDropsv3        inDrops version 3 chemistry
SCRBSeq          SCRB-Seq
SmartSeq3        Smart-seq3
SPLiT-seq        SPLiT-seq
SureCell         SureCell for ddSEQ
Visium           10x Visium Spatial Transcriptomics

执行kallisto bus命令，生成bus文件。

mkdir busdata
fq_dir=./rawdata
bus_dir=./busdata

cat  SRR_Acc_List.txt|while read id;do
kallisto bus  -x Bulk -i indices/Kallisto/transcripts.idx \
 -o $bus_dir/$id $fq_dir/$id.sra_1.fastq.gz $fq_dir/$id.sra_2.fastq.gz; done

在各自的目录下生成四个文件：

• matrix.ec：表达量的矩阵等价类文件；
• output.bus：原始的BUS(Barcode,UMI,Set format）格式比对结果文件；
• run_info.json：比对结果的概要信息；
• transcripts.txt：定量生成的转录本文件；
  目录结构如下：

SRR9169420_out
├── matrix.ec
├── output.bus
├── run_info.json
└── transcripts.txt

3. scRNAseq的定量

bustools分三个步骤对bus文件进行定量

1. 校正（correct）
2. 排序（sort）
3. 定量（count）

1.校正
在店小二的笔记中，可以通过10X的白名单对barcode进行校正。但笔者没有找到Smart-seq2的类似文件，故第一步暂时跳过。

bustools correct -w 10xv2_whitelist.txt -p SRR9169420_out/output.bus

2. 排序

bustools sort  -t 10 -p -o output.sort.bus out.bus

3. 定量

bustools count -o genecount -g ~/ref/homo_sapiens/transcripts_to_genes.txt" \
-e matrix.ec -t transcripts.txt --genecounts output.sort.bus 

定量完成后会生成一系列的结果文件，最重要的表达量文件是cells_x_genes.barcodes.txt、cells_x_genes.genes.txt、cells_x_genes.mtx(类似cellranger的结果)。

4. 批量运行

cat SRR_Acc_List.txt | while read id; do
  cd bus_data/$id
  bustools sort  -t 10 -T -p -o output.sort.bus output.bus
  bustools count -o genecount -g ~/ref/homo_sapiens/transcripts_to_genes.txt  -e matrix.ec -t transcripts.txt --genecounts output.sort.bus 
  cd ../..
done

参考内容

1. https://links.jianshu.com/go?to=https%3A%2F%2Fgithub.com%2Fpachterlab%2Fkallisto-transcriptome-indices%2Freleases
2. smart-seq2 实践 - 简书 (jianshu.com)
3. Releases · pachterlab/kallisto-transcriptome-indices (github.com)
4. Manual (bustools.github.io)
5. kallisto比对参考转录组 - 云+社区 - 腾讯云 (tencent.com)
6. scRNA-seq数据预处理 | 生信拾光 (renyx.top)