1

2

3

samtools sort -n yourfile.bam > yourfile_name.bam

htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/reference/hisat2_reference/Homo_sapiens.GRCh38.86.chr_patch_hapl_scaff.gtf > counts.txt

1

2

3

4

5

6

7

8

9

10

11

# 命令参数

-f | --format default: sam 设置输入文件的格式，该参数的值可以是sam或bam。

-r | --order default: name 设置sam或bam文件的排序方式，该参数的值可以是name或pos。前者表示按read名进行排序，后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序，当输入sam/bam文件是按pos方式排序的时候，两端reads的比对结果在sam/bam文件中一般不是紧邻的两行，程序会将reads对的第一个比对结果放入内存，直到读取到另一端read的比对结果。因此，选择pos可能会导致程序使用较多的内存，它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上，也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的，但HTSeq推荐使用name排序，且一般比对软件的默认输出结果也是按name进行排序的。

-s | --stranded default: yes 设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序；若是单端测序，yes表示read比对到了基因的正义链上；若是双末端测序，yes表示read1比对到了基因正义链上，read2比对到基因负义链上；reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解，一般情况下双末端链特异性测序，该参数的值应该选择reverse（本人暂时没有测试该参数）。

-a | --a default: 10 忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。

-t | --type default: exon 程序会对该指定的feature（gtf/gff文件第三列）进行表达量计算，而gtf/gff文件中其它的feature都会被忽略。

-i | --idattr default: gene_id 设置feature ID是由gtf/gff文件第9列那个标签决定的；若gtf/gff文件多行具有相同的feature ID，则它们来自同一个feature，程序会计算这些features的表达量之和赋给相应的feature ID。

-m | --mode default: union 设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知，对于原核生物，推荐使用intersection-strict模式；对于真核生物，推荐使用union模式。

-o | --samout 输出一个sam文件，该sam文件的比对结果中多了一个XF标签，表示该read比对到了某个feature上。

-q | --quiet 不输出程序运行的状态信息和警告信息。

-h | --help 输出帮助信息。

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

head counts.txt

ENSG00000000003 0

ENSG00000000005 0

ENSG00000000419 1171

ENSG00000000457 563

ENSG00000000460 703

ENSG00000000938 0

ENSG00000000971 1

ENSG00000001036 925

ENSG00000001084 1468

ENSG00000001167 2997

tail count.txt

ENSG00000283696 18

ENSG00000283697 0

ENSG00000283698 1

ENSG00000283699 0

ENSG00000283700 0

__no_feature    3469791

__ambiguous 630717

__too_low_aQual 1346501

__not_aligned   520623

__alignment_not_unique  2849422

1

2

gfold count -ann hg19Ref.gtf -tag sample1.sam -o sample1.read_cnt

gfold count -ann hg19Ref.gtf -tag sample2.sam -o sample2.read_cnt

1

2

3

4

5

6

7

8

9

10

11

12

13

14

GeneSymbol:

For BED file, this is the 4'th column. For GPF file, this is the first column. For GTF format, this corresponds to 'gene_id' if it exists, 'NA' otherwise.

GeneName:

For BED file, this is always 'NA'. For GPF file, this is the 12'th column. For GTF format, this corresponds to 'gene_name' if it exists, 'NA' otherwise.

Read Count:

The number of reads mapped to this gene.

Gene exon length:

The length sum of all the exons of this gene.

RPKM:（#这里需要注意但是双端测序技术还未普及，这里未使用FPKM，况且RPKM和FPKM也不是能很好的代表基因表达水平 ）

The expression level of this gene in RPKM.

1

2

3

4

5

6

7

8

9

10

11

head example.read_cnt

ENSG00000000003 TSPAN6  0   4535    0

ENSG00000000005 TNMD    0   1610    0

ENSG00000000419 DPM1    1588    1207    27.4411

ENSG00000000457 SCYL3   1344    6883    4.07267

ENSG00000000460 C1orf112    1334    5967    4.66292

ENSG00000000938 FGR 0   3474    0

ENSG00000000971 CFH 2   8145    0.0051215

ENSG00000001036 FUCA2   1427    2793    10.6564

ENSG00000001084 GCLC    2462    8463    6.06767

ENSG00000001167 NFYA    5123    3811    28.0378

1

2

wget https://github.com/arq5x/bedtools2/releases/download/v2.26.0/bedtools-2.26.0.tar.gz

tar zxvf bedtools-2.26.0.tar.gz

1

bedtools multicov -bams 1.bam 2.bam 3.bam 4.bam-bed file.bed > read.count.txt

1

2

3

4

5

# 注意，此处的bed文件需要自己处理得到的，需要四列，第一列为chrN，第二列第三列为基因位置，第四列为基因名。类似于：

chr1 0   10000   ivl1

chr1 10000   20000   ivl2

chr1 20000   30000   ivl3

chr1 30000   40000   ivl4