HTSeq-count是用于NGS数据分析一个Python脚本包的一部分，但是它的使用不需要Python的任何知识。Hseq-count需要比对文件为SAM/BAM格式，基因组注释文件为GFF/GTF文件。请注意，为了使Read的比对位置与基因组特征匹配，比对文件和注释文件必须具有相同的染色体名称。HTSeq-count先找出reads比对上的外显子，然后根据GTF文件中外显子的基因ID统计外显子水平的表达量。这就要求基因的所有外显子都具有相同的基因ID。Ensembl数据库中的GTF文件遵循这个规则，但是UCSC数据库中的GTF文件把转录本ID作为基因ID。这会导致Htseq-count计数错误，因为它不知道哪些转录本属于同一个基因，因此会分开计算不同名的转录本的表达量。Ensembl数据库的GTF文件可在

http://www.ensembl.org/info/data/ftp/index.html

获得，选择目的物种的“GTF”文件类型进行下载即可。在下面的例子中，我们使用双端测序比对文件，所以我们要下载人的GTF文件:

wget ftp://ftp.ensembl.org/pub/release74/gtf/homo_sapiens/Homo_sapiens.GRCh37.74.gtf.gz

然后解压文件：

gunzip Homo_sapiens.GRCh37.74.gtf.gz

HTSeq-count默认双端测序比对数据按照reads名称进行排序，以便双端比对数据在文件中前后跟随。也可以按照基因组的位置对比对文件进行排序（使用option-order = pos），但是这样要占用更多的内存。以下命令根据read名称对BAM进行排序，并生成文件hits_namesorted.bam

samtools sort –n accepted_hits.bam hits_namesorted

用htseq-count命令对比对文件进行计数的命令如下所示：

htseq-count –f bam --stranded=no hits_namesorted.bam Homo_sapiens.GRCh37.74.gtf > counts.txt

这里-f bam表示输入格式是BAM。默认参数下，先计算与GTF文件（--type = exon）中外显子位置相匹配的reads数，再统计属于同一基因（--idattr = gene_id）的外显子的reads数。Htseq-count默认数据是链特异性建库的，因此只计算落在与基因同一条DNA链上的reads。由于示例数据并非链特异性的，因此我们必须添加--stranded = no参数，以便在计算时包括比对到相反链的reads数。默认的计数模式是联合，但是你可以用--mode选项来改变它。你也可以设置一个比对质量阈值（例如-a 30）来过滤低比对质量的reads，该默认值是10。输出的counts.txt是一个包含每个基因的计数值的表格。在文件末尾，有五行列出了因为一下情况而不被包括到任何基因的reads数目。这五种情况是：

a. 根据BAM文件中的NH标签，reads被比对到参考基因组的多个位置上（alignment_not_unique，比对位置不单一）;

b. 根本无法被比对上（not_aligned）;

c. 其比对质量低于用户指定的阈值（too_low_aQual）;

d. 比对上的区域与多个基因重叠（ambigous）;

e. 比对上的区域没有与任何基因重叠（no_feature）。

可以用Linux命令join将不同样品的count.txt合并成一个表格：

join counts1.txt counts2.txt >count_table.txt

最后，在进行差异表达分析前，需要去掉最后5行，下面这条命令能去掉最后5行：

head -n -5 count_table.txt > genecounts.txt

欢迎关注

TCGA | 小工具 | 数据库 |组装| 注释 |   基因家族  |  Pvalue

基因预测  |bestorf |  sci | NAR | 在线工具 | 生存分析 | 热图

 生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos

 舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 |  进化 | 测序简史