extensions/net.sf.basedb.reggie/trunk/src/net/sf/basedb/reggie/vcf/VcfParser.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
4620 16 Nov 17 nicklas 3 import java.io.IOException;
4620 16 Nov 17 nicklas 4 import java.io.InputStream;
4620 16 Nov 17 nicklas 5 import java.util.ArrayList;
5692 29 Oct 19 nicklas 6 import java.util.Date;
4620 16 Nov 17 nicklas 7 import java.util.LinkedHashMap;
4620 16 Nov 17 nicklas 8 import java.util.List;
4620 16 Nov 17 nicklas 9 import java.util.Map;
4620 16 Nov 17 nicklas 10 import java.util.regex.Pattern;
4620 16 Nov 17 nicklas 11
4620 16 Nov 17 nicklas 12 import net.sf.basedb.reggie.vcf.FormatFactory.Format;
4620 16 Nov 17 nicklas 13 import net.sf.basedb.util.Values;
4620 16 Nov 17 nicklas 14 import net.sf.basedb.util.parser.FlatFileParser;
5692 29 Oct 19 nicklas 15 import net.sf.basedb.util.parser.FlatFileParser.Data;
6450 21 Oct 21 nicklas 16 import net.sf.basedb.util.parser.FlatFileParser.Line;
4620 16 Nov 17 nicklas 17 import net.sf.basedb.util.parser.FlatFileParser.LineType;
4620 16 Nov 17 nicklas 18 import net.sf.basedb.util.parser.Mapper;
4620 16 Nov 17 nicklas 19
4620 16 Nov 17 nicklas 20 /**
4620 16 Nov 17 nicklas 21   Parser implementation for parsing VCF files. A single
4620 16 Nov 17 nicklas 22   parser can be used to parse one or more files.
4620 16 Nov 17 nicklas 23
4620 16 Nov 17 nicklas 24   @author nicklas
4620 16 Nov 17 nicklas 25   @since 4.13
4620 16 Nov 17 nicklas 26 */
4620 16 Nov 17 nicklas 27 public class VcfParser 
4620 16 Nov 17 nicklas 28 {
4620 16 Nov 17 nicklas 29
4620 16 Nov 17 nicklas 30   
6445 20 Oct 21 nicklas 31   private final Map<String, SnpData> snpById;
6445 20 Oct 21 nicklas 32   private final Map<String, SnpData> snpByLoc;
5692 29 Oct 19 nicklas 33   private boolean useLineNoAsId;
5729 18 Nov 19 nicklas 34   private InfoFactory infoFactory;
6444 19 Oct 21 nicklas 35   private boolean hasRef;
4620 16 Nov 17 nicklas 36   
4620 16 Nov 17 nicklas 37   /**
4620 16 Nov 17 nicklas 38     Creates a new parser instance.
4620 16 Nov 17 nicklas 39   */
4620 16 Nov 17 nicklas 40   public VcfParser()
4620 16 Nov 17 nicklas 41   {
6445 20 Oct 21 nicklas 42     this.snpById = new LinkedHashMap<>();
6445 20 Oct 21 nicklas 43     this.snpByLoc = new LinkedHashMap<>();
4620 16 Nov 17 nicklas 44   }
4620 16 Nov 17 nicklas 45   
4620 16 Nov 17 nicklas 46   /**
6444 19 Oct 21 nicklas 47     Parse a reference VCF used for defining all variants of interest. 
6444 19 Oct 21 nicklas 48     This may be called multiple times to add more variants. After this method
6444 19 Oct 21 nicklas 49     has been called, the parse() method will only return data for the 
6444 19 Oct 21 nicklas 50     defined variants.
6444 19 Oct 21 nicklas 51     @since 4.34
6444 19 Oct 21 nicklas 52   */
6444 19 Oct 21 nicklas 53   public void parseRef(InputStream in, String name)
6444 19 Oct 21 nicklas 54     throws IOException
6444 19 Oct 21 nicklas 55   {
6444 19 Oct 21 nicklas 56     parseInternal(in, name, null, true);
6444 19 Oct 21 nicklas 57     hasRef = true;
6444 19 Oct 21 nicklas 58   }
6444 19 Oct 21 nicklas 59   
6444 19 Oct 21 nicklas 60   /**
5729 18 Nov 19 nicklas 61     Set this to TRUE to use line number as the id. Useful if there
5692 29 Oct 19 nicklas 62     are no ID values in the file, but note that this means that 
5692 29 Oct 19 nicklas 63     different VCF files can't be compared.
5729 18 Nov 19 nicklas 64     @since 4.24
5692 29 Oct 19 nicklas 65   */
5692 29 Oct 19 nicklas 66   public void setUseLineNoAsId(boolean useLineNoAsId)
5692 29 Oct 19 nicklas 67   {
5692 29 Oct 19 nicklas 68     this.useLineNoAsId = useLineNoAsId;
5692 29 Oct 19 nicklas 69   }
5692 29 Oct 19 nicklas 70   
5692 29 Oct 19 nicklas 71   /**
5729 18 Nov 19 nicklas 72     Set the factory to use for extracting information from
5729 18 Nov 19 nicklas 73     the INFO column. If not set, the INFO column is not parsed.
5729 18 Nov 19 nicklas 74     @since 4.24
5729 18 Nov 19 nicklas 75   */
5729 18 Nov 19 nicklas 76   public void setInfoFactory(InfoFactory infoFactory)
5729 18 Nov 19 nicklas 77   {
5729 18 Nov 19 nicklas 78     this.infoFactory = infoFactory;
5729 18 Nov 19 nicklas 79   }
5729 18 Nov 19 nicklas 80   
5729 18 Nov 19 nicklas 81   /**
4620 16 Nov 17 nicklas 82     Parse the given input stream as a VCF file. It is expected
4620 16 Nov 17 nicklas 83     that the file only contains data for a single sample.
4620 16 Nov 17 nicklas 84   */
4620 16 Nov 17 nicklas 85   public VcfData parse(InputStream in, String fileName)
5040 19 Oct 18 nicklas 86       throws IOException
5040 19 Oct 18 nicklas 87   {
5040 19 Oct 18 nicklas 88     return parse(in, fileName, null);
5040 19 Oct 18 nicklas 89   }
6444 19 Oct 21 nicklas 90
5040 19 Oct 18 nicklas 91   /**
5040 19 Oct 18 nicklas 92     Parse the given input stream as a VCF file. It is expected
5040 19 Oct 18 nicklas 93     that the file only contains data for a single sample.
5040 19 Oct 18 nicklas 94   */
5040 19 Oct 18 nicklas 95   public VcfData parse(InputStream in, String fileName, VcfDataFilter filter)
4620 16 Nov 17 nicklas 96     throws IOException
4620 16 Nov 17 nicklas 97   {
6444 19 Oct 21 nicklas 98     return parseInternal(in, fileName, filter, false);
6444 19 Oct 21 nicklas 99   }
6444 19 Oct 21 nicklas 100
6444 19 Oct 21 nicklas 101   
6444 19 Oct 21 nicklas 102   /**
6444 19 Oct 21 nicklas 103     Parse the given input stream as a VCF file. It is expected
6444 19 Oct 21 nicklas 104     that the file only contains data for a single sample.
6444 19 Oct 21 nicklas 105   */
6444 19 Oct 21 nicklas 106   private VcfData parseInternal(InputStream in, String fileName, VcfDataFilter filter, boolean isRef)
6444 19 Oct 21 nicklas 107     throws IOException
6444 19 Oct 21 nicklas 108   {
5040 19 Oct 18 nicklas 109     if (filter == null) filter = new VcfDataFilter(){}; // Pass all
4620 16 Nov 17 nicklas 110     FlatFileParser ffp = createFfp();
4620 16 Nov 17 nicklas 111     ffp.setInputStream(in, "UTF-8");
4620 16 Nov 17 nicklas 112     
4620 16 Nov 17 nicklas 113     LineType headerLine = ffp.parseHeaders();
4620 16 Nov 17 nicklas 114     int lineNo = ffp.getParsedLines();
4620 16 Nov 17 nicklas 115     if (headerLine != LineType.DATA_HEADER)
4620 16 Nov 17 nicklas 116     {
4620 16 Nov 17 nicklas 117       throw new IOException("File '" + fileName + "' line " + lineNo + ": Could not find header line starting with '#CHROM{tab}POS...'");
4620 16 Nov 17 nicklas 118     }
7422 14 Nov 23 nicklas 119     List<String> headers = ffp.getColumnHeaders();
7422 14 Nov 23 nicklas 120
7422 14 Nov 23 nicklas 121     boolean isPaired = headers.size() > 10;
7422 14 Nov 23 nicklas 122     VcfData vcf = new VcfData(isPaired);
6450 21 Oct 21 nicklas 123     for (int lNo = 0; lNo < ffp.getLineCount(); lNo++)
5729 18 Nov 19 nicklas 124     {
6450 21 Oct 21 nicklas 125       Line l = ffp.getLine(lNo);
6450 21 Oct 21 nicklas 126       if (l.type() == LineType.HEADER)
5729 18 Nov 19 nicklas 127       {
6450 21 Oct 21 nicklas 128         if (l.line().contains("INFO"))
6450 21 Oct 21 nicklas 129         {
6450 21 Oct 21 nicklas 130           vcf.addInfoHeader(l.name(), l.value());
6450 21 Oct 21 nicklas 131           if (infoFactory != null) infoFactory.addInfoHeader(l.name(), l.value());
6450 21 Oct 21 nicklas 132         }
6450 21 Oct 21 nicklas 133         else if (l.line().contains("FORMAT"))
6450 21 Oct 21 nicklas 134         {
6450 21 Oct 21 nicklas 135           vcf.addFormatHeader(l.name(), l.value());
6450 21 Oct 21 nicklas 136         }
5729 18 Nov 19 nicklas 137       }
5729 18 Nov 19 nicklas 138     }
5729 18 Nov 19 nicklas 139     
4620 16 Nov 17 nicklas 140     Mapper chromMapper = ffp.getMapper("\\#CHROM\\");
4620 16 Nov 17 nicklas 141     Mapper posMapper = ffp.getMapper("\\POS\\");
5692 29 Oct 19 nicklas 142     Mapper idMapper = useLineNoAsId ? new LineNoMapper() : ffp.getMapper("\\ID\\");
4620 16 Nov 17 nicklas 143     Mapper refMapper = ffp.getMapper("\\REF\\");
4620 16 Nov 17 nicklas 144     Mapper altMapper = ffp.getMapper("\\ALT\\");
7413 10 Nov 23 nicklas 145     Mapper filterMapper = ffp.getMapper("\\FILTER\\");
6444 19 Oct 21 nicklas 146     Mapper infoMapper = infoFactory == null ? null : ffp.getMapper("\\INFO\\");
6444 19 Oct 21 nicklas 147     Mapper formatMapper = isRef ? null : ffp.getMapper("\\FORMAT\\");
7422 14 Nov 23 nicklas 148     Mapper tumorMapper = isRef ? null : ffp.getMapper("\\" + (headers.size()-1)+"\\");
7422 14 Nov 23 nicklas 149     Mapper normalMapper = isRef || !isPaired ? null : ffp.getMapper("\\" + (headers.size()-2)+"\\");
4620 16 Nov 17 nicklas 150         
4620 16 Nov 17 nicklas 151     FormatFactory ff = new FormatFactory();
4620 16 Nov 17 nicklas 152     
6444 19 Oct 21 nicklas 153     int numSnp = 0;
4620 16 Nov 17 nicklas 154     while (ffp.hasMoreData())
4620 16 Nov 17 nicklas 155     {
4620 16 Nov 17 nicklas 156       FlatFileParser.Data line = ffp.nextData();
6444 19 Oct 21 nicklas 157       numSnp++;
6445 20 Oct 21 nicklas 158       
6445 20 Oct 21 nicklas 159       String chrom = chromMapper.getString(line);
6445 20 Oct 21 nicklas 160       if (!chrom.startsWith("chr")) chrom = "chr"+chrom;
6445 20 Oct 21 nicklas 161       String pos = posMapper.getString(line);
6445 20 Oct 21 nicklas 162       String ref = refMapper.getString(line);
6445 20 Oct 21 nicklas 163       String alt = altMapper.getString(line);
6445 20 Oct 21 nicklas 164       String loc = chrom+":"+pos+":"+ref+">"+alt;
7413 10 Nov 23 nicklas 165       String flt = filterMapper.getString(line);
6445 20 Oct 21 nicklas 166       
6445 20 Oct 21 nicklas 167       String id = idMapper.getString(line);
6445 20 Oct 21 nicklas 168       if (id == null || ".".equals(id)) id = loc;
4620 16 Nov 17 nicklas 169
4620 16 Nov 17 nicklas 170       // Add the SNP information if we have not seen this one before
6445 20 Oct 21 nicklas 171       SnpData snp = snpById.get(id);
6455 25 Oct 21 nicklas 172       boolean swapRefAlt = false;
6455 25 Oct 21 nicklas 173       if (snp != null && !snp.hasRefAlt(ref, alt)) 
6455 25 Oct 21 nicklas 174       {
6455 25 Oct 21 nicklas 175         if (snp.hasRefAlt(alt, ref))
6455 25 Oct 21 nicklas 176         {
6455 25 Oct 21 nicklas 177           swapRefAlt = true;
6455 25 Oct 21 nicklas 178         }
6455 25 Oct 21 nicklas 179         else
6455 25 Oct 21 nicklas 180         {
6455 25 Oct 21 nicklas 181           snp = null;
6455 25 Oct 21 nicklas 182         }
6455 25 Oct 21 nicklas 183       }
6445 20 Oct 21 nicklas 184       if (snp == null) snp = snpByLoc.get(loc);
5040 19 Oct 18 nicklas 185       if (snp == null)
4620 16 Nov 17 nicklas 186       {
6444 19 Oct 21 nicklas 187         if (hasRef && !isRef) continue; // Ignore this SNP and continue with next line
6444 19 Oct 21 nicklas 188         
5040 19 Oct 18 nicklas 189         snp = new SnpData(id);
7413 10 Nov 23 nicklas 190         snp.setChromosome(chrom);
7413 10 Nov 23 nicklas 191         snp.setPosition(Values.getLong(pos));
7413 10 Nov 23 nicklas 192         snp.setRef(ref);
7413 10 Nov 23 nicklas 193         snp.setAlt(alt);
7413 10 Nov 23 nicklas 194         snp.setFilter(flt);
5729 18 Nov 19 nicklas 195         if (infoFactory != null)
5729 18 Nov 19 nicklas 196         {
5729 18 Nov 19 nicklas 197           snp.setInfo(infoFactory.getInfo(infoMapper.getString(line), snp));
5729 18 Nov 19 nicklas 198         }
7413 10 Nov 23 nicklas 199         if (filter.acceptSnp(snp))
7413 10 Nov 23 nicklas 200         {
7413 10 Nov 23 nicklas 201           snpById.put(id, snp);
7413 10 Nov 23 nicklas 202           snpByLoc.put(loc, snp);
7413 10 Nov 23 nicklas 203         }
4620 16 Nov 17 nicklas 204       }
6455 25 Oct 21 nicklas 205
6444 19 Oct 21 nicklas 206       if (!isRef)
6444 19 Oct 21 nicklas 207       {
6444 19 Oct 21 nicklas 208         String format = formatMapper.getString(line);
7422 14 Nov 23 nicklas 209         String tumorData = tumorMapper.getString(line);  
7422 14 Nov 23 nicklas 210         if (format == null || tumorData == null) continue;
6444 19 Oct 21 nicklas 211         
6455 25 Oct 21 nicklas 212         Format fmt = ff.getFormat(format, swapRefAlt);
7422 14 Nov 23 nicklas 213         GtData gt = fmt.extractTo(tumorData, new GtData(snp.getId()));
7422 14 Nov 23 nicklas 214         GtData normalGt = normalMapper != null ? fmt.extractTo(normalMapper.getString(line), new GtData(snp.getId())) : null;
7413 10 Nov 23 nicklas 215         gt.setFilter(flt);
7422 14 Nov 23 nicklas 216         if (filter.accept(gt, snp)) vcf.addGt(normalGt, gt);
6444 19 Oct 21 nicklas 217       }
4620 16 Nov 17 nicklas 218     }
6444 19 Oct 21 nicklas 219     vcf.setTotalCount(numSnp);
4620 16 Nov 17 nicklas 220     return vcf;
4620 16 Nov 17 nicklas 221   }
4620 16 Nov 17 nicklas 222   
4620 16 Nov 17 nicklas 223   /**
6387 15 Sep 21 nicklas 224     Get the number of different SNP:s that has been parsed.
6389 15 Sep 21 nicklas 225     @since 4.32
6387 15 Sep 21 nicklas 226   */
6387 15 Sep 21 nicklas 227   public int getSnpCount()
6387 15 Sep 21 nicklas 228   {
6445 20 Oct 21 nicklas 229     return snpById.size();
6387 15 Sep 21 nicklas 230   }
6387 15 Sep 21 nicklas 231   
6387 15 Sep 21 nicklas 232   /**
4620 16 Nov 17 nicklas 233     Get all SNP data we have seen so far.
4620 16 Nov 17 nicklas 234   */
4620 16 Nov 17 nicklas 235   public List<SnpData> getSnpData()
4620 16 Nov 17 nicklas 236   {
6445 20 Oct 21 nicklas 237     return new ArrayList<>(snpById.values());
4620 16 Nov 17 nicklas 238   }
4620 16 Nov 17 nicklas 239   
4620 16 Nov 17 nicklas 240   /**
4620 16 Nov 17 nicklas 241     Get data for the SNP with the given id.
4620 16 Nov 17 nicklas 242   */
4620 16 Nov 17 nicklas 243   public SnpData getSnp(String id)
4620 16 Nov 17 nicklas 244   {
6445 20 Oct 21 nicklas 245     return snpById.get(id);
4620 16 Nov 17 nicklas 246   }
4620 16 Nov 17 nicklas 247   
4620 16 Nov 17 nicklas 248   /**
4620 16 Nov 17 nicklas 249     We look for the header line starting with: #CHROM{tab}POS
4620 16 Nov 17 nicklas 250     and then split all data lines on {tab}. Lines starting with
4620 16 Nov 17 nicklas 251     two # are ignored.
4620 16 Nov 17 nicklas 252   */
4620 16 Nov 17 nicklas 253   private FlatFileParser createFfp()
4620 16 Nov 17 nicklas 254   {
4620 16 Nov 17 nicklas 255     FlatFileParser ffp = new FlatFileParser();
4620 16 Nov 17 nicklas 256     ffp.setDataHeaderRegexp(Pattern.compile("#CHROM\\tPOS.*"));
4620 16 Nov 17 nicklas 257     ffp.setDataSplitterRegexp(Pattern.compile("\\t"));
6450 21 Oct 21 nicklas 258     ffp.setHeaderRegexp(Pattern.compile("##(?:INFO|FORMAT)\\=\\<ID\\=(\\w+),(.*)\\>"));
4620 16 Nov 17 nicklas 259     ffp.setIgnoreRegexp(Pattern.compile("##.*"));
4620 16 Nov 17 nicklas 260     return ffp;
4620 16 Nov 17 nicklas 261   }
4620 16 Nov 17 nicklas 262   
5692 29 Oct 19 nicklas 263   static class LineNoMapper
5692 29 Oct 19 nicklas 264     implements Mapper
5692 29 Oct 19 nicklas 265   {
5692 29 Oct 19 nicklas 266     @Override
5692 29 Oct 19 nicklas 267     @Deprecated
5692 29 Oct 19 nicklas 268     public String getValue(Data line) 
5692 29 Oct 19 nicklas 269     {
5692 29 Oct 19 nicklas 270       return getString(line);
5692 29 Oct 19 nicklas 271     }
5692 29 Oct 19 nicklas 272     
5692 29 Oct 19 nicklas 273     @Override
5692 29 Oct 19 nicklas 274     public String getString(Data line) 
5692 29 Oct 19 nicklas 275     {
5692 29 Oct 19 nicklas 276       return Integer.toString(line.dataLineNo());
5692 29 Oct 19 nicklas 277     }
5692 29 Oct 19 nicklas 278     
5692 29 Oct 19 nicklas 279     @Override
5692 29 Oct 19 nicklas 280     public Long getLong(Data line)
5692 29 Oct 19 nicklas 281     {
5692 29 Oct 19 nicklas 282       return (long)line.dataLineNo();
5692 29 Oct 19 nicklas 283     }
5692 29 Oct 19 nicklas 284     
5692 29 Oct 19 nicklas 285     @Override
5692 29 Oct 19 nicklas 286     public Integer getInt(Data line) 
5692 29 Oct 19 nicklas 287     {
5692 29 Oct 19 nicklas 288       return line.dataLineNo();
5692 29 Oct 19 nicklas 289     }
5692 29 Oct 19 nicklas 290     
5692 29 Oct 19 nicklas 291     @Override
5692 29 Oct 19 nicklas 292     public Float getFloat(Data line) 
5692 29 Oct 19 nicklas 293     {
5692 29 Oct 19 nicklas 294       return (float)line.dataLineNo();
5692 29 Oct 19 nicklas 295     }
5692 29 Oct 19 nicklas 296     
5692 29 Oct 19 nicklas 297     @Override
5692 29 Oct 19 nicklas 298     public Double getDouble(Data line)
5692 29 Oct 19 nicklas 299     {
5692 29 Oct 19 nicklas 300       return (double)line.dataLineNo();
5692 29 Oct 19 nicklas 301     }
5692 29 Oct 19 nicklas 302     
5692 29 Oct 19 nicklas 303     @Override
5692 29 Oct 19 nicklas 304     public Date getDate(Data line)
5692 29 Oct 19 nicklas 305     {
5692 29 Oct 19 nicklas 306       return null;
5692 29 Oct 19 nicklas 307     }
5692 29 Oct 19 nicklas 308     
5692 29 Oct 19 nicklas 309   }
4620 16 Nov 17 nicklas 310 }

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