提取最长转录本的代码

来源：互联网发布：软件产品设计文档模板编辑：程序博客网时间：2024/05/17 03:05

1.最长转录本：就是gene的id相同，但是序列的长度不一样，应该挑选出序列最长作为后面的分析

下面就是一个转录本的id文件，都是来自同一个gene,但是转录本的id和长度均不相同

分析：发现除了第4列表示的基因id相同以外，其他的列均不相同。

2.过滤最长转录本的编程思路：

考虑到绵羊的转录本的id的复杂性，主要按一下几个步骤来实现：

step1: 读取转录本序列的fasta文件，并统计序列的长度，生成一个ID文件，该ID文件比原来的转录本id多出一列，即变成13列，而第13列内容为序列的长度

>ENSOART00000000006 mt_genbank_import:known chromosome:Oar_v3.1:MT:2745:3699:1 gene:ENSOARG00000000006 gene_biotype:protein_coding transcript_biotype:protein_coding 955

step2: 根据上面的ID文件的第4列信息(gene:XXXX)和最后一列的长度信息，过滤得到最长的转录本的ID_filter文件

step3: 通过awk 将ID_filter的最后一列过滤掉，生成最终的ID_final因为如果不过滤掉就和最初输入的转录本序列对应不上

step4:通过存一个哈希，将ID_final中的id序列挑选出来

然后考虑用一个*.sh脚本将其全部串起来

问题：

能够将三个程序写成一个程序，而实现上面4步的功能？

除了运用perl的语言存哈希以外，还有其他的方法吗？

后续会继续更新

0 0