首页
学习
活动
专区
工具
TVP
发布

番茄转录组分析实战:表达定量(1)-HTSeq

在上一讲番茄转录组分析实战:表达定量概述中,我们提到HTSeq可以用于基因的表达定量。下面我就简单介绍一下这个软件。

HTSeq是一个python包,所以可以直接整合到python流程中,但今天我们直接在命令行中使用其中的可执行文件。如下图所示,HTSeq有多种功能,但今天我们只用其中的htseq-count来进行基因定量。

这里需要注意的是,htseq-count用法是,基于比对后的测序数据和基因feature列表,数有多少reads分别比对到每个feature上。

feature这里指的是染色体中的一段区域或多段区域。RNA-Seq中一般是指genes。下面这张图直观的解释了HTSeq是如何分配reads到各个基因的。首先,HTSeq有3种模式,分别为union/intersection_strict/intersection_nonempty,默认是模式,也是大部分情况下常用的模式。如果这三种模式还不能满足你的需求,可以自己编写模式,具体看官网上的文档。

但是,htseq的几个模式都有一点问题,就是上图的最后三行。当一条read遇到重叠基因或者说这个read能比对到基因组上的多个位置时,HTSeq就糊涂了,它干脆就将这些reads标为ambiguous。对待ambiguous reads,htseq可以通过选项,选择不同的处理方法。一是干脆将ambiguous reads都删去,都不计数,二是将ambiguous reads能够map到的每个基因的计数都加1。那这个参数究竟怎么选呢?作者给出的观点是默认参数,也就是直接将ambiguous reads都删除。

htseq-count使用

This script takes one or more alignment files in SAM/BAM format and a featurefile in GFF format and calculates for each feature the number of reads mappingto it. See http://htseq.readthedocs.io/en/master/count.html for details.

htseq-count的主要参数如下:

-fdefault: sam输入文件格式,sam或者bam-rdefault: namesam或者bam输入文件的排序方式,参数可以是name或者pos,name表示按read名进行排序,pos表示按比对的参考基因组位置进行排序,如果按name排序,则read pair出现在临近位置,能够节省内存,因为程序会将其中一个read放入内存中,直到找到另一端的read-sdefault: yes设置是否为链特异性测序,值有yes、no、reverse。"no"表示不管read是比对到基因所在链还是相反的链,对进行计数;"yes"表示,(1)在单端测序中,read需要比对到与基因相同的链,双末端测序中,第一条read需要比对到与基因相同的链,另一条read需要比对到与基因相反的链;"no"与"yes"相反;这里需要注意的是,一般的应该选"no"-----------------------------adefault: 10忽略低于指定值的比对结果-tdefault: exongtf或者gff第三列的feature(如exon、gene等)进行表达量计算,指定后其他feature将被忽略,在RNA-Seq中,对ensembl的gtf文件最好使用"exon";-idefault: gene_idgff或者gtf中可以用作feature ID的属性--additional-attrdefault:none除了-i属性外,还可以设置另一个feature ID-mdefault: union设置表达量计算模式,union/intersection_strict/intersection_nonempty--nonuniquedefault: none如何计算比对到多个特征的reads,具体见上文-o输出所有alignment的reads到一个sam文件中,通过一个可选的sam标签 ‘ XF ’来标注每一行对应的单位和计数,可以不设置

HTSEQ_COUNT=/public/ptbus/home/zhumy/software/anaconda3/bin/htseq-count

GTF="$WORKING_DIR/genes.gtf"

MAPPING=$WORKING_DIR/result/mapping

COUNT=$WORKING_DIR/result/count

mkdir -p$COUNT

cd$MAPPING

fori in`ls$MAPPING`;

do

filename=${i%.*}

samtools view -bS$isamtools sort - -n -o${filename}.sorted_n.bam

$HTSEQ_COUNT-f bam -s no -t gene${filename}.sorted_n.bam$GTF>$COUNT/${i}.htseq.count

done

上面是htseq输出文件中,最后的一部分统计量,每个含义如下图所示:

HTSeq的官网是『http://htseq.readthedocs.io/en/release_0.9.1/』,发表的文章是『HTSeq — A Python framework to work with high-throughput sequencing data』下面是HTSeq官网给出的使用者经常遇到的问题,红色框部分尤其重要。

  • 发表于:
  • 原文链接http://kuaibao.qq.com/s/20180201G18QSE00?refer=cp_1026
  • 腾讯「腾讯云开发者社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。
  • 如有侵权,请联系 cloudcommunity@tencent.com 删除。

扫码

添加站长 进交流群

领取专属 10元无门槛券

私享最新 技术干货

扫码加入开发者社群
领券