碱基数据处理中的算法研究

本文讲述了作者在处理大量生物医学数据,特别是DNA碱基序列互补转换过程中,从原始的线性处理方式优化到流式读写、内存分块的高效算法,大幅提升了处理速度,最终实现大文件在合理时间内完成处理。

最近搞了一些对生物医学数据的处理,实际的处理逻辑是相当简单的,其实处理起来还是处理量大,时间效率的问题

1.需求

医学数据中会有碱基对序列数据,他们的要求是碱基序列转换为和它互补的碱基序列,将其做成一个网页端的,可以直接上传文件处理大文件数据,也可以直接在网页中贴上碱基序列处理完展示出来。

2.需求分析

当了解了一些基本的生物碱基序列的基本知识后,他们就是处理一些DNA碱基序列,其中的规则就是A-T,C-G,也就是想数据中的A转化成T,T转化成A,C转化成G,G转化成C,然后将处理完的数据进行逆序,就可以得到互补的碱基序列了。整体的逻辑还是比较简单的,困难的点还是当处理大一点的数据量的时候能不能快速的处理好。

3.算法设计经历

经过分析需求,感觉还是比较简单的,快速的搭建了一个网页端的页面,技术选型也是用了最原生的servlet进行操作,这些简单的操作,没必要用到框架之类的。

整个思路就是对应转换字符,然后将字符串逆序的过程。很快第一版代码就出来了。

 public String handle(String text){
        StringBuffer sb = new StringBuffer();

        for (int i = 0; i < text.length(); i++) {
            if (text.charAt(i)=='A'){
                sb.append("T");
            }else if (text.charAt(i)=='C'){
                sb.append("G");
            }else if(text.charAt(i)=='G'){
                sb.append("C");
            }else if (text.charAt(i)=='T'){
                sb.append("A");
            }
        }
        return sb.reverse().toString();
    }

一个非常简单的小程序,拿到字符串以后遍历这上面所有的字符,然后放到一个StringBuffer中去,然后逆序。逻辑其实还是非常简单的,我测试了一下确实逻辑都是正确的,我以为大功告成的了,就开始准备写一个上传文件的代码,处理一个文件上的数据。

这是我才发现文件中的数据是一行一行的并且文件的第一行有个基因标识行,效果如下:
在这里插入图片描述

再次询问需求后,他们要求将第一行的标识信息不变,每一行的数据还是在同一行上,然后将其做成互补的碱基序列

其实这个也很好处理,我只要判断出第一行来,把他单独拿出来,然后将后面的一行一行数据处理成互补的碱基序列,我最后将每一行的数据倒序用\r\n拼接输出来,最后在头部放上标识信息不就可以啦,看样子也不是很难处理嘛

    public String handle(String text){
        String[] textSplit = text.split("\n");
        //处理每一行数据
        for (int i = 1; i < textSplit.length; i++) {
            textSplit[i] = matchHandle(textSplit[i]);
        }
        String result = textSplit[textSplit.length-1];
        //倒序拼接起来
        for (int i = textSplit.length-2; i >= 1; i--) {
            result = result +"\r\n"+textSplit[i];
        }
        //拼接上第一行数据
        result = textSplit[0] +"\r\n"+result;
        return result;
    }

    public String matchHandle(String text){
        StringBuffer sb = new StringBuffer();

        for (int i = 0; i < text.length(); i++) {
            if (text.charAt(i)=='A'){
                sb.append("T");
            }else if (text.charAt(i)=='C'){
                sb.append("G");
            }else if(text.charAt(i)=='G'){
                sb.append("C");
            }else if (text.charAt(i)=='T'){
                sb.append("A");
            }
        }
        return sb.reverse().toString();

    }

写完以后我拿着几十行的数据测试了算法的正确性,发现没有什么问题,那就加大文件呗,看看速度怎么样,我就上传了一个14M的碱基序列,然后我整个人都不好了,我发现处理这么小的数据程序整整跑了30分钟,也要是再大一点那不要时间更长了。这第一版代码是有问题的,我赶快的分析一下代码,看看是哪里时间消耗这么大。

分析一顿后发现这里面有消耗时间的就是我需要一行一行的读并加了字符串的拼接操作,并且还有一个逆序操作,这些都是一些比较耗时间的操作,我发现14M中的数据一共要有25万条,所有我一点一点的尝试修改一下。(这里其实字符串的拼接操作JVM底层已经做了相应的优化)

首先我将StringBuffer中的逆序操作修改掉了,不在用StringBuffer处理字符串,改用ArrayList处理,最后在Arraylist中逆序生成字符串

  public String matchHandle(String text){
        List list = new ArrayList();

        for (int i = 0; i < text.length(); i++) {
            if (text.charAt(i)=='A'){
                list.add("T");
            }else if (text.charAt(i)=='C'){
                list.add("G");
            }else if(text.charAt(i)=='G'){
                list.add("C");
            }else if (text.charAt(i)=='T'){
                list.add("A");
            }
        }
        Collections.reverse(list);
        return String.join("",list);
    }

经过这样子处理以后,我测试了一下发现,效果竟然好了很多,现在处理14M的数据需要10分钟就可以,但是时间还是有点长。既然这里是逆序的操作占用了很多的时间为什么不直接不让他逆序了,我可以在处理完转换后,直接从后向前输出出来不就好了,既然采用这用思路,我在处理串的时候就不用在分成切分成数组了,我直接把他当成一整个字符串来处理,将中的回车符也添加处理起来,这样他写进文件的时候,还照样的回车。

就这样我采用了一种空间换上时间的操作,将一个时间复杂度O(n2)的代码转化成了O(n)的代码,代码如下:

  public String handle1(String text){
        StringBuffer sb = new StringBuffer();

        String desc = text.substring(0,text.indexOf("\n"));

        for (int i = text.length()-1; i > desc.length() ; i--) {
            if (i == text.length()-1){
                sb.append(desc+"\n");
            }
            sb.append(mathchCharHandle(text.charAt(i)));
        }
        return sb.toString();
    }

    public char mathchCharHandle(char c){
        if (c=='A'){
            return 'T';
        }else if(c=='T'){
            return 'A';
        }else if (c=='C'){
            return 'G';
        }else if (c=='G'){
            return 'C';
        }else{
            return c;
        }
    }

我采用了这种做法以后,经过测试我发现14M的数据,不用1秒中就可以处理完毕这些数据,这与第一次的算法处理时间效率直接大大提升啊。

经过这几步的优化,我想进一步测试一下的数据,我用了一个3g的数据进行测试,程序直接报错了。

5-Dec-2020 15:46:42.820 SEVERE [http-apr-8080-exec-1] org.apache.catalina.core.StandardWrapperValve.invoke Servlet.service() for servlet [com.test.web.servlet.UploadHandleServlet] in context with path [/base_pair_match] threw exception
 java.lang.NegativeArraySizeException
	at org.apache.commons.fileupload.disk.DiskFileItem.get(DiskFileItem.java:329)
	at org.apache.commons.fileupload.disk.DiskFileItem.getString(DiskFileItem.java:379)
	at com.test.web.servlet.UploadHandleServlet.doPost(UploadHandleServlet.java:72)
	at javax.servlet.http.HttpServlet.service(HttpServlet.java:648)
	at javax.servlet.http.HttpServlet.service(HttpServlet.java:729)
	at org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:292)
	at org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:207)
	at org.apache.tomcat.websocket.server.WsFilter.doFilter(WsFilter.java:52)
	at org.apache.catalina.core.ApplicationFilterChain.internalDoFilter(ApplicationFilterChain.java:240)
	at org.apache.catalina.core.ApplicationFilterChain.doFilter(ApplicationFilterChain.java:207)
	at org.apache.catalina.core.StandardWrapperValve.invoke(StandardWrapperValve.java:212)
	at org.apache.catalina.core.StandardContextValve.invoke(StandardContextValve.java:106)

然后我尝试的用sublime text打开它,我的机器直接内存爆满,打不开。所以,报错应该是文件中的数据不能直接读到内存中去,所以对于这么大的数据我应该采用读一部分就处理一部分,并直接写出到文件中去,采用内外存交换这样子,才可以。这里处理文件的数据的逻辑还是不变的,刚刚做的写的算法也可以用上了。

现在的问题就是如何将大文件进行切分并处理好一些数据。

大文件就牵扯到文件内容的分块与合并,我一共想到了两个解决方案,第一种就是我从最后一行开始读数据将,一次读一块(块的大小自己设定),然后在从尾部插入的新的文件,我经过尝试后发现这个倒序的读文件处理的操作很难控制一块一块的数据处理,所以我就尝试了我的第二种方案,就是采用从前往后读文件,然后分块处理,然后写到新文件中是头部开始插入,我很快尝试做了一个按行的分块逻辑的程序的测试代码。

  public static void main(String[] args) throws IOException {
        long t1 = System.currentTimeMillis();
        String path = "D://Caenorhabditis_elegans.dna.chromosome.fa.txt";
        String newPath = "D://qq.txt";
        RandomAccessFile br = new RandomAccessFile(path, "r");
        String str = null;
        StringBuilder app = new StringBuilder();
        long i = 0;
        String readText = null;
        long lineCount = 0;
        //一块多少行
        long subCount = 250000; //Caen   10000 30188ms   100 52795ms 100000  30295ms
        //读行数
        File file = new File(path);
        long fileLength = file.length();
        LineNumberReader lnf = new LineNumberReader(new FileReader(path));
        lnf.skip(fileLength);
        lineCount = lnf.getLineNumber();
        //块的数量
        long blockNum = lineCount / subCount;

        long blockCount = 0;

        String firstLine = null;
        while ((str = br.readLine()) != null) {

            if (i == 0) {
                //处理第一块
                firstLine = str;
                i++;
            } else {

                i++;

                if (i % subCount == 0 || i == lineCount+1) { //处理每一块中的内容或者是不够一块的内容
                    app.append(str + "\n\r");
                    
                    readText = handle2(app.toString());

                    readText = readText.substring(2, readText.length());
                    appendFileHeader((readText + "\r\n").getBytes(), newPath);
                    app.delete(0, app.length());
                    blockCount++;
                } else {
                    app.append(str + "\n\r");
                }
            }
        }
        appendFileHeader((firstLine + "\r\n").getBytes(), newPath);
        long t2 = System.currentTimeMillis();
        System.out.println(t2 - t1 + "ms");
        br.close();

    }
    
    public static String handle2(String text) {

        StringBuffer sb = new StringBuffer();

        for (int i = text.length() - 1; i >= 0; i--) {
            sb.append(mathchCharHandle(text.charAt(i)));
        }
        return sb.toString();
    }

    public static char mathchCharHandle(char c) {
        if (c == 'A') {
            return 'T';
        } else if (c == 'T') {
            return 'A';
        } else if (c == 'C') {
            return 'G';
        } else if (c == 'G') {
            return 'C';
        } else {
            return c;
        }
    }

    /**
     * 向src文件添加header
     *
     * @param header
     * @param srcPath
     * @throws Exception
     */
    private static void appendFileHeader(byte[] header, String srcPath) {
        try {
            RandomAccessFile src = new RandomAccessFile(srcPath, "rw");
            int srcLength = (int) src.length();
            byte[] buff = new byte[srcLength];
            src.read(buff, 0, srcLength);
            src.seek(0);
            src.write(header);
            src.seek(header.length);
            src.write(buff);
            src.close();
        } catch (Exception e) {
            e.printStackTrace();
        }

    }

这个算法中难点就在于他在处理最后剩下的一部分,也就是分块后,不够一块的部分,我们该判断出来并进行处理,我采用的当他处理完每一块以后,我要将剩下的所有行也一并放进一个StringBuffer中,直接进行处理。实际逻辑与处理块中的数据是一样的,所以只需多加一个判断就可以了

然后进行测试发现确实代码不会报错了,不过我拿14M的数据测试了一下,发现处理完以后要30秒才可以完成,要是按这个时间的话,处理3G的数据估计要很久很久,这是不合适的啊。其实我们以前处理14M的文件的时候,不用1秒就可以的啊,现在怎么慢了这么多,应该是我分块的算法搞的不合适。我发现我不应该按行的逻辑去分块,这里面要读取每一行的操作实际上的非常占用时间的,所以我就直接改用流的方式进行读入处理。

这个代码同样也会遇到最后不到一块的情况,他的这个处理方式相对于刚刚行的处理可能要复杂一点,因为我要将每一块的数据读到一个Buffer中,这个buffer的大小对于最后一块要控制好,而且每一块中读多少个字节也是需要好好确认好的,要不然会读到某一行的一半,数据在逆序拼接起来的时候就会乱了,所以我确认了我处理的数据一行一共有61个字节,所以我这个就以61个字节的倍数作为buffer的大小。

 long t1 = System.currentTimeMillis();
        String path = "D://GRCh38_latest_genomic.fna";

        String   newPath = path.substring(0,path.lastIndexOf(".")) + "-" +new Date().getTime()+".txt";

        RandomAccessFile br = new RandomAccessFile(path, "r");
        
        StringBuilder app = new StringBuilder();
        long i = 0;
		//每一块的大小
        int block = 61 * 2000000;

        
        File file = new File(path);
        long fileLength = file.length(); //字节
		//块数
        long blockNum = fileLength / block;

        int resLength = (int) (fileLength - blockNum * block);
        //不能整除就要块数+1
        if (fileLength % block ==0){
            blockNum = blockNum;
        }else {
            blockNum++;
        }

        String firstLine = br.readLine();

        byte [] bytes = new byte[block];

        while (i< blockNum){
            if (i == blockNum-1){
                bytes = new byte[resLength];
                br.read(bytes,0,resLength);
            }else{
                br.read(bytes,0,bytes.length);
            }

            String text = new String(bytes).trim();

            text = handle2(text);
            appendFileHeader((text + "\r\n").getBytes(), newPath);
            i++;
        }
        appendFileHeader((firstLine + "\r\n").getBytes(), newPath);
        long t2 = System.currentTimeMillis();
        System.out.println(t2 - t1 + "ms");
        br.close();

    }


    public static String handle2(String text) {

        StringBuffer sb = new StringBuffer();

        for (int i = text.length() - 1; i >= 0; i--) {
            sb.append(mathchCharHandle(text.charAt(i)));
        }

        return sb.toString();
    }

    public static char mathchCharHandle(char c) {
        if (c == 'A') {
            return 'T';
        } else if (c == 'T') {
            return 'A';
        } else if (c == 'C') {
            return 'G';
        } else if (c == 'G') {
            return 'C';
        }else if(c=='\n'){
            return '\r';
        }else if(c=='\r'){
            return '\n';
        } else{
            return c;
        }
    }

    /**
     * 向src文件添加header
     *
     * @param header
     * @param srcPath
     * @throws Exception
     */
    private static void appendFileHeader(byte[] header, String srcPath) {
        try {
            RandomAccessFile src = new RandomAccessFile(srcPath, "rw");
            int srcLength = (int) src.length();
            byte[] buff = new byte[srcLength];
            src.read(buff, 0, srcLength);

            src.seek(0);
            src.write(header);

            src.seek(header.length);
            src.write(buff);

            src.close();
        } catch (Exception e) {
            e.printStackTrace();
        }
    }

经过这样的处理,我测试了一下14M的数据基本不用1秒就就可以解决,那2G左右的数据一分半就可以完成,这完全是在可以接受的范围内的。

4.项目总结

算法设计好了,我马上搭了一个web的页面,将自己分析好的算法放上去集成测试一下,发现没啥大问题,现在的处理满足他们医学数据的处理应该没有问题了,最后达成war给甲方部署上就大功告成啦。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值