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（質量分析）。php

這裏須要注意的是HTSeq做爲read counts的計數軟件，承接的是上游比對軟件對於clean data給出的比對結果即bam文件（由sam文件sort獲得），和HTSeq能行使一樣做用的還有相似於GFold，bedtools等軟件，我會在最後作一個基本的結果比對。css

附manual

附油管視頻講解

HTSeq的安裝

 1 # 建立存放文件夾  2 mkdir ~/biosoft/HTseq && cd ~/biosoft/HTseq  3 
 4 # download並解壓  5 wget https://pypi.python.org/packages/fd/94/b7c8c1dcb7a3c3dcbde66b8d29583df4fa0059d88cc3592f62d15ef539a2/HTSeq-0.9.1.tar.gz#md5=fc71e021bf284a68f5ac7533a57641ac
 6 tar zxvf HTSeq-0.9.1.tar.gz  7 cd HTSeq-0.9.1/
 8 
 9 #使用python命令安裝，此處注意，install --user參數最好用上，除非你能夠獲取root權限 10 python setup.py build 11 python setup.py install --user 12 
13 # add bin/ to your PATH 14 vim .bashrc 15 PATH=/home/path_to/.local/bin:$PATH 16 source .bashrc

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

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

# 命令參數 
-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 對照示意圖能夠選擇本身須要的模式html

我這裏使用intersection_nonempty

mode

HTSeq的輸出

HTSeq將Count結果輸出到標準輸出，其結果示例以下：

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的計數。python

安裝

gfold.V1.1.4.tar.gzdownload解壓後便可使用ios

使用

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文件包含五列：git

#說明很詳細，這裏再也不翻譯

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文件示例：github

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的輸入文件一致，可是結果出入仍是較大的，此處僅做說明，不加以推薦。express

Bedtools ：再一個count matrix的提取工具

bedtools是一個極其老牌的數據處理軟件了，由猶他大學一個實驗室開發，我也是看了生信菜鳥團Jimmy的一篇文章才知道也能夠用來計數的。vim

安裝

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

使用

bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bed file.bed > read.count.txt

# 注意，此處的bed文件須要本身處理獲得的，須要四列，第一列爲chrN，第二列第三列爲基因位置，第四列爲基因名。相似於：
chr1 0   10000   ivl1
chr1 10000   20000   ivl2
chr1 20000   30000   ivl3
chr1 30000   40000   ivl4

RNA-seq分析htseq-count的使用