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!