对基因进行定量,首先需要计算出比对到各个基因的read counts,得到没有校正的表达矩阵。
在基因层面的count计数的工具:Htseq-count,feature-count
1、Htseq-count
官网
https://htseq.readthedocs.io/en/release_0.11.1/count.html
它数count的模式
从它的第七种情况看出来:它没有办法解决基于转录本的定量,转录本存在可变剪接,这种overlap的情况很容易出现。
安装Htseq-count
conda install htseq
htseq-count -h #查看是否安装成功
调用Htseq-count
htseq-count -f bam -r pos --max-reads-in-buffer 1000000 \
--stranded no --minaqual 10 --type exon \
--idattr gene_id --mode union \
--nonunique none --secondary-alignments ignore \
--supplementary-alignments ignore \
--counts_output htseq-count_out/test_count.tsv \
tophat_out/accepted_hits.bam \
ref/gft/Homo_sapiens.GRCh38.108.chr.gtf
# -f 输入文件格式sam,bam
# -r 输入文件的排序方式
# --stranded no 链特异性
# --minaqual 10 错误率为1%
# --type exon 指定GTF注释中的要素类型
# --idattr gene_id 输出文件的行名
# --mode union 数count的模式
# --nonunique none 对未唯一对齐或分配不明确的reads的计数,全部不要
# --secondary-alignments ignore 对第二优的对齐的计数忽略
# --supplementary-alignments ignore 是否进行补充比对
#--counts_output 输出文件
#bam文件
#gtf文件
2、featureCounts
featureCounts: a ultrafast and accurate read summarization program
官网
https://subread.sourceforge.net/featureCounts.html
featureCounts是一个高效的通用reads统计程序,它对基因组特征如基因、外显子、启动子进行read counts计数。它可用用于对RNA-seq和DNA-seq的reads进行计数。
安装featureCounts
conda install subread
featureCounts -h #查看是否安装成功
调用featureCounts
featureCounts -t exon -g gene_id -Q 10 --primary \
-s 0 -p -T 1 \
-a ref/gft/Homo_sapiens.GRCh38.108.chr.gtf \
-o featurecounts_out/featurecounts.tsv \
tophat_out/accepted_hits.bam
# -t exon 指定GTF注释中的要素类型,这里指定的是外显子
# -g gene_id 输出文件的行名
# -Q 10 错误率10%
# --primary 只数最优的
# -s 0 链特异性
# -T 1 CPU
# -a gtf文件
# -o 输出文件
# 可以输入多个bam文件,会直接整合到一个表达矩阵
表达矩阵
可以提取我们需要的数据,gene_id 和count数
featureCounts比htseq-count快很多,另外featureCounts可以一次性对多个Bam文件进行操作。