尽管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文件进行定量
- 校正（correct）
- 排序（sort）
- 定量（count）
1.校正
在店小二的笔记中，可以通过10X的白名单对barcode进行校正。但笔者没有找到Smart-seq2的类似文件，故第一步暂时跳过。
bustools correct -w 10xv2_whitelist.txt -p SRR9169420_out/output.bus
- 排序
bustools sort -t 10 -p -o output.sort.bus out.bus
- 定量
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的结果)。
- 批量运行
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
参考内容
- https://links.jianshu.com/go?to=https%3A%2F%2Fgithub.com%2Fpachterlab%2Fkallisto-transcriptome-indices%2Freleases
- smart-seq2 实践 - 简书 (jianshu.com)
- Releases · pachterlab/kallisto-transcriptome-indices (github.com)
- Manual (bustools.github.io)
- kallisto比对参考转录组 - 云+社区 - 腾讯云 (tencent.com)
- scRNA-seq数据预处理 | 生信拾光 (renyx.top)
网友评论