1111import htsjdk .samtools .util .CloseableIterator ;
1212
1313public class BAMUtilities {
14- public static double calculateStandardizationRatio (File BAM ) throws IOException {
14+
15+ public static double calculateStandardizationRatio (File BAM , int read ) throws IOException {
1516 SamReader inputSam = SamReaderFactory .makeDefault ().open (BAM );
1617 AbstractBAMFileIndex bai = (AbstractBAMFileIndex ) inputSam .indexing ().getIndex ();
17- double counter = 0 ;
1818 double totalAligned = 0 ;
1919 double totalGenome = 0 ;
2020
@@ -23,20 +23,34 @@ public static double calculateStandardizationRatio(File BAM) throws IOException
2323 totalAligned += inputSam .indexing ().getIndex ().getMetaData (x ).getAlignedRecordCount ();
2424 totalGenome += seq .getSequenceLength ();
2525 }
26+
27+ double READ1 = 0 ;
28+ double READ2 = 0 ;
29+ double MID = 0 ;
2630 CloseableIterator <SAMRecord > iter = inputSam .iterator ();
2731 while (iter .hasNext ()) {
2832 SAMRecord sr = iter .next ();
29- if (sr .getReadPairedFlag ()) {
30- if (sr .getSecondOfPairFlag ()) { counter ++; } //count read 2 to remove from aligned reads
33+ if (!sr .getReadUnmappedFlag ()) { //Test for mappability
34+ if (sr .getReadPairedFlag ()) { //Test for paired-end status
35+ if (sr .getSecondOfPairFlag ()) { READ2 ++; } //count read 2
36+ else if (sr .getFirstOfPairFlag ()) { READ1 ++; } // count read 1
37+ if (sr .getProperPairFlag () && sr .getFirstOfPairFlag ()) { MID ++; } //count properly paired reads
38+ } else { //If the read is mapped but not paired-end, default to read 1
39+ READ1 ++;
40+ }
3141 }
3242 }
43+ iter .close ();
3344 inputSam .close ();
3445 bai .close ();
35- iter .close ();
3646
37- totalAligned -= counter ;
47+ //System.out.println("Genome Size: " + totalGenome + "\nTotal tags: " + totalAligned + "\nDetected Read 1: " + READ1 + "\nDetected Read 2: " + READ2 + "\nDetected Midpoints: " + MID);
48+ if (read == 0 ) { totalAligned = READ1 ; }
49+ else if (read == 1 ) { totalAligned = READ2 ; }
50+ else if (read == 2 ) { totalAligned = READ1 + READ2 ; }
51+ else if (read == 3 ) { totalAligned = MID ; }
52+
3853 if (totalAligned > 0 ) { return (totalAligned / totalGenome ); }
3954 else { return 1 ; }
4055 }
41-
4256}
0 commit comments