估计这个坑没有太多人能踩到。
前两日,有 GSAman 用户校正完基因结构注释后,导出 gff3 文件,用于 htseq-count 会失败。大体报错如下:
很明显,确实是报错。其中输入的是 .gtf 文件。需要明确的是无论是 GSAman 还是 TBtools,只会输出 GFF3 格式,从信息量来说,GSAman 输出的和 TBtools 输出的是等价的。
按照“常见的软件遇到问题首先找到不可信的软件作者原则”,绝大多数人会认为 gffread 没有问题,所以有问题就一定会优先认为是 GSAman 的问题。这个也正常,你不习惯也得习惯。情况就是这么个情况。
但是在我看来 GSAman 输出是肯定没问题的,信息也是全面的。如果不全面,那就用 TBtools 的 GXF Fix 处理一下。但是对方处理了,还是不行。
?到底问题出现在哪里。无奈之下,也只能要来文件,看看是否可以重现云云。
于是得到结果如下:
显而易见,从软件实现或者设计的角度来说 gffreads 是有问题的,同样的信息输入,如果行顺序有变化,那么就不能良好兼容。这块在 TBtools 里面是不存在的。因为对于 TBtools 来说,输入文件顺序是乱的,这个可太正好了,用户本来就可能手动去修改一些行列顺序,补充或者删除。但是 gffread 可没管这么多。
前述为了让排序更好看,TBtools GXF Fix 输出的行信息,或先输出 exon / CDS ,然后还是 mRNA / gene。我自认为这个很OK,比如 IGV 或者 stringtie 这些软件都没事。就 htseq-count 有事?
但事实就是这样。
索性,还是花了一点时间,对 TBtools GXF Fix 功能输出文件的排序逻辑重新构筑,确保输出行顺序可以更好符合比如 gffread 这些软件的需要,同时也是符合 htseq-count 这些软件的需要。
写在最好
很多时候,感觉写个软件,似乎总是需要“自证清白”?没什么意思,在 GSAman 上是这样,在 TBtools 上更是如此,正如以前。只是我现在想想觉得就这样吧,that's it。
为什么生信分析阻碍重重?因为很少人去care鲁棒性,“能干就行”便是原罪。
加油!想写就写,不想就别写。