/*
* 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.util.DBRefUtils;
import jalview.util.MapList;
import jalview.ws.SequenceFetcher;
import jalview.ws.seqfetcher.ASequenceFetcher;
import java.util.ArrayList;
import java.util.List;
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
{
/*
* A sub-class that ignores Parent attribute when comparing sequence
* features. This avoids 'duplicate' CDS features that only
* differ in their parent Transcript ids.
*/
class MySequenceFeature extends SequenceFeature
{
private SequenceFeature feat;
MySequenceFeature(SequenceFeature sf)
{
this.feat = sf;
}
@Override
public boolean equals(Object o)
{
return feat.equals(o, true);
}
}
/**
* 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.getDBRefs());
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.getDBRefs());
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.getDBRefs());
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].getDBRefs());
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 seqs
* sequences whose xrefs are being retrieved
* @param dna
* true if sequences are nucleotide
* @param source
* @param al
* alignment to search for cross-referenced sequences (and possibly
* add to)
* @return products (as dataset sequences)
*/
public static Alignment findXrefSequences(SequenceI[] seqs,
final boolean dna, final String source, AlignmentI al)
{
AlignmentI dataset = al.getDataset() == null ? al : al.getDataset();
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.getDBRefs());
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.getDBRefs());
// 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());
}
}
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)
{
updateDbrefMappings(dna, seq, xrfs, retrieved, cf);
SequenceIdMatcher matcher = new SequenceIdMatcher(
dataset.getSequences());
List copiedFeatures = new ArrayList();
CrossRef me = new CrossRef();
for (int rs = 0; rs < retrieved.length; rs++)
{
// TODO: examine each sequence for 'redundancy'
DBRefEntry[] dbr = retrieved[rs].getDBRefs();
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)
{
SequenceI matched = matcher
.findIdMatch(map.getTo());
if (matched != null)
{
/*
* already got an xref to this sequence; update this
* map to point to the same sequence, and add
* any new dbrefs to it
*/
for (DBRefEntry ref : map.getTo().getDBRefs())
{
matched.addDBRef(ref); // add or update mapping
}
map.setTo(matched);
}
else
{
matcher.add(map.getTo());
}
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
&& ms.getSequenceAsString().equals(
dss.getSequenceAsString()))
// && mappedrg.getSequenceAsString().equals(
// loc.getSequenceAsString()))
{
String msg = "Mapping updated from "
+ ms.getName()
+ " to retrieved crossreference "
+ dss.getName();
System.out.println(msg);
// method to update all refs of existing To on
// retrieved sequence with dss and merge any props
// on To onto dss.
map.setTo(dss);
/*
* copy sequence features as well, avoiding
* duplication (e.g. same variation from 2
* transcripts)
*/
SequenceFeature[] sfs = ms
.getSequenceFeatures();
if (sfs != null)
{
for (SequenceFeature feat : sfs)
{
/*
* we override SequenceFeature.equals here (but
* not elsewhere) to ignore Parent attribute
* TODO not quite working yet!
*/
if (!copiedFeatures
.contains(me.new MySequenceFeature(
feat)))
{
dss.addSequenceFeature(feat);
copiedFeatures.add(feat);
}
}
}
cf.addMap(retrieved[rs].getDatasetSequence(),
dss, map.getMap());
}
else
{
cf.addMap(retrieved[rs].getDatasetSequence(),
map.getTo(), map.getMap());
}
} 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)
{
ral = new Alignment(rseqs.toArray(new SequenceI[rseqs.size()]));
if (cf != null && !cf.isEmpty())
{
ral.addCodonFrame(cf);
}
}
return ral;
}
/**
* Updates any empty mappings in the cross-references with one to a compatible
* retrieved sequence if found, and adds any new mappings to the
* AlignedCodonFrame
*
* @param dna
* @param mapFrom
* @param xrefs
* @param retrieved
* @param acf
*/
static void updateDbrefMappings(boolean dna, SequenceI mapFrom,
DBRefEntry[] xrefs, SequenceI[] retrieved, AlignedCodonFrame acf)
{
SequenceIdMatcher matcher = new SequenceIdMatcher(retrieved);
for (DBRefEntry xref : xrefs)
{
if (!xref.hasMap())
{
String targetSeqName = xref.getSource() + "|"
+ xref.getAccessionId();
SequenceI[] matches = matcher.findAllIdMatches(targetSeqName);
if (matches == null)
{
return;
}
for (SequenceI seq : matches)
{
MapList mapping = null;
if (dna)
{
mapping = AlignmentUtils.mapCdnaToProtein(seq, mapFrom);
}
else
{
mapping = AlignmentUtils.mapCdnaToProtein(mapFrom, seq);
if (mapping != null)
{
mapping = mapping.getInverse();
}
}
if (mapping != null)
{
xref.setMap(new Mapping(seq, mapping));
if (dna)
{
AlignmentUtils.computeProteinFeatures(mapFrom, seq, mapping);
}
if (dna)
{
acf.addMap(mapFrom, seq, mapping);
}
else
{
acf.addMap(seq, mapFrom, mapping.getInverse());
}
continue;
}
}
}
}
}
/**
* 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.getDBRefs(), 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;
}
/**
* 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