Gene expression levels analysis -- HTSeq-count

前言

HTSeq 是一款针对于有参考基因组的,分析高通量测序数据的Python库。
官网: HTSeq


安装

HTSeq作为Python库推荐采用pip安装方式:

1
pip install htseq

还能采用源码安装、conda安装、apt-get等

使用

1
htseq-count [options] <alignment_file> <gtf_file>

这里的 alignment file 可以是 bam 或者 sam 文件

上游准备

  1. 已排序的比对文件: sam 或者 bam 文件
  2. 排序要求建议按 reads name 排序
  3. 计数模型:union, intersection-strict and intersection-nonempty
    Count_model
  4. 测序的建库工作是否为链特异性建库,一般不是
  5. 参数 type 和 idattr 的确定,非 Ensemble 的注释文件这两个参数尚待明确,知道的麻烦告知一下,提前感谢!

    下游适配

    HTSeq 的计数结果没有进行标准化,因为我下游分析是采用 DESeq2,不需要标准化

    常用参数

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    # 命令参数
    -f | --format default: sam 设置输入文件的格式,该参数的值可以是sam或bam。
    -r | --order default: name 设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
    -s | --stranded default: yes 设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。
    -a | --a default: 10 忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。
    -t | --type default: exon 程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
    -i | --idattr default: gene_id 设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
    -m | --mode default: union 设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
    -o | --samout 输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
    -q | --quiet 不输出程序运行的状态信息和警告信息。
    -h | --help 输出帮助信息。

批量化处理

借助bash:

1
2
3
4
5
6
7
8
9
#!/bin/bash

for i in 10 11; do
{
htseq-count -f bam -s no in/E${i}.sorted.bam\
anno/Gallus_gallus.Gallus_gallus-5.0.94.gtf > out/E${i}_count.txt
};
done
wait

此处参考基因组来自于 Ensemble,而 -t 和 -i 在采用 Ensemble 的注释文件时默认为:exon 和 gene_id

输出

结果文件为:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
head counts.txt
ENSG00000000003 0
ENSG00000000005 0
ENSG00000000419 1171
ENSG00000000457 563
ENSG00000000460 703
ENSG00000000938 0
ENSG00000000971 1
ENSG00000001036 925
ENSG00000001084 1468
ENSG00000001167 2997
 
tail count.txt
ENSG00000283696 18
ENSG00000283697 0
ENSG00000283698 1
ENSG00000283699 0
ENSG00000283700 0
__no_feature    3469791
__ambiguous 630717
__too_low_aQual 1346501
__not_aligned   520623
__alignment_not_unique  2849422

下游分析需要将其合并到一个文件中,一个行为基因名,列为样本名,中间为 count 的行列式文件。这里也可以将其中 count 为零的行删除。所以应该是先合并、再删除计数全为零的行。

参考

  1. RNA-seq分析htseq-count的使用
  2. RNA-seq流程分析比较之半小时得到差异基因
  3. 转录组入门(6): reads计数
  4. 官网
-------------本文结束 感谢您的阅读-------------
暖一下
ZJohnson wechat
扫一扫,领红包!
0%