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