已經(jīng)在多個(gè)地方看到有人吐槽htseq-count了我纪，原因就是它對(duì)bam的排序有嚴(yán)格要求慎宾，必須是按照reads name來排序才能運(yùn)行丐吓。本著”實(shí)踐是檢驗(yàn)真理的唯一標(biāo)準(zhǔn)“的理念，我自己又試了一下趟据。
#SRR3286802.s.bam已經(jīng)按坐標(biāo)排好序了
htseq-count -s no -r pos -f bam /map/SRR3286802.s.bam /ref/gene27655.gff \
> /count/SRR3286802.count 2> /count/SRR3286802.err
結(jié)果真的報(bào)錯(cuò)了
Error occured when processing SAM input (record #30072152 in file /ifs1/Grp3/huangsiyuan/learn_rnaseq/mRNA-seq/map/SRR3286802.s.bam):
Maximum alignment buffer size exceeded while pairing SAM alignments.
[Exception type: ValueError, raised in __init__.py:686]
我在SEQanswers上找到了一個(gè)解答券犁，提供了兩種解決方案，我都試了一遍汹碱。
'
Two options: （換成reads name排序方式粘衬；換軟件）
You can either re-sort your files lexographically/alphabetically and run htseq-counts with the lexigraphic order (in other words, NOT --order=pos)
or
You can switch to the much faster and generally more user friendly (easier options, easier output file format) program "featureCounts". Part of the "subread" package. [http://subread.sourceforge.net/](http://subread.sourceforge.net/)
I would switch to subread. htseq-counts is incredibly slow and very poor at handling large files without crashing.
'
1. 重新排序
$ cat 2.sh
ls SRR328680{2..7}.bam | while read id
do
samtools sort -n ${id} -o ${id%.*}.namesort.bam -@ 2 &
done
$ cat 3.sh
ls *.namesort.bam | while read id
do
htseq-count -s no -r name -f bam ${id} /ref/gene27655.gff > /count/${id%%.*}.count 2> /count/${id%%.*}.err &
done
理論上是可行的，以前試過咳促。但這次又出現(xiàn)了新問題稚新，跑出來的結(jié)果都是空文件。
$ cat SRR3286807.count
__no_feature 45205743
__ambiguous 0
__too_low_aQual 0
__not_aligned 1336845
__alignment_not_unique 2045753
這時(shí)我猜測跪腹，htseq-count要求bam(chr1,chr2...)和gff(1,2...)中的染色體名稱相同褂删，不然無法判斷，更改gff的chr名稱后冲茸，再按同樣的命令運(yùn)行一次屯阀。
awk '{ print "chr"$0 }' gene27655.gff > gene27655_re_name_chr.gff
結(jié)果還是一樣空的?! —— (運(yùn)行過程中提示W(wǎng)arning: No features of type 'exon' found.)我只能猜測是自己前面提取gene區(qū)間至gff（為了直接統(tǒng)計(jì)gene區(qū)間內(nèi)的reads數(shù)）這一步有些多此一舉了。
于是又拿完整的gff文件試了下轴术，Arabidopsis_thaliana.TAIR10.42.gff3难衰，又說某行第九列注釋信息中沒有gene_id信息，看了下確實(shí)沒有膳音。
1 araport11 exon 3631 3913 . + . \
Parent=transcript:AT1G01010.1;\
Name=AT1G01010.1.exon1;\
constitutive=1;\
ensembl_end_phase=1;\
ensembl_phase=-1;\
exon_id=AT1G01010.1.exon1;\
rank=1
又去NCBI下載了一個(gè)gtf注釋文件（基因組版本是一樣的）, GCF_000001735.4_TAIR10.1_genomic.gtf.gz召衔，這回終于有g(shù)ene_id了。更改染色體名稱后祭陷，重新用htseq-count計(jì)數(shù)苍凛。
總結(jié)一下
通過上面的折騰，對(duì)htseq-count計(jì)數(shù)有了新的認(rèn)識(shí)兵志。
- bam中的reads根據(jù)名稱排序醇蝴；
- 基因組feature都保留，不要多此一舉專門提出feature為gene的記錄想罕；
- gff/gtf的第9列注釋要有g(shù)ene_id悠栓；
- bam和gff/gtf的染色體名稱保持一致。
其中按价，第2惭适、3點(diǎn)可以自己另外設(shè)置，見htseq-count計(jì)數(shù)
2. 更換軟件
nohup wget https://nchc.dl.sourceforge.net/project/subread/subread-1.6.3/subread-1.6.3-source.tar.gz &
tar -zxvf subread-1.6.3-source.tar.gz
cd src/
make -f Makefile.Linux #編譯成功后楼镐，可執(zhí)行文件存儲(chǔ)在../bin/中
Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ...
nohup ~/mysoft/subread-1.6.3-source/bin/featureCounts -F GTF -t gene -g gene_id -T 4 \
-a ~/learn_rnaseq/mRNA-seq/ref/gene27655.gff -o ~/learn_rnaseq/mRNA-seq/count/test \
~/learn_rnaseq/mRNA-seq/map/SRR328680*.s.bam &
上面這一步運(yùn)行完之后會(huì)多出test.summary癞志，test這兩個(gè)文件
nohup ~/mysoft/subread-1.6.3-source/bin/featureCounts -F GTF -t gene -g gene_id -T 4 -p \
-a ~/learn_rnaseq/mRNA-seq/ref/gene27655.gff -o ~/learn_rnaseq/mRNA-seq/count/all \
~/learn_rnaseq/mRNA-seq/map/SRR328680*.s.bam &
上面這行只多了-p選項(xiàng)，為了看看在雙端測序中統(tǒng)計(jì)fragments和reads有什么區(qū)別框产。運(yùn)行完了多出all.summary凄杯，all兩個(gè)文件错洁。
從all和text兩個(gè)矩陣文件中，看不出什么明顯的差別戒突，但是加-p之后屯碴，運(yùn)行時(shí)間長了很多。但不管怎么樣膊存，比htseq-count快多了导而。
