1717import htsjdk .variant .vcf .VCFHeader ;
1818import htsjdk .variant .vcf .VCFHeaderLineType ;
1919import htsjdk .variant .vcf .VCFInfoHeaderLine ;
20+ import org .apache .commons .lang3 .ArrayUtils ;
2021import org .apache .commons .lang3 .StringUtils ;
2122import org .apache .commons .lang3 .math .NumberUtils ;
23+ import org .apache .commons .math3 .stat .descriptive .rank .Median ;
2224import org .apache .log4j .Logger ;
2325import org .biojava3 .core .sequence .DNASequence ;
2426import org .biojava3 .core .sequence .compound .AmbiguityDNACompoundSet ;
@@ -238,7 +240,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
238240 File ampliconDepths = new File (outputDir , "ampliconDepths.txt" );
239241 try (CSVWriter writer = new CSVWriter (IOUtil .openFileForBufferedUtf8Writing (ampliconDepths ), '\t' , CSVWriter .NO_QUOTE_CHARACTER ))
240242 {
241- writer .writeNext (new String []{"AmpliconName" , "Length" , "TotalDepth" , "AvgDepth" });
243+ writer .writeNext (new String []{"AmpliconName" , "Length" , "TotalDepth" , "AvgDepth" , "MedianDepth" });
242244 File primerFile = getPipelineCtx ().getSequenceSupport ().getCachedData (ampliconBedId );
243245 try (CSVReader reader = new CSVReader (Readers .getReader (primerFile ), '\t' ))
244246 {
@@ -256,6 +258,7 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
256258 String name = line [3 ];
257259
258260 int totalDepth = 0 ;
261+ List <Double > depths = new ArrayList <>();
259262 for (int i = start ;i <=end ;i ++)
260263 {
261264 String key = contig + ":" + i ;
@@ -266,11 +269,14 @@ public Output performAnalysisPerSampleRemote(Readset rs, File inputBam, Referenc
266269 }
267270
268271 totalDepth += gatkDepth .get (key );
272+ depths .add (gatkDepth .get (key ).doubleValue ());
269273 }
270274
271275 int length = (end - start + 1 );
272276 Double avg = (double )totalDepth / length ;
273- writer .writeNext (new String []{name , String .valueOf (length ), String .valueOf (totalDepth ), decimal .format (avg )});
277+ Median median = new Median ();
278+ double medianValue = median .evaluate (ArrayUtils .toPrimitive (depths .toArray (new Double [0 ])));
279+ writer .writeNext (new String []{name , String .valueOf (length ), String .valueOf (totalDepth ), decimal .format (avg ), decimal .format (medianValue )});
274280 }
275281 }
276282 }
0 commit comments