Counting relevant entries in a large bioinformatics file

Question

I have a working example of a piece of code which opens up a file, gathers information about the contents, and outputs a map which contains the information.

The file

The file type is an in-house creation, called an Xsam file. For those interested, it's based on the sam file, which is used commonly in bioinformatics. Each files starts with a header section, of which each line starts with "@" and can be safely ignored by this -> there are usually no more than 1000 lines in the header. The rest of the file is made up of read-pairs. Each read takes up a single line, and the lines are always in pairs. An example of a read pair:

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    

The challenge

These tab-delimited lines must be read, and the xm:A:... field should be interrogated to find the value. This value can either be u, r or x. There are many possible combinations, but we're interested in just a few. For example:

• ux - first read u, second read x.
• rx - first read r, second read x.
• xx - x - and x.

If the lines are ux or rx, the x will always be the second line.

After that, we introduce another character to the end of the sequence. For example ure or urd This represents a comparison of the third field mm01_24611438_24616233 in the strings. e denotes the fields must match, d denotes they must be difference and a denotes anything.

For the above pair: the second field matches, so it ends in e. Both xm:A fields are of u type. so the correct combination would be uue

In the example below, a p denotes the read can either be a u or an r but not an x.

The code

Below is a working snippet:

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

Benchmarking

I ran this code on an 11.8GB file of these reads, and the overall time of execution was ~112s. I have also read through the same file to see how long a BufferedReader would take to read the file without doing anything to the lines. This took ~28s. So the potential for time saving is quite large.

The 112s may not seem like a long time, but we run files up to 200GB, and this code must execute before the rest of the program can run.

If you have any questions, please ask. Apologies for the long post!

Answer 1

In order to speed this up, you'll need to avoid as many string creation operations as possible, because they are expensive. Especially the split operation is expensive. Not only does this create many new strings, it does this mostly unnecessarily, because you don't need all the substrings. Instead you need to do some low level searching using only positions in the string as pointers:

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
}

For rName it similarly: Look for the indexes of the second and third tab characters in the string, and compare the characters between them one for one, to see if they are equal:

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

Answer 2

If you are in for raw performance, try to avoid repeating potentially cost-intensive operations.

In this case, you split the lines twice with the same parameter, which repeatedly applies a regular expression under the hood. Instead of

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

Split once and reuse:

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

Apart from that I don't see much potential for time-saving. Curious to see the measurement after that change... :-)

Answer 3

Possible bug:

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

Based on that switch including the case "x" I would expect that "xx" is also a possible combination, which means your comment earlier on is wrong. Do you increment both the xxCount and rreCount/rrdCount ? Is this intentional?

---

I'm not convinced you really have that much potential to speed things up as you think. On each iteration of that while loop you're actually looping through the tab-separated strings of both lines which sounds like a lot of "work" for the computer to do. My gut feeling says you're not going to get close to the 28s runtime.

The only "obvious" thing I could find was that you're splitting each line twice:

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

It might help if you store the result of line.split("\t") into a variable and use it for both those statements.

If you're using a profiler to see where your code is taking most time it might help if you put these lines into a separate method:

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

You can also use that method for both mapping1 and mapping2 when passing in the list of Strings from the split.

---

Another smaller optimisation would be to to loop over the strings in a simple for loop instead of using a stream. The stream causes extra overhead.

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

Although I have no idea how much this would gain.

Answer 4

Not really related to your code, but thought I'd mention it since several commenters discussed it already.

If your file takes 28 seconds to read for an 11.8 GB file, that's roughly 431 MB per second. This is roughly the speed of a SATA SSD, maybe slightly lower, so I assume you're using that.

If at all possible, I would recommend reading your file from a PCIe M.2 SSD. There are multiple vendors who sell SSDs of this variety that according to UserBenchmark get speeds around 2200 MB per second sequential read. That is roughly 5 times your current read speed. Theoretically, your 11.8 GB file will then only take around 5.6 seconds to read. Once you are going to use 200 GB files, what would before take roughly 475 seconds will take roughly 93 seconds, saving you about 6 minutes per file. PCIe SSDs with a storage capacity of 500 GB are relatively cheap, around 200 USD. the same SSD with a storage capacity of 256 GB costs around 100 USD.

Answer 5

To figure out, how long the i/o takes, remove everything between within the while((line = fileReader.readLine()) != null){ loop. So you know where you can focus first. Or where it doesn't make any sense.

I'm reading 829MB's in around two seconds on my machine(Using your way to load the file, without the magic within the while loop). If I use BufferedInputStream->FileInputStream, it's half a second, using a 8*1024 buffer. Well, that works as long as you don't work with fancy charsets.

Here's an interesting article, one guy did a test a lot of ways to load data from a stream: https://stackoverflow.com/questions/309424/read-convert-an-inputstream-to-a-string.

Here's an interesting article about buffer sizes: https://stackoverflow.com/questions/236861/how-do-you-determine-the-ideal-buffer-size-when-using-fileinputstream

If nothing helps: If you're in charge of the part which writes the file: Write different files on different harddisks and have a thread for every file and merge the results together.

Regards, slowy