/* * Jalview - A Sequence Alignment Editor and Viewer ($$Version-Rel$$) * Copyright (C) $$Year-Rel$$ The Jalview Authors * * This file is part of Jalview. * * Jalview is free software: you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation, either version 3 * of the License, or (at your option) any later version. * * Jalview is distributed in the hope that it will be useful, but * WITHOUT ANY WARRANTY; without even the implied warranty * of MERCHANTABILITY or FITNESS FOR A PARTICULAR * PURPOSE. See the GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with Jalview. If not, see . * The Jalview Authors are detailed in the 'AUTHORS' file. */ package jalview.analysis; import jalview.datamodel.AlignedCodonFrame; import jalview.datamodel.Alignment; import jalview.datamodel.AlignmentI; import jalview.datamodel.DBRefEntry; import jalview.datamodel.DBRefSource; import jalview.datamodel.Mapping; import jalview.datamodel.Sequence; import jalview.datamodel.SequenceFeature; import jalview.datamodel.SequenceI; import jalview.io.gff.SequenceOntology; import jalview.schemes.ResidueProperties; import jalview.util.DBRefUtils; import jalview.util.MapList; import jalview.util.MappingUtils; import jalview.util.StringUtils; import jalview.ws.SequenceFetcher; import jalview.ws.seqfetcher.ASequenceFetcher; import java.util.ArrayList; import java.util.Collections; import java.util.LinkedHashMap; import java.util.List; import java.util.Map.Entry; import java.util.Vector; /** * Functions for cross-referencing sequence databases. user must first specify * if cross-referencing from protein or dna (set dna==true) * * @author JimP * */ public class CrossRef { /** * Select just the DNA or protein references for a protein or dna sequence * * @param fromDna * if true, select references from DNA (i.e. Protein databases), else * DNA database references * @param refs * a set of references to select from * @return */ public static DBRefEntry[] findXDbRefs(boolean fromDna, DBRefEntry[] refs) { return DBRefUtils.selectRefs(refs, fromDna ? DBRefSource.PROTEINDBS : DBRefSource.DNACODINGDBS); // could attempt to find other cross // refs here - ie PDB xrefs // (not dna, not protein seq) } /** * @param dna * true if seqs are DNA seqs * @param seqs * @return a list of sequence database cross reference source types */ public static String[] findSequenceXrefTypes(boolean dna, SequenceI[] seqs) { return findSequenceXrefTypes(dna, seqs, null); } /** * Indirect references are references from other sequences from the dataset to * any of the direct DBRefEntrys on the given sequences. * * @param dna * true if seqs are DNA seqs * @param seqs * @return a list of sequence database cross reference source types */ public static String[] findSequenceXrefTypes(boolean dna, SequenceI[] seqs, AlignmentI dataset) { String[] dbrefs = null; List refs = new ArrayList(); for (SequenceI seq : seqs) { if (seq != null) { SequenceI dss = seq; while (dss.getDatasetSequence() != null) { dss = dss.getDatasetSequence(); } DBRefEntry[] rfs = findXDbRefs(dna, dss.getDBRef()); if (rfs != null) { for (DBRefEntry ref : rfs) { if (!refs.contains(ref.getSource())) { refs.add(ref.getSource()); } } } if (dataset != null) { // search for references to this sequence's direct references. DBRefEntry[] lrfs = CrossRef.findXDbRefs(!dna, seq.getDBRef()); List rseqs = new ArrayList(); CrossRef.searchDatasetXrefs(seq, !dna, lrfs, dataset, rseqs, null); // don't need to specify codon frame for mapping here for (SequenceI rs : rseqs) { DBRefEntry[] xrs = findXDbRefs(dna, rs.getDBRef()); if (xrs != null) { for (DBRefEntry ref : xrs) { if (!refs.contains(ref.getSource())) { refs.add(ref.getSource()); } } } // looks like copy and paste - change rfs to xrs? // for (int r = 0; rfs != null && r < rfs.length; r++) // { // if (!refs.contains(rfs[r].getSource())) // { // refs.add(rfs[r].getSource()); // } // } } } } } if (refs.size() > 0) { dbrefs = new String[refs.size()]; refs.toArray(dbrefs); } return dbrefs; } public static boolean hasCdnaMap(SequenceI[] seqs) { // TODO unused - remove? String[] reftypes = findSequenceXrefTypes(false, seqs); for (int s = 0; s < reftypes.length; s++) { if (reftypes.equals(DBRefSource.EMBLCDS)) { return true; // no map } } return false; } public static SequenceI[] getCdnaMap(SequenceI[] seqs) { // TODO unused - remove? Vector cseqs = new Vector(); for (int s = 0; s < seqs.length; s++) { DBRefEntry[] cdna = findXDbRefs(true, seqs[s].getDBRef()); for (int c = 0; c < cdna.length; c++) { if (cdna[c].getSource().equals(DBRefSource.EMBLCDS)) { System.err .println("TODO: unimplemented sequence retrieval for coding region sequence."); // TODO: retrieve CDS dataset sequences // need global dataset sequence retriever/resolver to reuse refs // and construct Mapping entry. // insert gaps in CDS according to peptide gaps. // add gapped sequence to cseqs } } } if (cseqs.size() > 0) { SequenceI[] rsqs = new SequenceI[cseqs.size()]; cseqs.copyInto(rsqs); return rsqs; } return null; } /** * * @param dna * @param seqs * @return */ public static Alignment findXrefSequences(SequenceI[] seqs, boolean dna, String source) { return findXrefSequences(seqs, dna, source, null); } /** * * @param seqs * sequences whose xrefs are being retrieved * @param dna * true if sequences are nucleotide * @param source * @param dataset * alignment to search for product sequences. * @return products (as dataset sequences) */ public static Alignment findXrefSequences(SequenceI[] seqs, boolean dna, String source, AlignmentI dataset) { List rseqs = new ArrayList(); AlignedCodonFrame cf = new AlignedCodonFrame(); for (SequenceI seq : seqs) { SequenceI dss = seq; while (dss.getDatasetSequence() != null) { dss = dss.getDatasetSequence(); } boolean found = false; DBRefEntry[] xrfs = CrossRef.findXDbRefs(dna, dss.getDBRef()); if ((xrfs == null || xrfs.length == 0) && dataset != null) { System.out.println("Attempting to find ds Xrefs refs."); // FIXME should be dss not seq here? DBRefEntry[] lrfs = CrossRef.findXDbRefs(!dna, seq.getDBRef()); // less ambiguous would be a 'find primary dbRefEntry' method. // filter for desired source xref here found = CrossRef.searchDatasetXrefs(dss, !dna, lrfs, dataset, rseqs, cf); } for (int r = 0; xrfs != null && r < xrfs.length; r++) { DBRefEntry xref = xrfs[r]; if (source != null && !source.equals(xref.getSource())) { continue; } if (xref.hasMap()) { if (xref.getMap().getTo() != null) { SequenceI rsq = new Sequence(xref.getMap().getTo()); rseqs.add(rsq); if (xref.getMap().getMap().getFromRatio() != xref .getMap().getMap().getToRatio()) { // get sense of map correct for adding to product alignment. if (dna) { // map is from dna seq to a protein product cf.addMap(dss, rsq, xref.getMap().getMap()); } else { // map should be from protein seq to its coding dna cf.addMap(rsq, dss, xref.getMap().getMap().getInverse()); } /* * compute peptide variants from dna variants */ rsq.createDatasetSequence(); computeProteinVariants(seq, rsq, xref.getMap().getMap()); } found = true; } } if (!found) { // do a bit more work - search for sequences with references matching // xrefs on this sequence. if (dataset != null) { found |= searchDataset(dss, xref, dataset, rseqs, cf); // ,false,!dna); if (found) { xrfs[r] = null; // we've recovered seqs for this one. } } } } if (!found) { if (xrfs != null && xrfs.length > 0) { // Try and get the sequence reference... /* * Ideal world - we ask for a sequence fetcher implementation here if * (jalview.io.RunTimeEnvironment.getSequenceFetcher()) ( */ ASequenceFetcher sftch = new SequenceFetcher(); SequenceI[] retrieved = null; int l = xrfs.length; for (int r = 0; r < xrfs.length; r++) { // filter out any irrelevant or irretrievable references if (xrfs[r] == null || ((source != null && !source.equals(xrfs[r] .getSource())) || !sftch.isFetchable(xrfs[r] .getSource()))) { l--; xrfs[r] = null; } } if (l > 0) { System.out .println("Attempting to retrieve cross referenced sequences."); DBRefEntry[] t = new DBRefEntry[l]; l = 0; for (int r = 0; r < xrfs.length; r++) { if (xrfs[r] != null) { t[l++] = xrfs[r]; } } xrfs = t; try { retrieved = sftch.getSequences(xrfs, !dna); // problem here is we don't know which of xrfs resulted in which // retrieved element } catch (Exception e) { System.err .println("Problem whilst retrieving cross references for Sequence : " + seq.getName()); e.printStackTrace(); } if (retrieved != null) { for (int rs = 0; rs < retrieved.length; rs++) { // TODO: examine each sequence for 'redundancy' DBRefEntry[] dbr = retrieved[rs].getDBRef(); if (dbr != null && dbr.length > 0) { for (int di = 0; di < dbr.length; di++) { // find any entry where we should put in the sequence being // cross-referenced into the map Mapping map = dbr[di].getMap(); if (map != null) { if (map.getTo() != null && map.getMap() != null) { // should search the local dataset to find any existing // candidates for To ! try { // compare ms with dss and replace with dss in mapping // if map is congruent SequenceI ms = map.getTo(); int sf = map.getMap().getToLowest(); int st = map.getMap().getToHighest(); SequenceI mappedrg = ms.getSubSequence(sf, st); SequenceI loc = dss.getSubSequence(sf, st); if (mappedrg.getLength() > 0 && mappedrg.getSequenceAsString().equals( loc.getSequenceAsString())) { System.err .println("Mapping updated for retrieved crossreference"); // method to update all refs of existing To on // retrieved sequence with dss and merge any props // on To onto dss. map.setTo(dss); } } catch (Exception e) { System.err .println("Exception when consolidating Mapped sequence set..."); e.printStackTrace(System.err); } } } } } retrieved[rs].updatePDBIds(); rseqs.add(retrieved[rs]); } } } } } } Alignment ral = null; if (rseqs.size() > 0) { SequenceI[] rsqs = new SequenceI[rseqs.size()]; rseqs.toArray(rsqs); ral = new Alignment(rsqs); if (cf != null && !cf.isEmpty()) { ral.addCodonFrame(cf); } } return ral; } /** * find references to lrfs in the cross-reference set of each sequence in * dataset (that is not equal to sequenceI) Identifies matching DBRefEntry * based on source and accession string only - Map and Version are nulled. * * @param sequenceI * @param lrfs * @param dataset * @param rseqs * @return true if matches were found. */ private static boolean searchDatasetXrefs(SequenceI sequenceI, boolean dna, DBRefEntry[] lrfs, AlignmentI dataset, List rseqs, AlignedCodonFrame cf) { boolean found = false; if (lrfs == null) { return false; } for (int i = 0; i < lrfs.length; i++) { DBRefEntry xref = new DBRefEntry(lrfs[i]); // add in wildcards xref.setVersion(null); xref.setMap(null); found = searchDataset(sequenceI, xref, dataset, rseqs, cf, false, dna); } return found; } /** * search a given sequence dataset for references matching cross-references to * the given sequence * * @param sequenceI * @param xrf * @param dataset * @param rseqs * set of unique sequences * @param cf * @return true if one or more unique sequences were found and added */ public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf, AlignmentI dataset, List rseqs, AlignedCodonFrame cf) { return searchDataset(sequenceI, xrf, dataset, rseqs, cf, true, false); } /** * TODO: generalise to different protein classifications Search dataset for * DBRefEntrys matching the given one (xrf) and add the associated sequence to * rseq. * * @param sequenceI * @param xrf * @param dataset * @param rseqs * @param direct * - search all references or only subset * @param dna * search dna or protein xrefs (if direct=false) * @return true if relationship found and sequence added. */ public static boolean searchDataset(SequenceI sequenceI, DBRefEntry xrf, AlignmentI dataset, List rseqs, AlignedCodonFrame cf, boolean direct, boolean dna) { boolean found = false; SequenceI[] typer = new SequenceI[1]; if (dataset == null) { return false; } if (dataset.getSequences() == null) { System.err.println("Empty dataset sequence set - NO VECTOR"); return false; } List ds; synchronized (ds = dataset.getSequences()) { for (SequenceI nxt : ds) { if (nxt != null) { if (nxt.getDatasetSequence() != null) { System.err .println("Implementation warning: getProducts passed a dataset alignment without dataset sequences in it!"); } if (nxt != sequenceI && nxt != sequenceI.getDatasetSequence()) { // check if this is the correct sequence type { typer[0] = nxt; boolean isDna = jalview.util.Comparison.isNucleotide(typer); if ((direct && isDna == dna) || (!direct && isDna != dna)) { // skip this sequence because it is same molecule type continue; } } // look for direct or indirect references in common DBRefEntry[] poss = nxt.getDBRef(), cands = null; if (direct) { cands = jalview.util.DBRefUtils.searchRefs(poss, xrf); } else { poss = CrossRef.findXDbRefs(dna, poss); // cands = jalview.util.DBRefUtils.searchRefs(poss, xrf); } if (cands != null) { if (!rseqs.contains(nxt)) { rseqs.add(nxt); boolean foundmap = cf != null; // don't search if we aren't given a codon map object for (int r = 0; foundmap && r < cands.length; r++) { if (cands[r].hasMap()) { if (cands[r].getMap().getTo() != null && cands[r].getMap().getMap().getFromRatio() != cands[r] .getMap().getMap().getToRatio()) { foundmap = true; // get sense of map correct for adding to product // alignment. if (dna) { // map is from dna seq to a protein product cf.addMap(sequenceI, nxt, cands[r].getMap() .getMap()); } else { // map should be from protein seq to its coding dna cf.addMap(nxt, sequenceI, cands[r].getMap() .getMap().getInverse()); } } } } // TODO: add mapping between sequences if necessary found = true; } } } } } } return found; } /** * Computes variants in peptide product generated by variants in dna, and adds * them as sequence_variant features on the protein sequence. Returns the * number of variant features added. * * @param dnaSeq * @param peptide * @param dnaToProtein */ protected static int computeProteinVariants(SequenceI dnaSeq, SequenceI peptide, MapList dnaToProtein) { /* * map from peptide position to all variant features of the codon for it * LinkedHashMap ensures we add the peptide features in sequence order */ LinkedHashMap variants = new LinkedHashMap(); SequenceOntology so = SequenceOntology.getInstance(); SequenceFeature[] dnaFeatures = dnaSeq.getSequenceFeatures(); if (dnaFeatures == null) { return 0; } int[] lastCodon = null; int lastPeptidePostion = 0; /* * build a map of codon variations for peptides */ for (SequenceFeature sf : dnaFeatures) { int dnaCol = sf.getBegin(); if (dnaCol != sf.getEnd()) { // not handling multi-locus variant features continue; } if (so.isSequenceVariant(sf.getType())) { int[] mapsTo = dnaToProtein.locateInTo(dnaCol, dnaCol); if (mapsTo == null) { // feature doesn't lie within coding region continue; } int peptidePosition = mapsTo[0]; String[][] codonVariants = variants.get(peptidePosition); if (codonVariants == null) { codonVariants = new String[3][]; variants.put(peptidePosition, codonVariants); } /* * extract dna variants to a string array */ String alls = (String) sf.getValue("alleles"); if (alls == null) { continue; } String[] alleles = alls.split(","); /* * get this peptides codon positions e.g. [3, 4, 5] or [4, 7, 10] */ int[] codon = peptidePosition == lastPeptidePostion ? lastCodon : MappingUtils.flattenRanges(dnaToProtein.locateInFrom( peptidePosition, peptidePosition)); lastPeptidePostion = peptidePosition; lastCodon = codon; /* * save nucleotide (and this variant) for each codon position */ for (int codonPos = 0; codonPos < 3; codonPos++) { String nucleotide = String.valueOf(dnaSeq .getCharAt(codon[codonPos] - 1)); if (codon[codonPos] == dnaCol) { /* * record current dna base and its alleles */ String[] dnaVariants = new String[alleles.length + 1]; dnaVariants[0] = nucleotide; System.arraycopy(alleles, 0, dnaVariants, 1, alleles.length); codonVariants[codonPos] = dnaVariants; } else if (codonVariants[codonPos] == null) { /* * record current dna base only * (at least until we find any variation and overwrite it) */ codonVariants[codonPos] = new String[] { nucleotide }; } } } } /* * scan codon variations, compute peptide variants and add to peptide sequence */ int count = 0; for (Entry variant : variants.entrySet()) { int peptidePos = variant.getKey(); String[][] codonVariants = variant.getValue(); String residue = String.valueOf(peptide.getCharAt(peptidePos - 1)); // 0-based List peptideVariants = computePeptideVariants(codonVariants, residue); if (!peptideVariants.isEmpty()) { Collections.sort(peptideVariants); String desc = StringUtils.listToDelimitedString(peptideVariants, ", "); SequenceFeature sf = new SequenceFeature( SequenceOntology.SEQUENCE_VARIANT, desc, peptidePos, peptidePos, Float.NaN, null); peptide.getDatasetSequence().addSequenceFeature(sf); count++; } } return count; } /** * Returns a non-redundant list of all peptide translations generated by the * given dna variants, excluding the current residue value * * @param codonVariants * an array of base values for codon positions 1, 2, 3 * @param residue * the current residue translation * @return */ protected static List computePeptideVariants( String[][] codonVariants, String residue) { List result = new ArrayList(); for (String base1 : codonVariants[0]) { for (String base2 : codonVariants[1]) { for (String base3 : codonVariants[2]) { String codon = base1 + base2 + base3; // TODO: report frameshift/insertion/deletion // and multiple-base variants?! String peptide = codon.contains("-") ? "-" : ResidueProperties .codonTranslate(codon); if (peptide != null && !result.contains(peptide) && !peptide.equals(residue)) { result.add(peptide); } } } } return result; } /** * Computes a list of all peptide variants given dna variants * * @param dnaSeq * the coding dna sequence * @param codonVariants * variant features for each codon position (null if no variant) * @param residue * the canonical protein translation * @return */ protected static List computePeptideVariants(SequenceI dnaSeq, SequenceFeature[] codonVariants, String residue) { List result = new ArrayList(); int[][] dnaVariants = new int[3][]; for (int i = 0; i < 3; i++) { } // TODO Auto-generated method stub return null; } /** * precalculate different products that can be found for seqs in dataset and * return them. * * @param dna * @param seqs * @param dataset * @param fake * - don't actually build lists - just get types * @return public static Object[] buildXProductsList(boolean dna, SequenceI[] * seqs, AlignmentI dataset, boolean fake) { String types[] = * jalview.analysis.CrossRef.findSequenceXrefTypes( dna, seqs, * dataset); if (types != null) { System.out.println("Xref Types for: * "+(dna ? "dna" : "prot")); for (int t = 0; t < types.length; t++) { * System.out.println("Type: " + types[t]); SequenceI[] prod = * jalview.analysis.CrossRef.findXrefSequences(seqs, dna, types[t]); * System.out.println("Found " + ((prod == null) ? "no" : "" + * prod.length) + " products"); if (prod!=null) { for (int p=0; * p