盡管10Xgenomics已經(jīng)占據(jù)了單細(xì)胞測序的主流市場，Cellranger工具的開發(fā)也一站式解決了單細(xì)胞數(shù)據(jù)上游分析的難題。但是Smart-seq2仍然具有獨(dú)特的優(yōu)勢，并且在GEO數(shù)據(jù)庫中也沉淀了很多優(yōu)秀的Smart-seq2的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集，有時作者并不直接提供表達(dá)矩陣，因此，從上游開始獲取表達(dá)矩陣有時候成為唯一的方法。
但是相比Cellranger，處理Smart-seq2的工具卻十分有限，且網(wǎng)上的教程并不詳盡，容易引向歧途。
準(zhǔn)備
- 下載相關(guān)軟件
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
- 下載參考轉(zhuǎn)錄組
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
不同種屬的相關(guān)文件已經(jīng)有人整理并發(fā)布Releases · pachterlab/kallisto-transcriptome-indices
序列比對
Kallisto序列比對
Kallisto比對也稱為pseudo-alignment（假比對）。即它的比對方法和STAR不同，不是將完整reads比對到參考轉(zhuǎn)錄本上，而是將kmer（特定長度的序列）比對到參考轉(zhuǎn)錄本上。
1.SRR轉(zhuǎn)為fastq
ls SRR/SRR* | while read id;do \
(fastq-dump --gzip --split-3 -A `basename $id` -O rawdata $id &);done
2. 構(gòu)建參考基因組索引
mkdir indices/Kallisto
kallisto index -i indices/Kallisto/transcripts.idx ~/ref/transcripts.fa
只需要提供轉(zhuǎn)錄本的fasta格式的序列即可。-k參數(shù)指定kmer的長度，-i參數(shù)指定輸出的索引的名字，注意kallisto建立的索引為一個文件。或者可以通過kallisto-transcriptome-indices進(jìn)行下載。
2. Bus格式轉(zhuǎn)化
隨著kallisto版本的更新，開始支持越來越多的單細(xì)胞測序技術(shù)如10X genomic(V1-V3)、CELSEQ(V1-V3)、DROPSEQ、SMARTSEQ。你可以用如下的命令查看該軟件支持的所有測序技術(shù)。需要注意的是，在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
執(zhí)行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：表達(dá)量的矩陣等價類文件；
- output.bus：原始的BUS(Barcode,UMI,Set format）格式比對結(jié)果文件；
- run_info.json：比對結(jié)果的概要信息；
- transcripts.txt：定量生成的轉(zhuǎn)錄本文件；
目錄結(jié)構(gòu)如下：
SRR9169420_out
├── matrix.ec
├── output.bus
├── run_info.json
└── transcripts.txt
3. scRNAseq的定量
bustools分三個步驟對bus文件進(jìn)行定量
- 校正（correct）
- 排序（sort）
- 定量（count）
1.校正
在店小二的筆記中，可以通過10X的白名單對barcode進(jìn)行校正。但筆者沒有找到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
定量完成后會生成一系列的結(jié)果文件，最重要的表達(dá)量文件是cells_x_genes.barcodes.txt、cells_x_genes.genes.txt、cells_x_genes.mtx(類似cellranger的結(jié)果)。
- 批量運(yùn)行
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
參考內(nèi)容
- https://links.jianshu.com/go?to=https%3A%2F%2Fgithub.com%2Fpachterlab%2Fkallisto-transcriptome-indices%2Freleases
- smart-seq2 實(shí)踐 - 簡書 (jianshu.com)
- Releases · pachterlab/kallisto-transcriptome-indices (github.com)
- Manual (bustools.github.io)
- kallisto比對參考轉(zhuǎn)錄組 - 云+社區(qū) - 騰訊云 (tencent.com)
- scRNA-seq數(shù)據(jù)預(yù)處理 | 生信拾光 (renyx.top)
