de novo的转录组数据，比对的时候一般用的是自己组装好的trinity.fasta序列(挑选最长蛋白的转录本序列)来做参考，用bowtie2等工具直接将原始序列比对即可。所以比对 sam/bam文件本身就包含了参考序列的每一条转录本序列ID，直接对 sam/bam文件进行counts就知道每一个基因的表达量啦！

本来我是准备自己写脚本对sam文件进行counts就好，但是发现了samtools自带这样的工具：http://www.htslib.org/doc/samtools.html 

如果是针对基因组序列，那么这个功能用处不大，但是针对转录本序列，统计出来的就是我们想要的转录本表达量。

samtools idxstats tmp.bowtie2.sorted.bam |head
TR3|c0_g1_i1 1276 418 0
TR6|c0_g1_i1 1271 10 0
TR6|c0_g1_i2 944 5 0
TR6|c0_g1_i3 1281 4 0
TR6|c0_g1_i4 1224 53 0
TR6|c0_g1_i5 855 16 0
TR19|c0_g1_i2 1428 19 0
TR19|c0_g1_i3 2536 624 0
TR19|c0_g1_i4 3072 105 0
TR19|c0_g1_i5 1685 0 0

软件官网说明书，说的很清楚：

samtools idxstats in.sam|in.bam|in.cram

Retrieve and print stats in the index file corresponding to the input file. Before calling idxstats, the input BAM file must be indexed by samtools index.

The output is TAB-delimited with each line consisting of reference sequence name, sequence length, # mapped reads and # unmapped reads. It is written to stdout.

第三列，就是我们想要的表达量数据啦，比对到每个转录本序列的reads数量。

大家从我的转录本序列ID上面如果可以看出些什么问题，欢迎跟我交流，直接给我email就好了，jmzeng1314@163.com 

现在知道了每个转录本的表达量，把每个样本都做一下，就知道表达矩阵了，做差异分析就很简单了。但是得到的是差异转录本列表，不明白这些ID背后的意义，需要取注释，才能做下一步分析。

ls *sorted.bam |while read id
do
echo $id ${id%%.*}.t.counts
nohup samtools idxstats $id 1>${id%%.*}.t.counts 2>/dev/null  &
done