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给出的比对结果，而且必须为依据reads名称排序的sort bam文件。
所以在之前需应用samtools对bam文件进行排序。
nohup cat ../SRR_Acc_List.txt | while read line; do samtools sort -o $line\_sorted.bam -n -T sorted $line.bam; done &（依然是批量后台运行）
安装
使用conda进行安装 conda install -c bioconda htseq=1.99.2
定量
HTSeq计算counts数有3种模式，如下图所示（ambiguous表示该read比对到多个gene上；no_feature表示read没有比对到gene上）：
htseq-count对链特异性及非链特异性数据的处理是不同的因此在处理数据前一定要明确目标数据的建库方式（链特异性或非链特异性的）。明确这一关键问题之后就可以进行定量啦！
#htseq-count常用参数
-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
vim htseq.sh
cat SRR_Acc_List.txt | while read line
do
htseq-count -f bam -s no -m intersection-nonempty 2-mapping/$line\_sorted.bam ../../reference/homo/Homo_sapiens.GRCh38.104.gtf > 3-quantity/$line\_count.txt
done
#后台运行bash
nohup bash htseq.sh &
结果如下
接下来我们利用R语言对结果进行合并，先看一下文件内容
##看一下任意文件前十行
head SRR14760842_count.txt
ENSG00000000003 6157
ENSG00000000005 0
ENSG00000000419 2925
ENSG00000000457 791
ENSG00000000460 910
ENSG00000000938 0
ENSG00000000971 2192
ENSG00000001036 3245
ENSG00000001084 1427
ENSG00000001167 1897
#看一下任意文件后十行
tail SRR14760842_count.txt
ENSG00000288721 6
ENSG00000288722 166
ENSG00000288723 1
ENSG00000288724 0
ENSG00000288725 1
__no_feature 1980267
__ambiguous 1633742
__too_low_aQual 821454
__not_aligned 961622
__alignment_not_unique 1571191
观察数据我们可以发现，最后五行是不需要的因此需要去除，此外文件缺少表头信息，所以在处理时要注意标注，同时我们还需要对ensembleid进行注释。
上代码！
setwd("F:\\RNA-Seq\\GSE176393（SRP323246）\\3-quantity")
library(magrittr)
library(clusterProfiler)
library(org.Hs.eg.db)
library(limma)
#读取文件
files = list.files(pattern = "_count")
file <- list()
for (i in files) {
file[[i]]<-read.table(i,sep = "\t",check.names = F,col.names = c("gene_id",i))
}
rt = Reduce(function(x,y) merge(x,y,by="gene_id",all.x=TRUE),file,accumulate=FALSE)
#注释
rt <- rt[-1:-5,]
ensembl = rt[,1]
geneSymbol <- bitr(ensembl, fromType = "ENSEMBL", #fromType为输入数据ID类型
toType = "SYMBOL", #toType是指目标类型
OrgDb = org.Hs.eg.db)#Orgdb是指对应的注释包，这是人的注释包，相关注释包请更换
rt <- merge(geneSymbol,rt,by.x="ENSEMBL",by.y="gene_id")
#对基因名称重复的reads count进行处理
rt<-rt[,-1]
rt<-as.matrix(rt)
rownames(rt)<-rt[,1]
rt<-rt[,2:ncol(rt)]
rt<-avereps(rt)
rt1<-apply(rt,2,as.numeric)
rownames(rt1)<-rownames(rt)
rt1<-as.data.frame(rt1)
rt1<-cbind(row.names(rt1),rt1)
colnames(rt1)[1]<-"genes"
#输出
write.table(rt1,"rawread.txt",sep = "\t",row.names = F,col.names = T,quote = F)
最后就获得了合并的read count文件
genes SRR14760842_count.txt SRR14760843_count.txt SRR14760844_count.txt SRR14760845_count.txt SRR14760846_count.txt SRR14760847_count.txt SRR14760848_count.txt SRR14760849_count.txt SRR14760850_count.txt
TSPAN6 6157 5622 6525 4781 4797 4543 3796 3663 4326
TNMD 0 0 0 0 0 0 0 0 0
DPM1 2925 2901 2703 1977 2176 2214 2166 2153 2416
SCYL3 791 764 681 627 575 576 295 336 353
C1orf112 910 806 726 527 545 550 521 587 620
FGR 0 0 0 1 0 0 0 0 0
CFH 2192 2064 2219 1035 958 1036 245 260 317
FUCA2 3245 3049 3256 2229 2340 2095 2636 2942 3107
GCLC 1427 1428 1495 1390 1509 1281 3693 3543 3891
NFYA 1897 1747 1954 1605 1624 1551 1606 1438 1619
STPG1 407 395 403 545 540 545 489 421 569
NIPAL3 1570 1604 1839 2174 2263 2136 1309 1327 1575
LAS1L 2594 2444 2296 1770 1793 1916 1878 2103 2131
ENPP4 2015 2023 1998 2121 2189 2065 1867 1625 1916
SEMA3F 131 158 136 175 157 165 116 124 193
CFTR 310 341 331 161 119 126 39 35 38
