mRNA-seq学习（三）：htseq-count计数

htseq-count计数的相关内容前面在不同的学习阶段写过两次，分别是合并htseq-count的结果得到count matrix和htseq-count的一个坑，其中第二篇中关于“坑”的总结我觉得还是挺用的。

1. 基因表达定量的三个水平

• 基因
• 转录本
• 外显子

基因和外显子的定义明晰，统计起来相较于转录本简单；而由于不同的转录本往往存在外显子的重叠，统计起来就比较困难了。基因水平的定量常见。

2. 四种不同的reads计数思路

1. 当比对到有注释的基因组时，基于注释文件统计reads
2. 基于参考基因组的转录本组装时，如Cufflinks会提供注释文件, 且能够发现新的基因和转录本。这种情况下，也要结合软件给的注释文件计数
3. 比对到转录本序列可以直接计数，不借助注释文件
4. 重头组装出转录本序列，接下来同3

3. 哪些因素影响了feature内的reads数

• 测序深度
• feature长度
• feature复杂度
• GC偏好
• 序列特异偏好

常将前两者考虑在标准化之内

4. 关于HTSeq

4.1 如何处理多比对reads

HTSeq忽略掉这些多比对reads

4.2 HTSeq的计数模式

default: union

4.3 HTSeq的使用

usage: htseq-count [options] alignment_file gff_file

-f {sam,bam}  (default: sam)
-r {pos,name}  (default: name)
-s {yes,no,reverse}  (default: yes) #此处关于选项-s为我自己的认识，不一定对
    #数据是否来源于链特异性测序，链特异性是指在建库测序时，只测mRNA反转录出的cDNA序列，而不测该cDNA序列反向互补的另一条DNA序列；换句话说就是，链特异性能更准确反映出mRNA的序列信息
    #我们知道在gff/gtf中第7列是+-信息，+表示来源于参考基因组序列正链，-表示参考基因组序列的反向互补链
    #sam/bam文件的第2列是flag信息，也可以看出比对到正链还是负链
    #stranded=no，无链特异性，一条reads通过flag列知道比对到+还是-链后，不管是不是和gff/gtf相匹配，都算是这个feature中的
    #stranded=yes, 且se测序，要和gff/gtf相匹配，才算是这个feature中的
    #stranded=yes, 且pe测序，要和gff/gtf相匹配，才算是这个feature中的
    #stranded=reverse，是yes的相反，这时不是和gff/gtf相匹配了，而是恰好相反，可能源于另一种链特异性，只测cDNA序列反向互补的另一条DNA序列
-a MINAQUAL (default: 10)
    #忽略比对质量低于此值的比对结果
-t FEATURETYPE 
    #feature type (3rd column in GFF file) to be used, all features of other type are ignored (default, suitable for Ensembl GTF files: exon)
    #没想到这个还能自己设置
-i IDATTR
    #GFF attribute to be used as feature ID (default, suitable for Ensembl GTF files: gene_id)
-m {union,intersection-strict,intersection-nonempty} (default: union)

4.4 HTSeq输出结果

$ ls *count
SRR3286802.count  SRR3286803.count  SRR3286804.count  SRR3286805.count  SRR3286806.count  SRR3286807.count

#基于相同gff/gtf得到的计数文件，行数相同，第一列（基因名）相同
$ wc -l *count
  37889 SRR3286802.count
  37889 SRR3286803.count
  37889 SRR3286804.count
  37889 SRR3286805.count
  37889 SRR3286806.count
  37889 SRR3286807.count

#且最后5列统计了整个计数过程没有使用到的reads
$ tail -n 5 SRR3286802.count
__no_feature    237560
__ambiguous 1846779
__too_low_aQual 0
__not_aligned   1323985
__alignment_not_unique  2015872

• based on the NH tag in the BAM file, they aligned to more than one place in the reference genome (alignment_not_unique);
• they did not align at all (not_aligned);
• their alignment quality was lower than the user-specified threshold (too_low_aQual);
• their alignment overlapped with more than one gene (ambiguous);
• their alignment did not overlap any gene (no_feature).

4.5 什么情况下使用

因为是基于gff/gtf的feature来计数，所以比对策略应该是往参考基因组上比对。