DNA pattern的高效检索

赛题说明

grep 是一个非常常用的Linux命令,用它可以搜索文本文件的特定的字符串并将它们打印出来。对于核酸序列文件,比如描述参考基因组的fasta文件,或者描述测序下机数据的fastq文件,我们也希望有类似的软件,可以简便地搜索出目标模式子序列。但是对于核酸序列文件,我们有一些要求,比如它必须能够有细微的容错,比如错一到两个碱基。

 

本赛题要求的dnagrep就是这样的软件,它必须支持这样的命令:

dnagrep [options] <pattern> <input>

输入:

pattern是核酸序列,本软件只考虑处理DNA序列,所以它应该只有A/T/C/G/N,或者小写的a/t/c/g/n

input 是输入的文件,可以是fastq或者fasta类型的文件,如果该参数没给出,则从STDIN输入。

options如下:

-s --sub 允许待搜索序列与pattern相差的替换subsititution,默认为0

-i --indel 允许待搜索序列与pattern相差的插入删除indel,默认为0

-d --distance 允许待搜索序列与pattern相差的替换+插入删除,如果没有给出,则默认为 sub+indel

-t --thread 执行搜索的线程数,默认为1

输出:

结果输出都是到STDOUT,方便重定向,错误提示和警告等其他信息输出到STDERR

如果是fastq的匹配,则输出每一个匹配行<前序列>  <pattern 匹配的序列>   <后序列>

如果是fasta的匹配,则输出每一个匹配位<前10bp>  <pattern 匹配的序列>   <后10bp>

中间都是空格隔开

对于<pattern 匹配的序列>不是完美匹配的,采用以下方式:

替换用红色高亮

插入用橙色高亮

删除的位置使用橫杠 - 代替

 

示例命令

dnagrep -s2 -i1 -d2 AAAATTTTCCCCGGGG hg19.fa

表示从hg19.fa中搜索和AAAATTTTCCCCGGGG相似的序列,这些序列只需要最多不超过2步变换就可以AAAATTTTCCCCGGGG,这2步变换中最多只能有2个替换和1次插入删除

 

注意事项:

1,如果输入的文件是fasta文件则需要忽略其中的换行符,因为很多参考基因组fasta文件是每一行限宽的(比如70个字符)

2,需要支持查找反向互补的序列,即反向互补也认为是匹配

3,结合实际应用,同时也是为了简化问题,sub 的值限定为不超过2,indel的值限定为不超1,distance的值限定为不超过3。如果不是,则报错

4,序列比对都忽略大小写

5,pattern 序列中可以有N的存在,它可以匹配任意碱基。但是输入文件中的N只能匹配N本身。

6,如果是STDIN输入,则需要自动判断是fastq还是fasta,文件输入的可以通过扩展名进行判断

 

数据

fasta 文件:从UCSC 网站上下载http://hgdownload.cse.ucsc.edu/downloads.html

fastq 文件:OpenGene网站上提供了一对pair-end的数据供下载,见http://opengene.org/dataset.html 我们将使用其中的R1文件http://opengene.org/data/R1.fq.gz 的数据

 

提交结果

1、 源代码 (如果是C/C++源码目录下必须有Makefile文件,使用标准的make命令可以在程序目录下生成dnagrep可执行文件,如果是其他语言,请在文档中指明使用方式,并遵照同样的参数格式)

2、 技术文档(包括使用帮助文档,测试文档等)

评估方法

hg19.fasta文件,和已解压的R1.fq文件,使用一系列的pattern串和不同的参数进行dnagrep命令操作,并使用以下指标进行评分。

评估指标

1、 方法说明(5’)。考查方法的完整性和可用性。 

2、 环境的实用性或经济性(20’)。根据环境说明,配置相应的环境,评价环境可用性并打分,环境实用性和经济性很差的情况下,评定小组可放弃搭建。 

3、单线程执行的速度(50’)。 即使用 -t1 的执行速度,此项考查主程序的优化

4. 四线程执行的速度(20’)。  即使用 -t4 的执行速度,此项考查多线程的优化

5. 其他(5’) 

6. 加分项,其他有助于读写该压缩合适的附加功能给予适当加分,加分总分不能超过10分