这个问题很让人困惑,不少教程,先是STAR
比对,然后featureCounts
或HTSeq
再计算reads count
。那么我们看看,什么时候需要这样做,什么时候不需要这样做?
STAR比对可以直接输出reads count
STAR比对参数很多,其中有一个quantMode
,可以指定--quantMode GeneCounts
输出STAR
计算出的reads计数结果。(更多常用参数见 STAR比对线程也不是越多越好,多少是好?)
格式如下,有4列,各自解释如下:
代码语言:javascript复制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 iGenet${i}" >${i}/${i}.readsCount2
done &
用软件,看别人分享的教程是一个快捷的方式。但也要看软件手册,看看哪些参数适合自己的物种、数据,再进行分析,不能别人用什么、怎么用,自己完全跟随。
即便是看完手册后,你发现不需要额外设置参数,直接套教程的代码就行,那也需要看手册!
看了才有底气,指导你的数据用这个参数是合理的,而不是懵懵懂懂、人云亦云~~~