转录组分析要用比对后的bam或sam文件进行后续计数。也就是对基因/转录本的区域上比对到的reads进行计数。最后算出来某个转录本或者外显子上比对上多少个reads。显然，如果某个样本中比对到某个基因转录本上的reads多，那这个基因的表达量理论上就高。在不同样本之间只用计数来比对似乎比较简单。但考虑到不同基因长度不一样，还有管家基因等，在计算表达矩阵的时候要更复杂些。
用于计数的软件有：htseq-count，featurecounts(subread)等。
一、featurecounts
这个用于计数的软件时集成在subread里面的。可以用作RNAseq数据计数，也可以用于ChIPseq的数据计数分析哦。运行它需要两种文件：
1，BAM/SAM作为样本输入，里面有每个reads比对到染色体上的位置信息。
2，GFF/GTF/SAF文件作为注释文件，里面有feature 的属性，染色体位置起始点这类信息。
注意的是：第二类文件要和前面BAM/SAM获得用的注释fa文件的版本来源一致。比如用了UCSC的hg38的fa文件去比对，那现在就要用UCSC的gtr来注释feature。（查看day41复习）
此外对feature还要理解一下，feature就是基因组上的一段序列，在一条染色体上的某一片段，有明确起点和终点，有明确的特定功能——也就是它的大小性质已被定义。那么mapping到某个feature上的reads多少，就能反应样本的某个特定功能的强或弱。metafeature就是一些feature组成的集合。比如exon可以说是feature，一个基因所有exon组成了一个metafeature。
1.单端数据
参数
-g 默认为gene_id，就是需要提供一个id identifier 来将feature水平的统计汇总为meta-feature水平的统计。此参数必须为gtf上有的列名 。
-t 默认为exon，就是read只有落到这些feature上才会被统计到。此参数必须为gtf上有的列名 。可以设定为-t gene，就不按照exon了吧。
-f 设定后统计feature层面的数据，如exon-level，否则会统计meta-feature层面的数据
conda activate py2.7 #激活2.7环境
project=~/rnaseq #项目文件夹
#step5：count
mkdir $project/count
output_dir=$project/count
GTF=/home/cloudam/reference/subread_index/hg38.knownGene.gtf
featureCounts -a $GTF -o $output_dir/counts.txt $project/align/*.bam
用UCSC里面的knownGene.gtf做出来的注释，基因名完全看不懂。运行速度倒是很快，每个样本都不到1min，最后生成的count.txt文件几十M，而且是包含了四个样本的呢。
image.png
conda activate py2.7 #激活2.7环境
project=~/rnaseq #项目文件夹
#step5：count
mkdir $project/count
output_dir=$project/count
GTF=/home/cloudam/reference/subread_index/hg38.ncbiRefSeq.gtf
featureCounts -a $GTF -o $output_dir/counts2.txt $project/align/*.bam
用ncbiRefSeq.gtf做注释文件，就能得到熟悉的gene名啦。
image.png
2.双端数据
待续。。。。
网友评论