番茄转录组分析实战：表达定量（1）-HTSeq

在上一讲番茄转录组分析实战：表达定量概述中，我们提到HTSeq可以用于基因的表达定量。下面我就简单介绍一下这个软件。

HTSeq是一个python包，所以可以直接整合到python流程中，但今天我们直接在命令行中使用其中的可执行文件。如下图所示，HTSeq有多种功能，但今天我们只用其中的htseq-count来进行基因定量。

这里需要注意的是，htseq-count用法是，基于比对后的测序数据和基因feature列表，数有多少reads分别比对到每个feature上。

feature这里指的是染色体中的一段区域或多段区域。RNA-Seq中一般是指genes。下面这张图直观的解释了HTSeq是如何分配reads到各个基因的。首先，HTSeq有3种模式，分别为union/intersection_strict/intersection_nonempty，默认是模式，也是大部分情况下常用的模式。如果这三种模式还不能满足你的需求，可以自己编写模式，具体看官网上的文档。

但是，htseq的几个模式都有一点问题，就是上图的最后三行。当一条read遇到重叠基因或者说这个read能比对到基因组上的多个位置时，HTSeq就糊涂了，它干脆就将这些reads标为ambiguous。对待ambiguous reads，htseq可以通过选项，选择不同的处理方法。一是干脆将ambiguous reads都删去，都不计数，二是将ambiguous reads能够map到的每个基因的计数都加1。那这个参数究竟怎么选呢？作者给出的观点是默认参数，也就是直接将ambiguous reads都删除。

htseq-count使用

This script takes one or more alignment files in SAM/BAM format and a featurefile in GFF format and calculates for each feature the number of reads mappingto it. See http://htseq.readthedocs.io/en/master/count.html for details.

htseq-count的主要参数如下：

-fdefault: sam输入文件格式，sam或者bam-rdefault: namesam或者bam输入文件的排序方式，参数可以是name或者pos，name表示按read名进行排序，pos表示按比对的参考基因组位置进行排序,如果按name排序，则read pair出现在临近位置，能够节省内存，因为程序会将其中一个read放入内存中，直到找到另一端的read-sdefault: yes设置是否为链特异性测序，值有yes、no、reverse。"no"表示不管read是比对到基因所在链还是相反的链，对进行计数；"yes"表示，(1)在单端测序中，read需要比对到与基因相同的链，双末端测序中，第一条read需要比对到与基因相同的链，另一条read需要比对到与基因相反的链；"no"与"yes"相反；这里需要注意的是，一般的应该选"no"-----------------------------adefault: 10忽略低于指定值的比对结果-tdefault: exongtf或者gff第三列的feature(如exon、gene等)进行表达量计算，指定后其他feature将被忽略，在RNA-Seq中，对ensembl的gtf文件最好使用"exon"；-idefault: gene_idgff或者gtf中可以用作feature ID的属性--additional-attrdefault：none除了-i属性外，还可以设置另一个feature ID-mdefault: union设置表达量计算模式，union/intersection_strict/intersection_nonempty--nonuniquedefault: none如何计算比对到多个特征的reads，具体见上文-o输出所有alignment的reads到一个sam文件中，通过一个可选的sam标签 ‘ XF ’来标注每一行对应的单位和计数，可以不设置

HTSEQ_COUNT=/public/ptbus/home/zhumy/software/anaconda3/bin/htseq-count

GTF="$WORKING_DIR/genes.gtf"

MAPPING=$WORKING_DIR/result/mapping

COUNT=$WORKING_DIR/result/count

mkdir -p$COUNT

cd$MAPPING

fori in`ls$MAPPING`；

do

filename=${i%.*}

samtools view -bS$isamtools sort - -n -o${filename}.sorted_n.bam

$HTSEQ_COUNT-f bam -s no -t gene${filename}.sorted_n.bam$GTF>$COUNT/${i}.htseq.count

done

上面是htseq输出文件中，最后的一部分统计量，每个含义如下图所示：

HTSeq的官网是『http://htseq.readthedocs.io/en/release_0.9.1/』，发表的文章是『HTSeq — A Python framework to work with high-throughput sequencing data』下面是HTSeq官网给出的使用者经常遇到的问题，红色框部分尤其重要。