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


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

Файл

Тип файла находится в Доме творчества, называется файл 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 или Х. Есть много возможных комбинаций, но нас интересуют только некоторые из них. Например:

  • УБ - первых читать U, второе чтение х.
  • ГХ - сначала прочитайте R, второе чтение х.
  • ХХ - х - и X.

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

После этого мы введем еще один символ конца последовательности. Например ure или urd это сравнение третье поле mm01_24611438_24616233 в струнах. e обозначает поля должны совпадать, d означает, они должны быть разница и 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С. Так что потенциал для экономии времени достаточно большие.

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

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



Комментарии
4 ответа

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

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;
}

11
ответ дан 8 марта 2018 в 09:03 Источник Поделиться

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

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

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);

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

10
ответ дан 8 марта 2018 в 08:03 Источник Поделиться

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

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

На основе этого переключателя в том случае "х" я бы ожидал, что "ХХ" является также возможной комбинации, а значит ваш комментарий раньше на это не так. Вы приращение и 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?
}

Хотя я понятия не имею, сколько это будет получить.

6
ответ дан 8 марта 2018 в 08:03 Источник Поделиться

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

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

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

1
ответ дан 9 марта 2018 в 09:03 Источник Поделиться