转录组数据分析笔记（6）——HTseq计数定量

HTseq也是对有参考基因组转录数据进行表达量分析的，主要用于reads计数。这个软件功能就比较专一，不像stringtie还需要运行prepDE.py脚本进行数据转化，直接一步到位。那为什么我一开始不用HTseq呢？~~因为我遇到一个bug~~ 主要还是运算速度的问题，我比较了两种定量方式，HTseq定量虽然只有一步，但是速度远不如stringtie，也可能是我的问题，下面会说到。

1. HTseq定量获得raw count

vim一个新脚本，输入如下命令：

#!/bin/bash
for i in `seq 194 209`
do 
    htseq-count -f bam -s no \
    /media/sf_/data/fastq/bam/ERR1698"$i".bam \     # 输入bam文件
    /media/sf_/data/ref/Arabidopis_thaliana.gtf     # 参考基因组注释文件
    > /media/sf_/data/fastq/count/ERR1698"$i".count     # 自定义输出文件
done

参数详解

-f # 设置输入文件格式，可以是bam或者sam
-s # 设置是否是链特异性测序，设置no每一条reads都会和正义链和反义链进行比较

保存运行以后发现这个程序只能分配一个线程（也可能是我没找到分线程的方法），所以可以根据电脑内核数分几个批处理一起运行会快很多（不然就等着干瞪眼= =）。

还有一点非常重要！bam文件需要提前按照名称排序，不然会出现绝大部分reads mapping不到参考基因组，这种情况会在屏幕上输出提示信息，但是程序还是会继续跑……这时候就别犹豫了赶紧kill这个程序，就算跑完了数据都不能用。可以用samtools sort -n对bam文件进行名称排序，但是排序之后无法再用samtools index建立索引文件，这会导致HTseq运行速度比蜗牛还慢。暂时没找到更好的办法摊手。

经过漫长长长长长的时间等待，我们可以看看结果文件的head和tail（这里就放一张图吧）：

前面记录了基因名称和mapping上的reads数，最后5行对应不同的mapping情况，在不同的模式下意义不同，官网给出的区别如下图，默认是union模式：

计数结果也可以用multiqc合并，生成在线报告，这里可以直观地看到每个样品比对上的reads数百分比，这里16个样品的比对率都超过80%，说明计数结果都还不错。

2. HTseq结果文件处理

HTseq计数定量后得到的是每一个样品的每个基因reads数，我们需要合并每个样品定量数据，手动修改成DESeq2能识别的raw count表达矩阵，还需要再准备一个样本列表矩阵，才能进行后续的DESeq分析。参考一下stringtie最后生成的表达量矩阵文件，我们也需要将HTseq定量结果整理成csv格式（逗号作为分隔符），第一列是基因名，后面是按照样品序列的排序，中间是表达矩阵。

再来看一看HTseq定量生成的文件详情，同样第一列是基因名，后面是raw count数量，^I 表示两列数据是以制表符tab键分隔的，$为换行符。

我的方法比较笨比，除第一个ERR1698194.count文件保留外，其他所有count文件第一列删去并命名为cut.count，然后合并ERR1698194.count和其他所有cut.count文件，再将所有的制表符替换为逗号，最后加上第一行行名和改文件名。

用awk命令删除第一列，写入到新的cut.count文件中：

for i in `seq 195 209`
do 
    cat ERR1698"$i".count | awk '{$1 = ""; print $0}' > ERR1698"$i"cut.count
done

paste组合ERR1698194样品和其他cut.count文件到alldata.count:

$ paste ERR1698194.count *cut.count > alldata.count

看看alldata.count的数据格式，列数没有问题，但是awk删除列产生了空格：

用sed命令删除所有空格，替换所有制表符为逗号（两步可以合一步）：

$ sed 's/ //g' alldata.count > alldata1.count

$ sed 's/\t/,/g' alldata1.count > alldata2.count

这样就手动生成符合csv格式的文件了，只需加上第一行：

这里样本量比较少，我直接vim复制粘贴的方法加了第一行，重命名一下文件就完成了表达矩阵的制作，可以用于DESeq2分析了！

因为本人比较小白，上面处理过程就有些啰嗦了，总的思路就是改成csv格式文件的样式就可以。

样本列表矩阵的制作过程和stringtie一模一样，点击这里查看，本篇不再赘述。