转录组对基因表达量进行计数(Htseq、featureCounts)

对基因进行定量,首先需要计算出比对到各个基因的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文件进行操作。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容