查看原文
其他

转录组定量可以用替换featureCounts代替HTSeq-count

2018-01-10 anlan 生信技能树

在转录组定量分析时,如果采用的是alignment-based转录组定量策略,那么一般会使用的是HISAT2、STAR或者TopHat等比对软件。

接着则是对转录组进行定量,如果是基于基因水平的定量,我之前一般是采用HTSeq-count工具来获取每个基因上的count数。所谓count数,个人简单的理解为根据不同比对情况,将reads分配到各个基因上。HTSeq-count对于多重比对的reads则是采取舍弃策略。当HTSeq-count选择默认参数(-m 默认模式),那么reads是以下图所示的union的情况进行分配的




除了HTSeq-count工具外,其实也可以使用bedtools工具的multicov进行简单的基因水平定量。其需要一个所有基因的位置信息的bed文件,然后计算比对结果bam文件中的reads出现在基因interval上的个数,功能比较简单(说白了就是基于reads位置信息),适用性和准确性一般来说是没有HTSeq-count好。

以上两款在基因水平定量的软件有一个共同点,就是慢!如果样本数上升到几十及上百个样时,其所消耗的时间是以天记的,这时则可以考虑用featureCounts这款软件了,缺点还不知道,但是其优点已经听到过N遍了,就是快。。。

软件下载及安装

下载肯定去官网下载http://subread.sourceforge.net/

featureCounts: a ultrafast and accurate read summarization program

现在featureCounts已经整合在Subread里面了,粗略看了下简介和pdf文档,发现Subread是个功能很全面的软件,而且还有相对应的R包Rsubread,有机会再去瞧瞧。Subread有二进制版本,那么直接下载解压即可使用了

软件的使用

软件的作者认为其软件的优点在于(我就复制黏贴了):

  • It carries out precise and accurate read assignments by taking care of indels, junctions and structural variants in the reads

  • It takes only half a minute to summarize 20 million reads(真是快。。。)

  • It supports GTF and SAF format annotation

  • It supports strand-specific read counting

  • It can count reads at feature (eg. exon) or meta-feature (eg. gene) level

  • Highly flexible in counting multi-mapping and multi-overlapping reads. Such reads can be excluded, fully counted or fractionally counted(这点跟HTSeq-count不一样了,其对于多重比对的reads并不是只采用全部丢弃的策略,按照其说法是更加灵活的对待)

  • It gives users full control on the summarization of paired-end reads, including allowing them to check if both ends are mapped and/or if the fragment length falls within the specified range(可让使用者更加个性化的使用)

  • Reduce ambuiguity in assigning read pairs by searching features that overlap with both reads from the pair

  • It allows users to specify whether chimeric fragments should be counted(考虑的有点周到)

  • Automatically detect input format (SAM or BAM)

  • Automatically sort paired-end reads. Users can provide either location-sorted or namesorted bams files to featureCounts. Read sorting is implemented on the fly and it only incurs minimal time cost

featureCounts的参数有点多,可用-h先查看下

  1. ~/biosoft/subread/subread-1.6.0-Linux-x86_64/bin/featureCounts -h

按照官网教程简单使用下:

  1. ~/biosoft/subread/subread-1.6.0-Linux-x86_64/bin/featureCounts -T 6 -p -t exon -g gene_id -a ~/annotation/mm10/gencode.vM13.annotation.gtf -o SRR3589959_featureCounts222.txt SRR3589959.bam

主要的参数:

-a 输入GTF/GFF基因组注释文件

-p 这个参数是针对paired-end数据

-F 指定-a注释文件的格式,默认是GTF

-g 从注释文件中提取Meta-features信息用于read count,默认是gene_id

-t 跟-g一样的意思,其是默认将exon作为一个feature

-o 输出文件

-T 多线程数

其他参数介绍只能看文档了,不常用的话也是记不住的,要用 49 30754 49 15264 0 0 1126 0 0:00:27 0:00:13 0:00:14 2930再去翻就行

运行中和运行后有两张图可以看看,主要讲了其运行中的一些信息,如下:


featurecount2


从Multimapping reads和Multi-overlapping reads均为not counted可看出我的运行参数是不考虑multimapping的情况的,并且如果一个read在多个基因上有overlaps(multi-overlapping),则软件也是不计入count数的(只对RNA-SEQ数据,ChIP-seq数据则是例外)

Min overlapping bases为1则说明--minOverlap参数默认是1,因而只考虑全部overlap的reads

Successfully assigned fragments : 22730409 (59.9%)则说明只有59.9%的paired reads定量到了基因上,如果想具体看看那些没有分配上的reads是出于原因,则可看输出中的summary文件


featurecount3


上图的相关说明可以看软件的pdf文档,其实也蛮好理解的

Running time : 0.33 minutes只能说真是快。。。

SRR3589959_featureCounts222.txt文件包含了featureCount结果的信息,如果想了解每个基因上的count数,则只需要提取出第1列和第7列的信息

  1. cut -f 1,7 SRR3589959_featureCounts.txt |grep -v '^#' >feacturCounts.txt

然后统计下featureCounts结果总count数和HTSeq-count结果的总count数(这个之前转录组学习时已经跑过了)

  1. perl -alne '`$sum += $`F[1]; END {print $sum}' feacturCounts.txt

  2. #22730409

  3. grep -v '^_' HTSeq-count.txt |perl -alne '`$sum += $`F[1]; END {print $sum}'

  4. #21514247

从结果上看,两者数量相差还是蛮小的,看样子两者软件的默认参数的对待reads的方式还是差不多的(主要featureCounts的默认参数对待Multimapping reads是相同的,因为从图3可看出,MultiMapping的reads还是蛮多的)

最后用R作个散点图看看两者同个基因的reads数目大小的分布(取log2),并计算个相关系数

  1. library(ggplot2)

  2. data <- read.table(file = "count.txt", sep = "\t", header = T, row.names = 1, stringsAsFactors = F)

  3. data_matrix <- log2(data + 1)

  4. model <- lm(featurecounts ~ htseqcount, data_matrix)

  5. r2 <- summary(model)$r.squared

  6. r2 <- paste0("r^2=", r2)

  7. p <- ggplot(data_matrix, aes(x = featurecounts, y = htseqcount)) +

  8.  geom_point(colour = "grey60", size = 1) +

  9.  theme_light() +

  10.  stat_smooth(method = lm, se = FALSE) +

  11.  annotate("text", label = r2, x = 3.5, y = 1)

  12. p

从下图中可看出,两者的数据是呈正相关,两者绝大部分的counts数是非常接近的,在低表达量的那部分数据中HTSeq-count的值偏小于featureCounts。因而至少这个数据来说,我觉得featureCounts代替HTSeq-count进行表达量定量是完全没问题的。


featureCount_htseq-count



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存