# Counting relevant entries in a large bioinformatics file

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:

• 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
*/
//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 {
String line;
String line2;

if(!line.startsWith("@")){
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!

• Please do not update the code in your question to incorporate feedback from answers, doing so goes against the Question + Answer style of Code Review. This is not a forum where you should keep the most updated version in your question. Please see what you may and may not do after receiving answers. – Zeta Mar 8 '18 at 9:18
• Speed optimisation should always start by an analysis of what is taking time. Run your program through a cpu sampling tool (like jvisualvm, included in the oracle jvm). It will tell you what is taking time, so you can focus on optimizing that and nothing else (as speed optimization tends to decrease readability). – Thierry Mar 8 '18 at 9:49
• @Zeta Is that a standard template? Because it doesn't say anything about invalidating the posted (or work-in-progress) answers, which is the actual reason we don't edit code in questions. – Fund Monica's Lawsuit Mar 8 '18 at 23:44
• @NicHartley Somewhat. There are other variants, see codereview.meta.stackexchange.com/a/4954/21002. – Zeta Mar 9 '18 at 6:17
• If performance is key, consider next to code optimalization to use a faster hard disk (SSD) to read the file from. – Thomas Mar 9 '18 at 8:00

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
{

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

}

}


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
{

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

• Really impressive. I've implemented this change and run the same file. The execution time is 43.5s. It's definitely more code to maintain, but the time saving is massive. Thanks! – Sam Mar 8 '18 at 9:56
• Here's a perfect example of trading readability for performance. Seeing how your main concern is performance it's most likely worth it using this code. To reduce your fear to maintain it later on it's probably a good idea to add extra comments to explain why and what this does. And maybe rename variables that better fit how you would write it yourself. A +1 for this answer from me in any case. – Imus Mar 8 '18 at 13:20
• As an additional note, for further performance you could go multi-process fairly easily. Read in the file, and then allocate equal sized chunks of lines between several different processes. Each can then search it's own section. Depending on the number of cores your computer has this could be anything up to a 8x speed improvement (not including the file reading stage, as presumably this part is IO limited not search limited). – Drgabble Mar 8 '18 at 16:02
• Yes, reading the file is probably the slowest part now. Try running the programm without any processing, just the reading of the file. That's basically the fastest you can get. If you can try it on a computer with an SSD. – RoToRa Mar 8 '18 at 18:37

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... :-)

• Thanks for the suggestion, that's actually reduced the time down to 70 seconds. I had no idea splitting a line was so intensive! – Sam Mar 8 '18 at 8:46
• So 28s to read the file line by line. 40s gained from splitting only once, meaning 40s still need for that 1 split. You're probably not going to get better than that 70 seconds. – Imus Mar 8 '18 at 8:50
• Here's an easy optimization: Do Pattern regex = Pattern.compile("\t"); outside of the loop and then use regex.split(line) instead of line.split("\t"). This way, you compile the regex only once per file instead of twice per line. Going even further, you could convert it into a class variable, then the regex would be compiled only once for the lifetime of your application (generally). – xehpuk Mar 9 '18 at 8:25

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.

• Good catch on the bug. That would indeed lead to rrd being incremented when the lines are 'xx'. mtj also picked up on the second split, and it seems to have been a considerable decrease, down to 70s. I've implemented the parseMapping function, and it's got the time down to 63s, another 7 seconds off. Thanks – Sam Mar 8 '18 at 8:53
• Oh wow, I didn't expect it to help that much still. Great. – Imus Mar 8 '18 at 9:29
• As for the bug, there might also be a similar one in that else part. I suggest you create your own really small input file that has each of the cases and check if the total counts are correct. It's really easy to overlook a case in the code as you've seen. – Imus Mar 8 '18 at 9:33

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.

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