鹼基資料處理中的演算法研究
最近搞了一些對生物醫學資料的處理，實際的處理邏輯是相當簡單的，其實處理起來還是處理量大，時間效率的問題
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. 一個非常簡單的小程式，拿到字串以後遍歷這上面所有的字元，然後放到一個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給甲方部署上就大功告成啦。