Подсчет соответствующих записей в большом файле биоинформатики

У меня есть рабочий пример фрагмента кода, который открывает файл, собирает информацию о содержимом и выводит карту, содержащую информацию.

Файл

Тип файла - это собственное создание, называемое Xsam-файлом. Для тех, кто заинтересован, он основан на файле sam, который обычно используется в биоинформатике. Каждый файл начинается с раздела заголовка, из которого каждая строка начинается с «@» и может быть безопасно проигнорирована этим -> в заголовке обычно не более 1000 строк. Остальная часть файла состоит из пар считывания. Каждое чтение занимает одну строку, и строки всегда попарно. Пример пары чтения:

D43P8DQ1:194:H3W7GADXY:1:2104:5516:41310    99  mm01_24611438_24616266  2276    150 5S41M   =   2360    133 NNAGGTGAATAGAATTATACCATATCGTAGTCCTTTTTGTACAATA  ~~HHHFHHBGIJDHIFHGGGIIIJJGICGGCBHIIJJHIIEGHCGF  xl:i:2276   xr:i:2316   xs:i:41 xd:A:f  xm:A:u  xa:A:""     xL:i:2276   xR:i:2402   xS:i:127    xW:i:43 xP:i:0  xQ:i:0  xC:A:"" xD:A:"" PG:Z:novoalign  AS:i:72 UQ:i:72 NM:i:0  MD:Z:41 PQ:i:190    SM:i:150    AM:i:150    
D43P8DQ1:194:H3W7GADXY:1:2104:5516:41310    147 mm01_24611438_24616266  2360    150 43M2S   =   2276    -133    GTCATCATTGATATATTGTGAGTATATTGGTGAGTAGACCAAGAN   JIGJJJIIJJJJGJGJIJJIJJJJJJIEDJJJJJIHGEGEF?<F~   xl:i:2360   xr:i:2402   xs:i:43 xd:A:r  xm:A:u  xa:A:""     xL:i:2276   xR:i:2402   xS:i:127    xW:i:43 xP:i:0  xQ:i:0  xC:A:"" xD:A:"" PG:Z:novoalign  AS:i:118    UQ:i:118    NM:i:2  MD:Z:22G14G5    PQ:i:190    SM:i:150    AM:i:150    

Задача

Эти строки с разделителями табуляции должны быть прочитаны, и поле xm:A:... должно быть опрошено, чтобы найти значение. Это значение может быть как u, r, так и x. Есть много возможных комбинаций, но нас интересует только несколько. Например:

  • ux - сначала прочитайте u, затем прочитайте x.
  •    
  • rx - сначала прочитанное r, второе - x.
  •    
  • xx - x - и x.

Если строки ux или rx, x всегда будет второй строкой.

После этого мы вводим другой символ в конец последовательности. Например ure или urd Это представляет собой сравнение третьего поля mm01_24611438_24616233 в строках. e означает, что поля должны совпадать, d означает, что они должны быть разными, а a обозначает что-либо.

Для вышеуказанной пары: второе поле совпадает, поэтому оно заканчивается на e. Оба поля xm: A имеют тип u. поэтому правильная комбинация будет uue

В приведенном ниже примере p обозначает, что чтение может быть либо u или r, но не x.

Код

Ниже приведен рабочий фрагмент:

/** Loop through input file and pull out data from the file - types for Paired-end reads
 * @param inputSummary map of MappingTypePE to counts of that type
 * @param inputFile input Xsam file
 * @return Map of String (the mapping type, i.e AAA) to the number of counts for that type
 */
public static LinkedHashMap<String, Integer> mockPopulateWithIncrementingVariablesRestructureDirectStreams(LinkedHashMap<String, Integer> inputSummary, String inputFile) {
    //initialise map


    int aaaCount = 0;
    int paaCount = 0;
    int uueCount = 0;
    int uudCount = 0;
    int rreCount = 0;
    int rrdCount = 0;
    int ureCount = 0;
    int urdCount = 0;
    int uxCount = 0;
    int rxCount = 0;
    int xxCount = 0;


    try {
        BufferedReader fileReader = new BufferedReader(new FileReader(new File(inputFile)));
        String line;
        String line2;

        // /skip past the header
        while((line = fileReader.readLine()) != null){
            if(!line.startsWith("@")){
                if((line2 = fileReader.readLine()) != null){
                    if(percCount == 1000){
                        percCount = 0;
                    }

                    aaaCount++; //always increment anything

                    //get the rnames -> third field
                    String rName1 = line.split("\t")[2];
                    String rName2 = line2.split("\t")[2];

                    //get stats
                    Stream<String> s1 = Stream.of(line.split("\t"));
                    Stream<String> s2 = Stream.of(line2.split("\t"));

                    String mapping1 = s1.filter(d -> d.startsWith("xm"))
                                        .map(res -> res.substring(res.lastIndexOf(':') + 1))
                                        .findFirst()
                                        .get();

                    String mapping2 = s2.filter(d -> d.startsWith("xm"))
                            .map(res -> res.substring(res.lastIndexOf(':') + 1))
                            .findFirst()
                            .get();
                    //paa if first mapping type is not x
                    if(!mapping1.equals("x")){
                        paaCount++;
                    }

                    if(mapping1.equals(mapping2)){ // must be rr or uu
                        //E
                        if(rName1.equals(rName2)){
                            if(mapping1.equals("u")) uueCount++;
                            else rreCount++;
                        }else{
                            //D
                            if(mapping1.equals("u")) uudCount++;
                            else rrdCount++;
                        }
                    }else{ //must be ur or ru
                        if(rName1.equals(rName1)) ureCount++;
                        else urdCount++;
                    }
                    //x cases
                    if(mapping2.equals("x")){
                        switch (mapping1) {
                            case "x":
                                xxCount++;
                                break;
                            case "u":
                                uxCount++;
                                break;
                            default:
                                rxCount++;
                                break;
                        }
                    }

                    percCount++;
                }
            }
        }

        //add the variables to the map
        inputSummary.put("AAA", aaaCount);
        inputSummary.put("PAA", paaCount);
        inputSummary.put("UUE", uueCount);
        inputSummary.put("UUD", uudCount);
        inputSummary.put("RRE", rreCount);
        inputSummary.put("RRD", rrdCount);
        inputSummary.put("URE", ureCount);
        inputSummary.put("URD", urdCount);
        inputSummary.put("UX", uxCount);
        inputSummary.put("RX", rxCount);
        inputSummary.put("XX", xxCount);

    }catch (IOException ioe){
        System.out.println(ioe.getMessage());
    }

    return inputSummary;
}

Бенчмаркинг

Я запустил этот код в 11.8 ГБ файле этих чтений, а общее время выполнения - 112. Я также прочитал один и тот же файл, чтобы узнать, сколько времени потребуется BufferedReader для чтения файла, не делая ничего для строк. Это заняло ~ 28 с. Таким образом, потенциал экономии времени довольно большой.

112s могут показаться недолго, но мы запускаем файлы до 200 ГБ, и этот код должен выполняться до того, как остальная часть программы может работать.

Если у вас есть какие-либо вопросы, пожалуйста, спросите. Извинения за длинный пост!

11 голосов | спросил Sam 8 MaramThu, 08 Mar 2018 11:04:56 +03002018-03-08T11:04:56+03:0011 2018, 11:04:56

5 ответов


11

Чтобы ускорить это, вам нужно избежать как можно большего числа операций создания строк, потому что они дороги. Особенно операция разделения очень дорога. Это не только создает много новых строк, но и делает это в основном без необходимости, потому что вам не нужны все подстроки. Вместо этого вам нужно выполнить поиск низкого уровня, используя только позиции в строке в качестве указателей:

public static void main (String[] args) throws java.lang.Exception
{

    String sample = "D43P8DQ1:194:H3W7GADXY:1:2104:5516:41310\t99\tmm01_24611438_24616266\t2276\t150\t5S41M\t=\t2360\t133\tNNAGGTGAATAGAATTATACCATATCGTAGTCCTTTTTGTACAATA\t~~HHHFHHBGIJDHIFHGGGIIIJJGICGGCBHIIJJHIIEGHCGF\txl:i:2276\txr:i:2316\txs:i:41\txd:A:f\txm:A:u\txa:A:\"\"\txL:i:2276\txR:i:2402\txS:i:127\txW:i:43\txP:i:0\txQ:i:0\txC:A:\"\"\txD:A:\"\"\tPG:Z:novoalign\tAS:i:72\tUQ:i:72\tNM:i:0\tMD:Z:41\tPQ:i:190\tSM:i:150\tAM:i:150";
    char mapping1 = find_xmA(sample);

    System.out.println(mapping1);
}

public static char find_xmA(String sample) {
    int charPos = findPosAfter(sample, "\txm:A:");
    if (charPos == -1) {
        return '\0'; // return NULL character if not found.
    }
    return sample.charAt(charPos);
}

public static int findPosAfter(String haystack, String needle) {
    int hLen = haystack.length();
    int nLen = needle.length();
    int maxSearch = hLen - nLen;

    outer: for (int i = 0; i < maxSearch; i++) {
        for (int j = 0; j < nLen; j++) {
            if (haystack.charAt(i + j) != needle.charAt(j)) {
                continue outer;
            }
        }

        // If it reaches here, match has been found:
        return i + nLen;

    }

    return -1; // Not found
}

Для rName аналогично: найдите индексы второго и третьего символов табуляции в строке и сравните символы между ними один на один, чтобы убедиться, что они равны:

public static void main (String[] args) throws java.lang.Exception
{

    String sample1 = "D43P8DQ1:194:H3W7GADXY:1:2104:5516:41310\t99\tmm01_24611438_24616266\t2276\t150\t5S41M\t=\t2360\t133\tNNAGGTGAATAGAATTATACCATATCGTAGTCCTTTTTGTACAATA\t~~HHHFHHBGIJDHIFHGGGIIIJJGICGGCBHIIJJHIIEGHCGF\txl:i:2276\txr:i:2316\txs:i:41\txd:A:f\txm:A:u\txa:A:\"\"\txL:i:2276\txR:i:2402\txS:i:127\txW:i:43\txP:i:0\txQ:i:0\txC:A:\"\"\txD:A:\"\"\tPG:Z:novoalign\tAS:i:72\tUQ:i:72\tNM:i:0\tMD:Z:41\tPQ:i:190\tSM:i:150\tAM:i:150";
    String sample2 = "D43P8DQ1:194:H3W7GADXY:1:2104:5516:41310\t147\tmm01_24611438_24616266\t2360\t150\t43M2S\t=\t2276\t-133\tGTCATCATTGATATATTGTGAGTATATTGGTGAGTAGACCAAGAN\tJIGJJJIIJJJJGJGJIJJIJJJJJJIEDJJJJJIHGEGEF?<F~\txl:i:2360\txr:i:2402\txs:i:43\txd:A:r\txm:A:u\txa:A:\"\"\txL:i:2276\txR:i:2402\txS:i:127\txW:i:43\txP:i:0\txQ:i:0\txC:A:\"\"\txD:A:\"\"\tPG:Z:novoalign\tAS:i:118\tUQ:i:118\tNM:i:2\tMD:Z:22G14G5\tPQ:i:190\tSM:i:150\tAM:i:150";

    int pos1_1 = findXthChar(sample1, '\t', 2, 0) + 1;
    int pos1_2 = findXthChar(sample1, '\t', 1, pos1_1); // same as just sample1.indexOf('\t', pos1_1)

    int pos2_1 = findXthChar(sample2, '\t', 2, 0) + 1;
    int pos2_2 = findXthChar(sample2, '\t', 1, pos2_1); // same as just sample2.indexOf('\t', pos2_1)

    // Assuming no errors (return value -1) here 

    boolean rNameEqual = areEqualAt(sample1, pos1_1, pos1_2, sample2, pos2_1, pos2_2);

    System.out.println(rNameEqual);
}

private static int findXthChar(String sample, char c, int xth, int fromPos) {
    int pos = sample.indexOf(c, fromPos);
    if (pos == -1) {
        return -1;
    }
    if (xth == 1) {
        return pos;
    }
    return findXthChar(sample, c, xth - 1, pos + 1);
}

private static boolean areEqualAt(String s1, int p11, int p12, String s2, int p21, int p22) {
    int len = p12 - p11;
    if (len != p22 - p21) {
        // Not the same length
        return false;
    }

    for (int i = 0; i < len; i++) {
        if (s1.charAt(p11 + i) != s2.charAt(p21 + i)) {
            return false;
        }
    }

    return true;
}
ответил RoToRa 8 MarpmThu, 08 Mar 2018 12:35:52 +03002018-03-08T12:35:52+03:0012 2018, 12:35:52
10

Если вы работаете с необработанной производительностью, старайтесь не повторять потенциально дорогостоящие операции.

В этом случае вы разделяете строки дважды с тем же параметром, который повторно применяет регулярное выражение под капотом. Вместо

Sring rName1 = line.split("\t")[2];
String rName2 = line2.split("\t")[2];

Stream<String> s1 = Stream.of(line.split("\t"));
Stream<String> s2 = Stream.of(line2.split("\t"));

Разделить один раз и повторно использовать:

String[] splitLine1 = line.split("\t");
String[] splitLine2 = line2.split("\t");

Sring rName1 = splitLine1[2];
String rName2 = splitLine2[2];

Stream<String> s1 = Stream.of(splitLine1);
Stream<String> s2 = Stream.of(splitLine2);

Кроме того, я не вижу большого потенциала для экономии времени. Любопытно видеть измерение после этого изменения ...: -)

ответил mtj 8 MaramThu, 08 Mar 2018 11:37:22 +03002018-03-08T11:37:22+03:0011 2018, 11:37:22
6

Возможная ошибка:

if(mapping1.equals(mapping2)){ // must be rr or uu
...
 if(mapping2.equals("x")){
                    switch (mapping1) {
                        case "x":
                            xxCount++;
                            break;

Основываясь на этом переключателе, включая случай «x», я бы ожидал, что «xx» также является возможной комбинацией, что означает, что ваш комментарий ранее неверен. Вы увеличиваете и код xxCount и rreCount /rrdCount? Это намеренно?


Я не уверен, что у вас действительно есть такой потенциал, чтобы ускорить процесс, как вы думаете. На каждой итерации этого цикла while вы фактически просматриваете строки, разделенные вкладками обеих строк, которые звучат как «работа» для компьютера. Я чувствую, что ты не собираешься приближаться к 28-м рабочему времени.

Единственное «очевидное», что я мог найти, это то, что вы разделили каждую строку дважды:

                String rName1 = line.split("\t")[2];
                ...
                Stream<String> s1 = Stream.of(line.split("\t"));

Это может помочь, если вы сохраните результат line.split("\t") в переменную и используйте его для обоих этих операторов.

Если вы используете профилировщик, чтобы узнать, где ваш код занимает больше всего времени, это может помочь, если вы поместите эти строки в отдельный метод:

            Stream<String> s1 = Stream.of(line.split("\t"));
            String mapping1 = s1.filter(d -> d.startsWith("xm"))
                                .map(res -> res.substring(res.lastIndexOf(':') + 1))
                                .findFirst()
                                .get();

Вы также можете использовать этот метод для mapping1 и mapping2 при передаче в списке строк из split.


Еще одна меньшая оптимизация заключалась бы в том, чтобы перебирать строки в простом цикле for вместо использования потока. Поток создает дополнительные накладные расходы.

public static String parseMapping(String[] line){
    for(String word : line){
        if (word.startsWith("xm")) {
            return word.substring(word.lastIndexOf(':') + 1);
        }
    }
    return null; // handle wrong file? can't happen?
}

Хотя я понятия не имею, сколько это выиграет.

ответил Imus 8 MaramThu, 08 Mar 2018 11:42:33 +03002018-03-08T11:42:33+03:0011 2018, 11:42:33
1

На самом деле это не связано с вашим кодом, но я думал, что упомянул об этом, поскольку несколько комментаторов уже обсуждали его.

Если ваш файл занимает 28 секунд для чтения для 11,8 ГБ файла, это примерно 431 МБ в секунду. Это примерно скорость SSD SSD, возможно, немного ниже, поэтому я предполагаю, что вы используете это.

Если это вообще возможно, я бы рекомендовал прочитать ваш файл с SSD PCIe M.2. Есть несколько продавцов, которые продают SSD этого сорта, которые в соответствии с UserBenchmark получают скорость около 2200 МБ в секунду для последовательного чтения. Это примерно в 5 раз превышает текущую скорость чтения. Теоретически, ваш файл размером 11,8 ГБ займет всего около 5,6 секунд для чтения. После того, как вы собираетесь использовать 200-гигабайтные файлы, то до примерно 475 секунд займет примерно 93 секунды, сохранив около 6 минут на файл. PCIe SSD с объемом памяти 500 ГБ относительно дешевы, около 200 долларов США. тот же SSD с объемом памяти 256 ГБ стоит около 100 долларов США.

ответил Nzall 9 MarpmFri, 09 Mar 2018 12:11:35 +03002018-03-09T12:11:35+03:0012 2018, 12:11:35
0

Чтобы выяснить, сколько времени занимает ввод-вывод, удалите все между циклами while((line = fileReader.readLine()) != null){. Значит, вы знаете, где вы можете сосредоточиться в первую очередь. Или там, где это не имеет никакого смысла.

Я читаю 829 МБ примерно за две секунды на моей машине (используя свой способ загрузить файл без магии в цикле while). Если я использую BufferedInputStream-> FileInputStream, это составляет половину секунды, используя буфер 8 * 1024. Ну, это работает до тех пор, пока вы не работаете с фантастическими кодировками.

Вот интересная статья, один парень проверил множество способов загрузки данных из потока: https://stackoverflow.com/questions/309424/read-convert-an- InputStream к струне .

Вот интересная статья о размерах буфера: https://stackoverflow.com/questions/236861 /как-делать-вы-определения-The-идеально-буфера размера, когда-используя-FileInputStream

Если ничего не помогает: если вы отвечаете за часть, которая записывает файл: записывайте разные файлы на разных жестких дисках и создавайте поток для каждого файла и объединяйте результаты вместе.

С уважением, slowy

ответил slowy 14 MarpmWed, 14 Mar 2018 23:12:30 +03002018-03-14T23:12:30+03:0011 2018, 23:12:30

Похожие вопросы

Популярные теги

security × 330linux × 316macos × 2827 × 268performance × 244command-line × 241sql-server × 235joomla-3.x × 222java × 189c++ × 186windows × 180cisco × 168bash × 158c# × 142gmail × 139arduino-uno × 139javascript × 134ssh × 133seo × 132mysql × 132