标签：count 文件 sam seq RNA 测序 HTSeq bam
HTSeq作为一款可以处理高通量数据的python包，由Simon Anders, Paul Theodor Pyl, Wolfgang Huber等人携手推出HTSeq — A Python framework to work with high-throughput sequencing data。自发布以来就备受广大分析人员青睐，其提供了许多功能给那些熟悉python的大佬们去自信修改使用，同时也兼顾着给小白们提供了两个可以拿来可用的可执行文件 htseq-count（计数） 和 htseq-qa（质量分析）。
这里需要注意的是HTSeq作为read counts的计数软件，承接的是上游比对软件对于clean data给出的比对结果即bam文件（由sam文件sort得到），和HTSeq能行使同样作用的还有类似于GFold，bedtools等软件，我会在最后做一个基本的结果比对。
附manual
附油管视频讲解
HTSeq的安装
HTSeq安装
HTSeq使用注意事项
- HTSeq是对有参考基因组的转录组测序数据进行表达量分析的，其输入文件必须有SAM和GTF文件。
- 一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析，而不是一个样品内部基因的表达量比较。因此，HTSeq设置了-a参数的默认值10，来忽略掉比对到多个位置的reads信息，其结果有利于后续的差异分析。
- 输入的GTF文件中不能包含可变剪接信息，否则HTSeq会认为每个可变剪接都是单独的基因，导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous，从而不能计算到基因的count中。即使设置-i参数的值为transcript_id，其结果一样是不准确的，只是得到transcripts的表达量。
HTSeq的使用
#这里承接的是上游hisat2比对软件得到的bam文件，sort by pos, 所以需要重新sort
|
1
2
3
|
samtools sort -n yourfile.bam > yourfile_name.bam
htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/reference/hisat2_reference/Homo_sapiens.GRCh38.86.chr_patch_hapl_scaff.gtf > counts.txt
|
|
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 输出帮助信息。
|
htseq-count 的三种比对模式
union, intersection-strict and intersection-nonempty 对照示意图可以选择自己需要的模式
我这里使用intersection_nonempty
mode
HTSeq的输出
HTSeq将Count结果输出到标准输出，其结果示例如下：
|
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
|
GFold：另一个count matrix的提取工具
GFold是一款2012年同济大学的研究组发表在Bioinformatics 上的软件，旨在通过对于相对基因变化找出RNA-seq中表达差异的基因，同时也可以用作read count的计数。
安装
gfold.V1.1.4.tar.gzdownload解压后即可使用
使用
|
1
2
|
gfold count -ann hg19Ref.gtf -tag sample1.sam -o sample1.read_cnt
gfold count -ann hg19Ref.gtf -tag sample2.sam -o sample2.read_cnt
|
输出
output文件包含五列：
#说明很详细，这里不再翻译
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
|
GeneSymbol:
For BED file, this is the 4'th column. For GPF file, this is the first column. For GTF format, this corresponds to 'gene_id' if it exists, 'NA' otherwise.
GeneName:
For BED file, this is always 'NA'. For GPF file, this is the 12'th column. For GTF format, this corresponds to 'gene_name' if it exists, 'NA' otherwise.
Read Count:
The number of reads mapped to this gene.
Gene exon length:
The length sum of all the exons of this gene.
RPKM:（#这里需要注意但是双端测序技术还未普及，这里未使用FPKM，况且RPKM和FPKM也不是能很好的代表基因表达水平 ）
The expression level of this gene in RPKM.
|
output文件示例：
|
1
2
3
4
5
6
7
8
9
10
11
|
head example.read_cnt
ENSG00000000003 TSPAN6 0 4535 0
ENSG00000000005 TNMD 0 1610 0
ENSG00000000419 DPM1 1588 1207 27.4411
ENSG00000000457 SCYL3 1344 6883 4.07267
ENSG00000000460 C1orf112 1334 5967 4.66292
ENSG00000000938 FGR 0 3474 0
ENSG00000000971 CFH 2 8145 0.0051215
ENSG00000001036 FUCA2 1427 2793 10.6564
ENSG00000001084 GCLC 2462 8463 6.06767
ENSG00000001167 NFYA 5123 3811 28.0378
|
此处使用示例bam文件or sam文件和HTSeq的输入文件一致，但是结果出入还是较大的，此处仅作说明，不加以推荐。
Bedtools ：再一个count matrix的提取工具
bedtools是一个极其老牌的数据处理软件了，由犹他大学一个实验室开发，我也是看了生信菜鸟团Jimmy的一篇文章才知道也可以用来计数的。
安装
|
1
2
|
wget https://github.com/arq5x/bedtools2/releases/download/v2.26.0/bedtools-2.26.0.tar.gz
tar zxvf bedtools-2.26.0.tar.gz
|
使用
|
1
|
bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bed file.bed > read.count.txt
|
|
1
2
3
4
5
|
# 注意，此处的bed文件需要自己处理得到的，需要四列，第一列为chrN，第二列第三列为基因位置，第四列为基因名。类似于：
chr1 0 10000 ivl1
chr1 10000 20000 ivl2
chr1 20000 30000 ivl3
chr1 30000 40000 ivl4
|
输出
标签: linux, RNA-seq
好文要顶 已关注 收藏该文 
标签：count,文件,sam,seq,RNA,测序,HTSeq,bam
来源： https://blog.csdn.net/u010608296/article/details/90208115
本站声明：
1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。
