Hemberg-lab单细胞转录组数据分析（五）-CFANZ编程社区

使用STAR比对read

现在我们已经对测序原始数据进行了质控，获得了高质量的Clean data，下一步就是把它们比对到参考基因组。如果我们想定量基因表达或鉴定样本之间差异表达的基因，则通常需要某种形式的比对。

用于短序列比对的工具已经开发了很多（转录组分析工具哪家强？），但今天我们主要涉及两个。第一个工具是STAR。对于测序数据中的每条reads，STAR尝试找到能与参考基因组中一个或多个位置匹配的最长可能序列。例如，在下图中，有一个短序列（蓝色），它跨越两个外显子和一个选择性剪接点（紫色）。STAR能够发现短序列的第一部分与第一外显子的序列匹配，而第二部分与第二外显子中的序列匹配。因为STAR能够用这种方式识别剪接事件，所以它被称为splice aware的比对工具 (一般转录组分析的比对工具都需要有这个功能)。

通常STAR把短序列比对到参考基因组时允许检测新的剪接事件或染色体重排事件。然而，STAR的一个问题是它需要大量的内存，尤其是参考基因组很大（例如老鼠和人类，大约需要30 G内存）的时候。为了加速今天的分析，我们将使用STAR把reads比对到只包含2000个转录本的参考转录组上。请注意，这不是常用或推荐的做法，只是考虑时间问题才这样做。我们通常建议比对到参考基因组 (但过程与此类似)。

执行STAR比对需要两个步骤。在第一步中，用户向STAR提供参考基因组序列（FASTA）和注释（GTF）(见NGS基础 - 参考基因组和基因注释文件)，来创建基因组索引 (生信宝典注：建索引的目的是为了提高比对速度)。第二步，STAR将用户的短序列数据比对到基因组索引。

首先创建索引。请记住，由于时间的原因，我们练习时是比对到转录组而不是基因组，所以我们只需要向STAR提供用于比对的转录本序列即可，不需要基因注释文件。可以从 Ensembl (https://www.ensembl.org/info/data/ftp/index.html)中得到许多模式生物的转录本序列。

任务1：用下面命令创建索引

mkdir indices
mkdir indices/STAR
STAR --runThreadN 4 --runMode genomeGenerate --genomeDir indices/STAR --genomeFastaFiles Share/2000_reference.transcripts.fa

任务2：STAR用到的每个选项表示什么意思？提示：用STAR手册来帮助你 https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf

任务3：如果我们要比对到基因组而不是转录组，命令将有何不同？

# 实际应用时比对到基因组
# 命令如下
mkdir -p star_GRCh38
# --runThreadN 2: 指定使用2个线程
# --sjdbOverhang 100: 默认
STAR --runMode genomeGenerate --runThreadN 2 --genomeDir star_GRCh38 \
     --genomeFastaFiles GRCh38.fa --sjdbGTFfile GRCh38.gtf

# 输出文件如下，注意看下输出文件的大小，有无空文件，基因数是否对。
#
ls -sh star_GRCh38
#
# 总用量 2.1G
# 4.0K chrLength.txt      368K exonInfo.tab          1.5G SAindex
# 4.0K chrNameLength.txt   24K geneInfo.tab          204K sjdbInfo.txt
# 4.0K chrName.txt         64M Genome                204K sjdbList.fromGTF.out.tab
# 4.0K chrStart.txt       4.0K genomeParameters.txt  204K sjdbList.out.tab
# 732K exonGeTrInfo.tab   516M SA                    224K transcriptInfo.tab
#
# STAR解析后的基因数
wc -l star_GRCh38/geneInfo.tab
#
# 原始GTF的基因数
grep -cP '\tgene\t' GRCh38.gtf

现在完成了索引创建，执行比对步骤。

任务4：尝试写出用STAR比对reads到构建好的索引的命令。使用STAR手册来帮助你。如果你觉得知道了答案，在下一节中核对它是否与解决方案匹配并执行比对。

任务5：尝试了解比对输出的结果文件的意义。

STAR比对命令

mkdir results
mkdir results/STAR

STAR --runThreadN 4 --genomeDir indices/STAR --readFilesIn Share/ERR522959_1.fastq Share/ERR522959_2.fastq --outFileNamePrefix results/STAR/

# 更复杂些的比对命令和参数解释
# 生信宝典
mkdir -p trt_N061011
# --runThreadN 4: 使用4个线程
# --readFilesIn: 输入文件，左端和右端
# --readFilesCommand zcat：gzip压缩，指定解压方式
# --genomeDir：基因组索引目录的位置
# -S: 指定输出文件

# 动物一般写 1000000，植物一般写5000
max_intron_size=100000

# --genomeLoad LoadAndKeep : 共享内存
# 其他参数自己对着star的帮助手册看
star_p=" --outFilterType BySJout --outSAMattributes NH HI AS NM MD \
       --outFilterMultimapNmax 20 --alignSJoverhangMin 8 --alignSJDBoverhangMin 1 \
       --alignIntronMin 20 --alignIntronMax ${max_intron_size} \
                --alignMatesGapMax ${max_intron_size} \
                --outFilterMatchNminOverLread 0.66 --outFilterScoreMinOverLread 0.66 \
                --winAnchorMultimapNmax 70 --seedSearchStartLmax 45 \
                --outSAMattrIHstart 0 --outSAMstrandField intronMotif \
                --genomeLoad LoadAndKeep --outReadsUnmapped Fastx \
                --outSAMtype BAM Unsorted --quantMode TranscriptomeSAM GeneCounts"

# STAR比对单个样品                
STAR --runMode alignReads --runThreadN 4 \
        --readFilesIn trt_N061011_1.fq.gz trt_N061011_2.fq.gz \
        --readFilesCommand zcat --genomeDir genome/star_GRCh38 \
        --outFileNamePrefix trt_N061011/trt_N061011\. ${star_p}

# STAR比对单个样品                
STAR --runMode alignReads --runThreadN 4 \
        --readFilesIn trt_N061011_1.fq.gz trt_N061011_2.fq.gz \
        --readFilesCommand zcat --genomeDir genome/star_GRCh38 \
        --outFileNamePrefix trt_N061011/trt_N061011\. ${star_p}

# Aug 03 11:44:27 ..... started STAR run
# Aug 03 11:44:27 ..... loading genome
# Aug 03 11:44:30 ..... started mapping
# Aug 03 11:44:48 ..... finished successfully

# 输出文件

ls -sh trt_N061011/*

# trt_N061011.Aligned.out.bam: 比对到基因组的bam文件
# trt_N061011.Aligned.toTranscriptome.out.bam：比对到转录组的bam文件，供RSEM计算TPM使用 
# trt_N061011.Log.final.out: reads比对日志
# trt_N061011.Log.out: 运行参数和过程
# trt_N061011.Log.progress.out: 运行日志

# trt_N061011.ReadsPerGene.out.tab: 每个基因的reads count，链非特异性RNASeq选第2列.
# column 1: gene ID
# column 2: counts for unstranded RNA-seq
# column 3: counts for the 1st read strand aligned with RNA (htseq-count option -s yes)
# column 4: counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)
# 
# Select the output according to the strandedness of your data. Note, that if you have stranded data an
d choose one of the columns 3 or 4, the other column (4 or 3) will give you the count of antisense read
s. 

# trt_N061011.SJ.out.tab: Junction reads
# trt_N061011.Unmapped.out.mate1：未比对上的reads
# trt_N061011.Unmapped.out.mate2：未比对上的reads

Kallisto和Pseudo-Alignment

STAR是序列比对工具，而Kallisto是伪比对工具 [@bray_2016]。它们的区别是：比对工具是把reads比对回参考基因组或转录组，而伪比对工具是把k-mers比对到参考转录组。

什么是k-mer？

k-mer是来源于测序短序列中的长度为k的子序列。例如，假设有短序列ATCCCGGGTTAT，想从中获得7-mer。为此，我们将提取前七个碱基作为第一个7-mer，然后向下移动一个碱基获得第二个7-mer，以此类推。下图显示了从序列中可以得到的所有7-mers：

为什么比对k-mers而不是reads？

有两个主要原因：

伪比对工具利用算法技巧使得比对k-mers比比对reads速度快很多，具体可以参阅 [@bray_2016]了解详情。
在某些情况下，伪对齐工具可能比传统比对工具更好的处理测序错误问题。例如，假设序列上第一个碱基中存在测序错误，本来是T却测序成了A。对伪比对工具来说，只会影响第一个7-mer而不会影响后续7-mer的比对。

Kallisto的伪比对模式

Kallisto有一个为单细胞转录组特别设计的伪比对模式。和STAR不同，Kallisto比对到的是参考转录组而不是参考基因组，意味着Kallisto是将序列比对到剪接异构体而不是基因上，对单细胞来讲，这是有挑战性的：

单细胞RNA-seq的覆盖率低于bulk RNA-seq，这意味着可以从序列中获得的信息总量减少了。
许多单细胞RNA-seq方案具有3'覆盖偏好性，如果两种剪接异构体只在5’末端不同，则很难确定序列来源于哪个剪接异构体。
一些单细胞RNA-seq方案测序读长短，也难以区分来源于哪个剪接异构体。

Kallisto的伪模式采用了略微不同的伪比对方法。它不与剪接异构体比对，而是与等价类 (equivalence classes)比对。所以如果read比对到多个异构体，Kallisto会记录read比对到包含有所有异构体的等价类，因此可以使用等价类计数而非基因或转录本计数用于下游的聚类等分析。具体见下图：

今天我们只用一个细胞来展示伪比对，但是Kallisto其实可以同时对大量细胞进行伪比对，并且可以使用UMIs的信息。具体细节请参看https://pachterlab.github.io/kallisto/manual。

类似于STAR，伪比对前也需要为Kallisto生成一个索引。

任务6：用下面的命令生成Kallisto的索引，从Kallisto手册https://pachterlab.github.io/kallisto/manual中找到这个命令各个选项的含义。

mkdir indices/Kallisto
kallisto index -i indices/Kallisto/transcripts.idx Share/2000_reference.transcripts.fa

任务7: 从Kallisto手册找到实现伪比对的命令，如果你觉得知道了答案，可以在下一节核查它是否正确并执行。

Kallisto伪比对命令

mkdir results/Kallisto
kallisto pseudo -i indices/Kallisto/transcripts.idx -o results/Kallisto -b batch.txt

参考https://pachterlab.github.io/kallisto/manual构建batch.txt。

batch.txt文件格式如下：第一列是细胞的名字ID, 后面两列是对应的PE reads。所有#开头的行都忽略。如果是单端测序，命令行参数指定为--single -l -s，batch.txt只需提供一个fastq文件。

#id file1 file 2
cell1 cell1_1.fastq.gz cell1_2.fastq.gz
cell2 cell2_1.fastq.gz cell2_2.fastq.gz
cell3 cell3_1.fastq.gz cell3_2.fastq.gz
...

如果使用了--umi参数，batch.txt文件格式如下：

#id umi-file file-1
cell1 cell_1.umi cell_1.fastq.gz
cell2 cell_2.umi cell_2.fastq.gz
cell3 cell_3.umi cell_3.fastq.gz
...

umi-file是文本文件，每行是一个read的UMI序列，必须与fastq中reads的顺序一致。即便UMI数据是单端，也不会用到测序片段的长度信息。

TTACACTGAC
CCACTCTATG
CAGGAAATCG
...

kallisto在UMI模式下，不是计算每个equivalence class的reads count，而是使用测序reads做pseudoalign获得equivalence class, 计算每个equivalence class上的UMIs个数。

理解Kallisto的输出结果

上面命令会生成4个文件-matrix.cells, matrix.ec, matrix.tsv 和 run_info.json.

matrix.cell 为细胞ID列表，测试数据只用了一个细胞，这个文件应该只包含“ERR522959”
matrix.ec 含有所有用到的等价类信息。每一行的第一个数字是等价类ID，第二个数字对应等价类的转录文本ID。举个例子，10 1,2,3表示等价类10包含转录文本ID 1, 2和3。ID号对应于转录本在reference.transcripts.fa中出现的顺序。第一个出现的转录本ID为0，转录本ID 1, 2和3对应于reference.transcripts.fa中的第二，第三和第四条转录本。
matrix.tsv 表示每个等价类在不同细胞中的reads count信息。第一个数字是等价类ID，和matrix.ec中的定义一样。第二个数字是细胞ID，和matrix.cell文件中的细胞名字顺序一一对应。第三个数字是等价类在该细胞中的reads count。举个例子，5 1 3指细胞1中有3个reads比对到等价类5。细胞ID同样是零起始的索引，所以细胞1指matrix.cells中的第二个细胞。
run_info.json 含有关于Kallisto执行的信息，一般可忽略。