以前是没有想过用这个软件的，直到有一个我的htseq无法对比对的bam文件进行基因计数（后来我才发现htseq无法计数的原因是gtf版本不同导致坐标不同，而且gtf对染色体编号没有加上chr），我简单搜索了一下，发现bedtools multicov也有类似的功能，所以我选择它来试试看！
首先注意它需要sort的bam文件及bam的index
bedtools multicov depends upon index BAM files in order to count the number of overlaps in each BAM file. As such, each BAM file should be position sorted (samtool sort aln.bam aln.sort) and indexed (samtools index aln.sort.bam) with either samtools or bamtools.
首先安装它：
wget https://github.com/arq5x/bedtools2/releases/download/v2.23.0/bedtools-2.23.0.tar.gz
解压开后
Make clean
Make all
然后就可以看到它的bin目录下全部是程序啦
命令很简单的
bedtools multicov [OPTIONS] -bams BAM1 BAM2 BAM3 ... BAMn -bed <BED/GFF/VCF>
By default, multicov will report the count of alignments in each input BAM file that overlap.
例子：
bedtools multicov -bams aln1.bam aln2.bam aln3.bam -bed ivls-of-interest.bed
ivls-of-interest.bed这个文件是必须的，可能需要自己制作，其实用gtf文件也可以的
chr1 0 10000 ivl1
chr1 10000 20000 ivl2
chr1 20000 30000 ivl3
chr1 30000 40000 ivl4
输出结果前三列是坐标，第四列是基因名，跟我们的bed文件一样，只是最后三列是三个样本的计数，是添加上来的！
chr1 0 10000 ivl1 100 2234 0
chr1 10000 20000 ivl2 123 3245 1000
chr1 20000 30000 ivl3 213 2332 2034
chr1 30000 40000 ivl4 335 7654 0
同样是对gene的reads计数，bedtools的multicov工具与htseq-count的区别是
i'd guess it's due to the fact that htseq-count only reports one hit per aligned read assuming that read is aligned uniquely and does not overlap multiple features listed in your GTF. if an aligned read hits more than one feature in your GTF then it doesn't report that hit. bedtools gives you raw hits which includes every 1 hit for every intersection of every alignment with any features in the GTF no matter how many times it aligned or how many features it hit. you might think, "wow, htseq-count is dropping a lot of information". yes, it is! i've moved to using other tools to count hits to genes (RSEM/eXpress) since they disambiguate those ambiguous alignments and as a result you get counts for all of your aligned reads. in a genome with alternative splicing you lose too much data using htseq-count, in my opinion.
而且专门有个文献在讨论这个问题
http://barcwiki.wi.mit.edu/wiki/SOPs/rna-seq-diff-expressions
http://www.nature.com/nbt/journal/v31/n1/abs/nbt.2450.html
Differential analysis of gene regulation at transcript resolution with RNA-seq
下面我讲一个实际的例子
我的bam文件如下
bedtools multicov -bams 740WT1.bam 741WT2.bam 742KO1.bam 743KO2.bam -bed mm9.bed
得到的这个矩阵就可以去用DESeq包来进行差异分析啦！
