read计数

在开始计数前我们需要获得注释文件，人和小鼠的都下载吧。

wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M25/gencode.vM25.annotation.gtf.gz
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_34/gencode.v34.annotation.gtf.gz

然后就是正式的分析流程了

conda install htseq #先安装软件
mkdir -p RNA-seq/matrix

for ((i=59;i<=62;i++));do htseq-count -r pos -f bam ~/data/SRR35899${i}_sorted.bam gencode.vM25.annotation.gtf > ~/data/RNA-seq/matrix/SRR35899${i}.count;done

回过头再看看htseq-count怎么用吧，引用自hoptop

-f bam/sam：指定输入文件格式，默认SAM
-r name/pos: 你需要利用samtool sort对数据根据read name或者位置进行排序，默认是name
-s yes/no/reverse: 数据是否来自于strand-specific assay。DNA是双链的，所以需要判断到底来自于哪条链。如果选择了no，那么每一条read都会跟正义链和反义链进行比较。默认的yes对于双端测序表示第一个read都在同一个链上，第二个read则在另一条链上。
-a 最低质量，剔除低于阈值的read
-m 模式 union（默认）, intersection-strict and intersection-nonempty。一般而言就用默认的，作者也是这样认为的。
-i id attribute: 在GTF文件的最后一栏里，会有这个基因的多个命名方式（如下）， RNA-Seq数据分析常用的是gene_id，当然你可以写一个脚本替换成其他命名方式。

这一步大概花了6h，得到4个矩阵文件，得到4个.count文件，把它们拷贝到电脑上，到这里Linux的部分就结束了，我们转战R。