extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/vcf/VcfData.java

Code
Comments
Other
Rev Date Author Line
4620 16 Nov 17 nicklas 1 package net.sf.basedb.reggie.vcf;
4620 16 Nov 17 nicklas 2
4646 14 Dec 17 nicklas 3 import java.io.Serializable;
4620 16 Nov 17 nicklas 4 import java.util.ArrayList;
4620 16 Nov 17 nicklas 5 import java.util.LinkedHashMap;
4620 16 Nov 17 nicklas 6 import java.util.List;
4620 16 Nov 17 nicklas 7 import java.util.Map;
4620 16 Nov 17 nicklas 8
4620 16 Nov 17 nicklas 9 import org.json.simple.JSONObject;
4620 16 Nov 17 nicklas 10
4620 16 Nov 17 nicklas 11 /**
4620 16 Nov 17 nicklas 12   Summary of data found in a VCF.
4620 16 Nov 17 nicklas 13   @author nicklas
4620 16 Nov 17 nicklas 14   @since 4.13
4620 16 Nov 17 nicklas 15 */
4620 16 Nov 17 nicklas 16 public class VcfData 
4646 14 Dec 17 nicklas 17   implements Serializable
4620 16 Nov 17 nicklas 18 {
4620 16 Nov 17 nicklas 19
4646 14 Dec 17 nicklas 20   private static final long serialVersionUID = -3339783464081195474L;
7422 14 Nov 23 nicklas 21
7422 14 Nov 23 nicklas 22   private final boolean paired;
4646 14 Dec 17 nicklas 23   
6444 19 Oct 21 nicklas 24   private int totalCount;
4620 16 Nov 17 nicklas 25   private int homRefCount;
4620 16 Nov 17 nicklas 26   private int homAltCount;
4620 16 Nov 17 nicklas 27   private int hetCount;
4620 16 Nov 17 nicklas 28   
4620 16 Nov 17 nicklas 29   private int badCount;
7414 10 Nov 23 nicklas 30   private int phasedCount;
4620 16 Nov 17 nicklas 31   
6513 07 Dec 21 nicklas 32   private float sumHomGQ;
6513 07 Dec 21 nicklas 33   private float sumHetGQ;
4620 16 Nov 17 nicklas 34   
6513 07 Dec 21 nicklas 35   private int highGqCount;
6513 07 Dec 21 nicklas 36   private int lowGqCount;
4648 15 Dec 17 nicklas 37   
7422 14 Nov 23 nicklas 38   private Map<String, GtData> tumorGt;
7422 14 Nov 23 nicklas 39   private Map<String, GtData> normalGt;
6450 21 Oct 21 nicklas 40   private JSONObject infoHeaders;
6450 21 Oct 21 nicklas 41   private JSONObject formatHeaders;
6513 07 Dec 21 nicklas 42   private GenotypeQualityModel qualityModel;
4620 16 Nov 17 nicklas 43   
4712 22 Mar 18 nicklas 44   private int fileId;
4712 22 Mar 18 nicklas 45     
7422 14 Nov 23 nicklas 46   public VcfData(boolean paired)
4620 16 Nov 17 nicklas 47   {
7422 14 Nov 23 nicklas 48     this.paired = paired;
6513 07 Dec 21 nicklas 49     this.qualityModel = null;
7422 14 Nov 23 nicklas 50     this.tumorGt = new LinkedHashMap<>();
7422 14 Nov 23 nicklas 51     this.normalGt = new LinkedHashMap<>();
6450 21 Oct 21 nicklas 52     this.infoHeaders = new JSONObject();
6450 21 Oct 21 nicklas 53     this.formatHeaders = new JSONObject();
4620 16 Nov 17 nicklas 54   }
4620 16 Nov 17 nicklas 55   
7422 14 Nov 23 nicklas 56   /**
7422 14 Nov 23 nicklas 57     Is this VCF paired normal/tumor data or not.
7422 14 Nov 23 nicklas 58     @since 4.50
7422 14 Nov 23 nicklas 59   */
7422 14 Nov 23 nicklas 60   public boolean isPaired()
7422 14 Nov 23 nicklas 61   {
7422 14 Nov 23 nicklas 62     return paired;
7422 14 Nov 23 nicklas 63   }
7422 14 Nov 23 nicklas 64   
6444 19 Oct 21 nicklas 65   public void setTotalCount(int totalCount)
6444 19 Oct 21 nicklas 66   {
6444 19 Oct 21 nicklas 67     this.totalCount = totalCount;
6444 19 Oct 21 nicklas 68   }
6444 19 Oct 21 nicklas 69   
6444 19 Oct 21 nicklas 70   public int getTotalCount()
6444 19 Oct 21 nicklas 71   {
6444 19 Oct 21 nicklas 72     return totalCount;
6444 19 Oct 21 nicklas 73   }
6444 19 Oct 21 nicklas 74   
6450 21 Oct 21 nicklas 75   /**
6450 21 Oct 21 nicklas 76     Adds information about an ##INFO header.
6450 21 Oct 21 nicklas 77     @since 4.34
6450 21 Oct 21 nicklas 78   */
6450 21 Oct 21 nicklas 79   public void addInfoHeader(String id, String data)
6450 21 Oct 21 nicklas 80   {
6450 21 Oct 21 nicklas 81     infoHeaders.put(id, data);
6450 21 Oct 21 nicklas 82   }
6444 19 Oct 21 nicklas 83   
4620 16 Nov 17 nicklas 84   /**
6450 21 Oct 21 nicklas 85     Adds information about a ##FORMAT header.
6450 21 Oct 21 nicklas 86     @since 4.34
6450 21 Oct 21 nicklas 87   */
6450 21 Oct 21 nicklas 88   public void addFormatHeader(String id, String data)
6450 21 Oct 21 nicklas 89   {
6450 21 Oct 21 nicklas 90     formatHeaders.put(id, data);
6513 07 Dec 21 nicklas 91     if ("GQ".equals(id))
6513 07 Dec 21 nicklas 92     {
6513 07 Dec 21 nicklas 93       if (data.contains("GenCall Score"))
6513 07 Dec 21 nicklas 94       {
6513 07 Dec 21 nicklas 95         qualityModel = GenotypeQualityModel.GENCALL;
6513 07 Dec 21 nicklas 96       }
6513 07 Dec 21 nicklas 97       else if (data.contains("Genotype Quality"))
6513 07 Dec 21 nicklas 98       {
6513 07 Dec 21 nicklas 99         qualityModel = GenotypeQualityModel.HAPLOTYPECALLER;
6513 07 Dec 21 nicklas 100       }
6513 07 Dec 21 nicklas 101     }
6450 21 Oct 21 nicklas 102   }
6450 21 Oct 21 nicklas 103   
6450 21 Oct 21 nicklas 104   /**
6513 07 Dec 21 nicklas 105     Set the quality score model to use.
6513 07 Dec 21 nicklas 106     @since 4.34
6513 07 Dec 21 nicklas 107   */
6513 07 Dec 21 nicklas 108   public void setQualityModel(GenotypeQualityModel qualityModel)
6513 07 Dec 21 nicklas 109   {
6513 07 Dec 21 nicklas 110     this.qualityModel = qualityModel;
6513 07 Dec 21 nicklas 111   }
6513 07 Dec 21 nicklas 112   
6513 07 Dec 21 nicklas 113   /**
6513 07 Dec 21 nicklas 114     Get the quality score model in use.
6513 07 Dec 21 nicklas 115     @since 4.34
6513 07 Dec 21 nicklas 116   */
6513 07 Dec 21 nicklas 117   public GenotypeQualityModel getQualityModel()
6513 07 Dec 21 nicklas 118   {
6513 07 Dec 21 nicklas 119     return qualityModel;
6513 07 Dec 21 nicklas 120   }
6513 07 Dec 21 nicklas 121
6513 07 Dec 21 nicklas 122   /**
6450 21 Oct 21 nicklas 123     Get information about an ##INFO header.
6450 21 Oct 21 nicklas 124     @since 4.34
6450 21 Oct 21 nicklas 125   */
6450 21 Oct 21 nicklas 126   public String getInfoHeader(String id)
6450 21 Oct 21 nicklas 127   {
6450 21 Oct 21 nicklas 128     return (String)infoHeaders.get(id);
6450 21 Oct 21 nicklas 129   }
6450 21 Oct 21 nicklas 130   
6450 21 Oct 21 nicklas 131   /**
6450 21 Oct 21 nicklas 132     Get information about a ##FORMAT header.
6450 21 Oct 21 nicklas 133     @since 4.34
6450 21 Oct 21 nicklas 134   */
6450 21 Oct 21 nicklas 135   public String getFormatHeader(String id)
6450 21 Oct 21 nicklas 136   {
6450 21 Oct 21 nicklas 137     return (String)formatHeaders.get(id);
6450 21 Oct 21 nicklas 138   }
6450 21 Oct 21 nicklas 139   
6450 21 Oct 21 nicklas 140   /**
4620 16 Nov 17 nicklas 141     Adds genotype information for the given SNP. If
4620 16 Nov 17 nicklas 142     the data has no Genotype, it is ignored and only
4620 16 Nov 17 nicklas 143     the 'bad' counter is updated.
4620 16 Nov 17 nicklas 144     
4620 16 Nov 17 nicklas 145     @param id rsId for the SNP
4620 16 Nov 17 nicklas 146     @param gt Genotype information
4620 16 Nov 17 nicklas 147   */
4620 16 Nov 17 nicklas 148   public void addGt(GtData gt)
4620 16 Nov 17 nicklas 149   {
7422 14 Nov 23 nicklas 150     addGt(null, gt);
7422 14 Nov 23 nicklas 151   }
7422 14 Nov 23 nicklas 152   
7422 14 Nov 23 nicklas 153   /**
7422 14 Nov 23 nicklas 154     Adds genotype information for the given SNP. Statistics if
7422 14 Nov 23 nicklas 155     collected for the tumor genotype information. If
7422 14 Nov 23 nicklas 156     the data has no Genotype, it is ignored and only
7422 14 Nov 23 nicklas 157     the 'bad' counter is updated.
4648 15 Dec 17 nicklas 158     
7422 14 Nov 23 nicklas 159     @param normal Genotype information for the normal sample (or null if no normal)
7422 14 Nov 23 nicklas 160     @param tumor Genotype information for the tumor
7422 14 Nov 23 nicklas 161   */
7422 14 Nov 23 nicklas 162   public void addGt(GtData normal, GtData tumor)
7422 14 Nov 23 nicklas 163   {
7422 14 Nov 23 nicklas 164     GenoType genoType = tumor.getGenoType();
7422 14 Nov 23 nicklas 165     float gq = tumor.getGQ();
7422 14 Nov 23 nicklas 166     
4620 16 Nov 17 nicklas 167     if (genoType == GenoType.HOM_REF)
4620 16 Nov 17 nicklas 168     {
4620 16 Nov 17 nicklas 169       homRefCount++;
4648 15 Dec 17 nicklas 170       sumHomGQ += gq;
4620 16 Nov 17 nicklas 171     }
4620 16 Nov 17 nicklas 172     else if (genoType == GenoType.HOM_ALT)
4620 16 Nov 17 nicklas 173     {
4620 16 Nov 17 nicklas 174       homAltCount++;
4648 15 Dec 17 nicklas 175       sumHomGQ += gq;
4620 16 Nov 17 nicklas 176     }
4620 16 Nov 17 nicklas 177     else if (genoType == GenoType.HET)
4620 16 Nov 17 nicklas 178     {
4620 16 Nov 17 nicklas 179       hetCount++;
4648 15 Dec 17 nicklas 180       sumHetGQ += gq;
4620 16 Nov 17 nicklas 181     }
4620 16 Nov 17 nicklas 182     else
4620 16 Nov 17 nicklas 183     {
4620 16 Nov 17 nicklas 184       badCount++;
4620 16 Nov 17 nicklas 185       return;
4620 16 Nov 17 nicklas 186     }
7422 14 Nov 23 nicklas 187     if (tumor.isPhased()) phasedCount++;
4620 16 Nov 17 nicklas 188     
6513 07 Dec 21 nicklas 189     if (qualityModel != null)
4648 15 Dec 17 nicklas 190     {
6513 07 Dec 21 nicklas 191       if (qualityModel.isHigh(gq)) 
6513 07 Dec 21 nicklas 192       {
6513 07 Dec 21 nicklas 193         highGqCount++;
6513 07 Dec 21 nicklas 194       }
6513 07 Dec 21 nicklas 195       else if (qualityModel.isLow(gq))
6513 07 Dec 21 nicklas 196       {
6513 07 Dec 21 nicklas 197         lowGqCount++;
6513 07 Dec 21 nicklas 198       }
4648 15 Dec 17 nicklas 199     }
4648 15 Dec 17 nicklas 200     
7422 14 Nov 23 nicklas 201     this.tumorGt.put(tumor.getId(), tumor);
7422 14 Nov 23 nicklas 202     if (normal != null) normalGt.put(normal.getId(), normal);
4620 16 Nov 17 nicklas 203   }
4620 16 Nov 17 nicklas 204   
4620 16 Nov 17 nicklas 205   /**
4620 16 Nov 17 nicklas 206     Get the number of genotypes called as HOM_REF.
4620 16 Nov 17 nicklas 207   */
4620 16 Nov 17 nicklas 208   public int getHomRefCount()
4620 16 Nov 17 nicklas 209   {
4620 16 Nov 17 nicklas 210     return homRefCount;
4620 16 Nov 17 nicklas 211   }
4620 16 Nov 17 nicklas 212   
4620 16 Nov 17 nicklas 213   /**
4620 16 Nov 17 nicklas 214     Get the number of genotypes called as HOM_ALT.
4620 16 Nov 17 nicklas 215   */
4620 16 Nov 17 nicklas 216   public int getHomAltCount()
4620 16 Nov 17 nicklas 217   {
4620 16 Nov 17 nicklas 218     return homAltCount;
4620 16 Nov 17 nicklas 219   }
4620 16 Nov 17 nicklas 220   
4620 16 Nov 17 nicklas 221   /**
4620 16 Nov 17 nicklas 222     Get the number of genotypes called as HET.
4620 16 Nov 17 nicklas 223   */
4620 16 Nov 17 nicklas 224   public int getHetCount()
4620 16 Nov 17 nicklas 225   {
4620 16 Nov 17 nicklas 226     return hetCount;
4620 16 Nov 17 nicklas 227   }
4620 16 Nov 17 nicklas 228   
4620 16 Nov 17 nicklas 229   /**
4671 07 Feb 18 nicklas 230     Get the percentage of HET genotypes.
4671 07 Feb 18 nicklas 231   */
4671 07 Feb 18 nicklas 232   public float getHetPercentage()
4671 07 Feb 18 nicklas 233   {
7422 14 Nov 23 nicklas 234     return tumorGt.size() > 0 ? 100f * hetCount / tumorGt.size() : 0;
4671 07 Feb 18 nicklas 235   }
4671 07 Feb 18 nicklas 236   
4671 07 Feb 18 nicklas 237   /**
4620 16 Nov 17 nicklas 238     Get the total number of genotypes (not including bad ones).
4620 16 Nov 17 nicklas 239   */
4620 16 Nov 17 nicklas 240   public int getGtCount()
4620 16 Nov 17 nicklas 241   {
7422 14 Nov 23 nicklas 242     return tumorGt.size();
4620 16 Nov 17 nicklas 243   }
4620 16 Nov 17 nicklas 244   
4620 16 Nov 17 nicklas 245   /**
4620 16 Nov 17 nicklas 246     Get the number of 'bad' genotypes.
4620 16 Nov 17 nicklas 247   */
4620 16 Nov 17 nicklas 248   public int getBadCount()
4620 16 Nov 17 nicklas 249   {
4620 16 Nov 17 nicklas 250     return badCount;
4620 16 Nov 17 nicklas 251   }
4620 16 Nov 17 nicklas 252   
4620 16 Nov 17 nicklas 253   /**
7414 10 Nov 23 nicklas 254     Get the number of genotypes that are phased.
7414 10 Nov 23 nicklas 255     @since 4.50
7414 10 Nov 23 nicklas 256   */
7414 10 Nov 23 nicklas 257   public int getPhasedCount()
7414 10 Nov 23 nicklas 258   {
7414 10 Nov 23 nicklas 259     return phasedCount;
7414 10 Nov 23 nicklas 260   }
7414 10 Nov 23 nicklas 261   
7414 10 Nov 23 nicklas 262   /**
4620 16 Nov 17 nicklas 263     Get the sum of all GQ values for all HOM_REF or HOM_ALT genotypes.
4620 16 Nov 17 nicklas 264   */
6513 07 Dec 21 nicklas 265   public float getSumHomGQ()
4620 16 Nov 17 nicklas 266   {
4620 16 Nov 17 nicklas 267     return sumHomGQ;
4620 16 Nov 17 nicklas 268   }
4620 16 Nov 17 nicklas 269   
4620 16 Nov 17 nicklas 270   /**
4620 16 Nov 17 nicklas 271     Get the sum of all GQ values for all HET genotypes.
4620 16 Nov 17 nicklas 272   */
6513 07 Dec 21 nicklas 273   public float getSumHetGQ()
4620 16 Nov 17 nicklas 274   {
4620 16 Nov 17 nicklas 275     return sumHetGQ;
4620 16 Nov 17 nicklas 276   }
4620 16 Nov 17 nicklas 277   
4620 16 Nov 17 nicklas 278   /**
6513 07 Dec 21 nicklas 279     Get the number of genotypes that has a high GQ
6513 07 Dec 21 nicklas 280     value (eg. 99 or xx).
4648 15 Dec 17 nicklas 281   */
6513 07 Dec 21 nicklas 282   public int getHighGQCount()
4648 15 Dec 17 nicklas 283   {
6513 07 Dec 21 nicklas 284     return highGqCount;
4648 15 Dec 17 nicklas 285   }
4648 15 Dec 17 nicklas 286   
4648 15 Dec 17 nicklas 287   /**
6513 07 Dec 21 nicklas 288     Get the number of genotype that has low GQ value
6513 07 Dec 21 nicklas 289     (eg. below 50 or 0.2).
4648 15 Dec 17 nicklas 290   */
6513 07 Dec 21 nicklas 291   public int getLowGQCount()
4648 15 Dec 17 nicklas 292   {
6513 07 Dec 21 nicklas 293     return lowGqCount;
4648 15 Dec 17 nicklas 294   }
4648 15 Dec 17 nicklas 295   
4648 15 Dec 17 nicklas 296   /**
4620 16 Nov 17 nicklas 297     Get all genotypes found in this file.
4620 16 Nov 17 nicklas 298   */
4620 16 Nov 17 nicklas 299   public List<GtData> getGtData()
4620 16 Nov 17 nicklas 300   {
7422 14 Nov 23 nicklas 301     return new ArrayList<>(tumorGt.values());
4620 16 Nov 17 nicklas 302   }
4620 16 Nov 17 nicklas 303
4620 16 Nov 17 nicklas 304   /**
4620 16 Nov 17 nicklas 305     Get the genotype data for the SNP with the given rsId.
4620 16 Nov 17 nicklas 306   */
4620 16 Nov 17 nicklas 307   public GtData getGtData(String id)
4620 16 Nov 17 nicklas 308   {
7422 14 Nov 23 nicklas 309     return tumorGt.get(id);
4620 16 Nov 17 nicklas 310   }
4620 16 Nov 17 nicklas 311
4620 16 Nov 17 nicklas 312   /**
7422 14 Nov 23 nicklas 313     Get the genotype data for paired normal with the given SNP id.
7422 14 Nov 23 nicklas 314     @since 4.50
7422 14 Nov 23 nicklas 315   */
7422 14 Nov 23 nicklas 316   public GtData getNormalGtData(String id)
7422 14 Nov 23 nicklas 317   {
7422 14 Nov 23 nicklas 318     return normalGt.get(id);
7422 14 Nov 23 nicklas 319   }
7422 14 Nov 23 nicklas 320   
7422 14 Nov 23 nicklas 321   /**
4646 14 Dec 17 nicklas 322     Get the id of the File item in BASE this data is coming from.
4646 14 Dec 17 nicklas 323   */
4646 14 Dec 17 nicklas 324   public int getFileId()
4646 14 Dec 17 nicklas 325   {
4646 14 Dec 17 nicklas 326     return fileId;
4646 14 Dec 17 nicklas 327   }
4646 14 Dec 17 nicklas 328   
4646 14 Dec 17 nicklas 329   public void setFileId(int fileId)
4646 14 Dec 17 nicklas 330   {
4646 14 Dec 17 nicklas 331     this.fileId = fileId;
4646 14 Dec 17 nicklas 332   }
4646 14 Dec 17 nicklas 333   
4646 14 Dec 17 nicklas 334   /**
6513 07 Dec 21 nicklas 335     Is the genotype call of low quality?
6513 07 Dec 21 nicklas 336     @since 4.34
6513 07 Dec 21 nicklas 337   */
6513 07 Dec 21 nicklas 338   public boolean isLowQuality(GtData gt)
6513 07 Dec 21 nicklas 339   {
6513 07 Dec 21 nicklas 340     return qualityModel != null && qualityModel.isLow(gt.getGQ());
6513 07 Dec 21 nicklas 341   }
6513 07 Dec 21 nicklas 342   
6513 07 Dec 21 nicklas 343   /**
6513 07 Dec 21 nicklas 344     Is the genotype call of high quality?
6513 07 Dec 21 nicklas 345     @since 4.34
6513 07 Dec 21 nicklas 346   */
6513 07 Dec 21 nicklas 347   public boolean isHighQuality(GtData gt)
6513 07 Dec 21 nicklas 348   {
6513 07 Dec 21 nicklas 349     return qualityModel != null && qualityModel.isHigh(gt.getGQ());
6513 07 Dec 21 nicklas 350   }
6513 07 Dec 21 nicklas 351   
6513 07 Dec 21 nicklas 352   /**
4620 16 Nov 17 nicklas 353     Get all information as a JSON object. A new JSON object 
4620 16 Nov 17 nicklas 354     is created each time this method is called.
4620 16 Nov 17 nicklas 355   */
4620 16 Nov 17 nicklas 356   public JSONObject asJSONObject()
4620 16 Nov 17 nicklas 357   {
6451 22 Oct 21 nicklas 358     return asJSONObject(true);
6451 22 Oct 21 nicklas 359   }
6451 22 Oct 21 nicklas 360   
6451 22 Oct 21 nicklas 361   /**
6451 22 Oct 21 nicklas 362     Get all information as a JSON object. A new JSON object 
6451 22 Oct 21 nicklas 363     is created each time this method is called.
6451 22 Oct 21 nicklas 364     @since 4.34
6451 22 Oct 21 nicklas 365   */
6451 22 Oct 21 nicklas 366   public JSONObject asJSONObject(boolean includeHeaders)
6451 22 Oct 21 nicklas 367   {
4620 16 Nov 17 nicklas 368     JSONObject json = new JSONObject();
6444 19 Oct 21 nicklas 369     json.put("totalCount", getTotalCount());
4620 16 Nov 17 nicklas 370     json.put("homRefCount", getHomRefCount());
4620 16 Nov 17 nicklas 371     json.put("homAltCount", getHomAltCount());
4620 16 Nov 17 nicklas 372     json.put("hetCount", getHetCount());
4620 16 Nov 17 nicklas 373     json.put("gtCount", getGtCount());
4620 16 Nov 17 nicklas 374     json.put("badCount", getBadCount());
7414 10 Nov 23 nicklas 375     json.put("phasedCount", getPhasedCount());
4620 16 Nov 17 nicklas 376     json.put("sumHomGQ", getSumHomGQ());
4620 16 Nov 17 nicklas 377     json.put("sumHetGQ", getSumHetGQ());
6513 07 Dec 21 nicklas 378     if (getQualityModel() != null) json.put("qualityModel", getQualityModel().name());
6513 07 Dec 21 nicklas 379     json.put("highGQCount", getHighGQCount());
6513 07 Dec 21 nicklas 380     json.put("lowGQCount", getLowGQCount());
4646 14 Dec 17 nicklas 381     json.put("fileId", getFileId());
6452 22 Oct 21 nicklas 382     if (includeHeaders)
6451 22 Oct 21 nicklas 383     {
6451 22 Oct 21 nicklas 384       json.put("infoHeaders", infoHeaders);
6451 22 Oct 21 nicklas 385       json.put("formatHeaders", formatHeaders);
6451 22 Oct 21 nicklas 386     }
4620 16 Nov 17 nicklas 387     return json;
4620 16 Nov 17 nicklas 388   }
6451 22 Oct 21 nicklas 389
4620 16 Nov 17 nicklas 390 }

Generated on Tue Nov 21 03:35:56 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)