scFactory.captureUnmappedScores = true;
scFactory.losslessQS = true;
scFactory.preserveReadNames = true;
SequenceBaseProvider provider = new ByteArraySequenceBaseProvider(refBytes);
SAMFileReader reader = new SAMFileReader(bamFile);
SAMRecordIterator iterator = reader.query(seqName, 0, 0, false);
List<CramRecord> records = new ArrayList<CramRecord>();
PairedTemplateAssembler a = new PairedTemplateAssembler();
int counter = 0;
int bases = 0;
long size = 0;
Definition definition = new Definition();
definition.contentId = "input.bam".getBytes();
definition.formatMajor = 1;
definition.formatMinor = 0;
definition.magick = "CRAM".getBytes();
Format format = new FormatFactory().createFormat(definition);
SAMRecord samRecord = null;
CramHeader cramHeader = createCramHeader(reader.getFileHeader());
ByteArrayOutputStream hBaos= new ByteArrayOutputStream() ;
Container hC = new Container() ;
Block hBlock = new Block() ;
hBlock.contentId=0;
hBlock.contentType=0;
hBlock.method = CompressionMethod.GZIP.byteValue() ;
hC.blocks = new Block[]{hBlock} ;
hC.containers= new Container[0] ;
CramHeaderIO.write(cramHeader, hBaos) ;
hBlock.data = hBaos.toByteArray() ;
format.writeContainer(hC, hBaos) ;
System.out.println("Header size: " + hBaos.size());
for (int b = 0; b < maxRecords; b+=recordsPerBlock) {
bases = 0 ;
size = 0 ;
for (int i = 0; i < recordsPerBlock; i++) {
samRecord = iterator.next();
a.addSAMRecord(samRecord);
bases += samRecord.getReadLength();
}
for (int i = 0; i < recordsPerBlock; i++) {
while ((samRecord = a.nextSAMRecord()) != null) {
CramRecord cramRecord = scFactory.createCramRecord(samRecord);
records.add(cramRecord);
setPairingInformation(cramRecord, samRecord, a.distanceToNextFragment());
}
while ((samRecord = a.fetchNextSAMRecord()) != null) {
CramRecord cramRecord = scFactory.createCramRecord(samRecord);
records.add(cramRecord);
setPairingInformation(cramRecord, samRecord, a.distanceToNextFragment());
}
}
Container container = new Coder(recordsPerSlice).writeRecords(cramHeader, records, provider);
ByteArrayOutputStream baos = new ByteArrayOutputStream();
long time1 = System.nanoTime();
format.writeContainer(container, baos);
long time2 = System.nanoTime();
System.out.printf("Container written: bytes=%d in %.2f ms.\n", baos.size(), (time2 - time1) / 1000000f);
size += baos.size();
records.clear();
System.out.printf("Total records: %d, bases=%d, bytes=%d, b/b=%.2f\n", counter, bases, size, 8f * size
/ bases);
time1 = System.nanoTime();
Container c = format.readContainer(new ByteArrayInputStream(baos.toByteArray()));
time2 = System.nanoTime();
System.out.printf("Container read in %.2f ms.\n", (time2 - time1) / 1000000f);
time1 = System.nanoTime();
List<CramRecord> readRecords = new Coder(recordsPerSlice).readRecords(c, cramHeader, provider);
time2 = System.nanoTime();
System.out.printf("Cram records read in in %.2f ms.\n", (time2 - time1) / 1000000f);
System.out.println(readRecords.size());
RestoreBases rb = new RestoreBases(provider, "20");
time1 = System.nanoTime();
for (CramRecord r : readRecords)
rb.restoreReadBases(r);
time2 = System.nanoTime();
System.out.printf("Bases restored in in %.2f ms.\n", (time2 - time1) / 1000000f);
RestoreQualityScores rq = new RestoreQualityScores();
time1 = System.nanoTime();
for (CramRecord r : readRecords)
rq.restoreQualityScores(r);
time2 = System.nanoTime();
System.out.printf("Scores restored in in %.2f ms.\n", (time2 - time1) / 1000000f);
time1 = System.nanoTime();
restoreNamesAndSpots("name", 1, readRecords);
time2 = System.nanoTime();
System.out.printf("Pairing restored in in %.2f ms.\n", (time2 - time1) / 1000000f);
// for (int i = 0; i < 10; i++)
// System.out.println(readRecords.get(i));
// for (int i = readRecords.size() - 10; i < readRecords.size();
// i++)
// System.out.println(readRecords.get(i));
time1 = System.nanoTime();
SAMFileHeader header = Utils.cramHeader2SamHeader(cramHeader);
List<SAMRecord> samRecords = convert(readRecords, header, null, seqName);
for (SAMRecord r : samRecords)
Utils.calculateMdAndNmTags(r, refBytes, true, true);
time2 = System.nanoTime();
System.out.printf("SAMRecords restored in in %.2f ms.\n", (time2 - time1) / 1000000f);
// for (int i = 0; i < 10; i++)
// System.out.print(samRecords.get(i).getSAMString());
// for (int i = samRecords.size() - 10; i < samRecords.size(); i++)
// System.out.print(samRecords.get(i).getSAMString());
}
reader.close();
}