STAR直接就可以输出readsCount，为什么还需要featurecounts？

这个问题很让人困惑，不少教程，先是STAR比对，然后featureCounts或HTSeq再计算reads count。那么我们看看，什么时候需要这样做，什么时候不需要这样做？

STAR比对可以直接输出reads count

STAR比对参数很多，其中有一个quantMode，可以指定--quantMode GeneCounts输出STAR计算出的reads计数结果。(更多常用参数见 STAR比对线程也不是越多越好，多少是好？)

格式如下，有4列，各自解释如下：

trt_N061011.ReadsPerGene.out.tab: 每个基因的reads count，链非特异性RNASeq选第2列。

column 1: gene ID

column 2: counts for unstranded RNA-seq

column 3: counts for the 1st read strand aligned with RNA (htseq-count option -s yes)

column 4: counts for the 2nd read strand aligned with RNA (htseq-count option -s reverse)

N_unmapped    127    127    127
N_multimapping    3745    3745    3745
N_noFeature    21487    234292    234796
N_ambiguous    23935    5678    5571
ENSG00000178591    0    0    0
ENSG00000125788    0    0    0
ENSG00000088782    0    0    0
ENSG00000185982    0    0    0
ENSG00000125903    0    0    0
ENSG00000186458    0    0    0
ENSG00000272874    0    0    0
ENSG00000196476    71    33    38

这个结果与HTSeq的输出结果是完全一致的。所以我们如果是比对完之后未做转录本拼装，直接对已知基因（构建基因组索引时GTF中囊括的基因）进行定量时，完全不需要再次用featureCounts或HTSeq再计算reads count。

如果做了转录本拼装，对基因定量，需要HTSEQ/FEATURECOUNTS

假设我们用STAR比对后，做了转录本拼装，获得一个拼装比较后的注释文件assembeCompare2Ref.annotated.gtf，对新基因或有新转录本的基因进行定量时，就需要再次计算了。

下面代码显示的是htseq的计算方式，featureCounts改一下也适用。

for i in `tail -n +2 sampleFile | cut -f 1`; do 
    htseq-count -f bam -r pos -a 10 -t exon -s no -i gene_id -m union ${i}/${i}.Aligned.sortedByCoord.out.bam assembeCompare2Ref.annotated.gtf >${i}/${i}.readsCount
    grep -v '^__' ${i}/${i}.readsCount | sed "1 iGene\t${i}"  >${i}/${i}.readsCount2
done &

用软件，看别人分享的教程是一个快捷的方式。但也要看软件手册，看看哪些参数适合自己的物种、数据，再进行分析，不能别人用什么、怎么用，自己完全跟随。

即便是看完手册后，你发现不需要额外设置参数，直接套教程的代码就行，那也需要看手册！

看了才有底气，指导你的数据用这个参数是合理的，而不是懵懵懂懂、人云亦云~~~

• 给地球上的22.8万种植物建立图像索引

• 公众号标题变迁史 - 你有没有被套路？

• 生信学习学的是什么？常识！

• 学生信，不是贪多的，而是求精的！

• 生信学习入门常见错误可能的原因分类总结和求助指南

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集