001/** 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 * Created on 2013-05-28 021 * Created by Douglas Myers-Turnbull 022 * 023 * @since 3.0.6 024 */ 025package org.biojava.nbio.structure.io; 026 027import org.biojava.nbio.core.alignment.template.AlignedSequence; 028import org.biojava.nbio.core.alignment.template.SequencePair; 029import org.biojava.nbio.structure.*; 030import org.biojava.nbio.structure.align.model.AFPChain; 031import org.biojava.nbio.structure.align.util.AlignmentTools; 032import org.biojava.nbio.structure.align.xml.AFPChainXMLConverter; 033import org.biojava.nbio.core.exceptions.CompoundNotFoundException; 034import org.biojava.nbio.core.sequence.ProteinSequence; 035import org.biojava.nbio.core.sequence.compound.AminoAcidCompound; 036import org.biojava.nbio.core.sequence.compound.AminoAcidCompoundSet; 037import org.biojava.nbio.core.sequence.io.CasePreservingProteinSequenceCreator; 038import org.biojava.nbio.core.sequence.io.FastaReader; 039import org.biojava.nbio.core.sequence.io.GenericFastaHeaderParser; 040import org.biojava.nbio.core.sequence.io.template.SequenceCreatorInterface; 041import org.biojava.nbio.core.sequence.io.template.SequenceHeaderParserInterface; 042import org.biojava.nbio.core.sequence.template.Sequence; 043import org.biojava.nbio.core.util.SequenceTools; 044import org.slf4j.Logger; 045import org.slf4j.LoggerFactory; 046 047import java.io.File; 048import java.io.FileInputStream; 049import java.io.IOException; 050import java.io.InputStream; 051import java.util.*; 052 053/** 054 * A collection of static utilities to convert between {@link AFPChain AFPChains} and {@link FastaSequence FastaSequences}. 055 * 056 * @author dmyersturnbull 057 * @see StructureSequenceMatcher 058 * @see FastaStructureParser 059 * @see SeqRes2AtomAligner 060 */ 061public class FastaAFPChainConverter { 062 063 private final static Logger logger = LoggerFactory.getLogger(FastaAFPChainConverter.class); 064 065 066 public static AFPChain cpFastaToAfpChain(String first, String second, Structure structure, int cpSite) throws StructureException, CompoundNotFoundException { 067 ProteinSequence s1 = new ProteinSequence(first); 068 s1.setUserCollection(getAlignedUserCollection(first)); 069 ProteinSequence s2 = new ProteinSequence(second); 070 s2.setUserCollection(getAlignedUserCollection(second)); 071 return cpFastaToAfpChain(s1, s2, structure, cpSite); 072 } 073 074 /** 075 * Takes a structure and sequence corresponding to an alignment between a structure or sequence and itself (or even a structure with a sequence), where the result has a circular permutation site 076 * {@link cpSite} residues to the right. 077 * 078 * @param fastaFile A FASTA file containing exactly 2 sequences, the first unpermuted and the second permuted 079 * @param cpSite 080 * The number of residues from the beginning of the sequence at which the circular permutation site occurs; can be positive or negative; values greater than the length of the sequence 081 * are acceptable 082 * @throws IOException 083 * @throws StructureException 084 */ 085 public static AFPChain cpFastaToAfpChain(File fastaFile, Structure structure, int cpSite) throws IOException, StructureException { 086 InputStream inStream = new FileInputStream(fastaFile); 087 SequenceCreatorInterface<AminoAcidCompound> creator = new CasePreservingProteinSequenceCreator(AminoAcidCompoundSet.getAminoAcidCompoundSet()); 088 SequenceHeaderParserInterface<ProteinSequence, AminoAcidCompound> headerParser = new GenericFastaHeaderParser<ProteinSequence, AminoAcidCompound>(); 089 FastaReader<ProteinSequence, AminoAcidCompound> fastaReader = new FastaReader<ProteinSequence, AminoAcidCompound>(inStream, headerParser, creator); 090 LinkedHashMap<String, ProteinSequence> sequences = fastaReader.process(); 091 inStream.close(); 092 Iterator<ProteinSequence> iter = sequences.values().iterator(); 093 ProteinSequence first = iter.next(); 094 ProteinSequence second = iter.next(); 095 return cpFastaToAfpChain(first, second, structure, cpSite); 096 } 097 098 /** 099 * Takes a structure and sequence corresponding to an alignment between a structure or sequence and itself (or even a structure with a sequence), where the result has a circular permutation site 100 * {@link cpSite} residues to the right. 101 * 102 * @param first The unpermuted sequence 103 * @param second The sequence permuted by cpSite 104 * @param cpSite 105 * The number of residues from the beginning of the sequence at which the circular permutation site occurs; can be positive or negative; values greater than the length of the sequence 106 * are acceptable 107 * @throws StructureException 108 */ 109 public static AFPChain cpFastaToAfpChain(ProteinSequence first, ProteinSequence second, Structure structure, int cpSite) 110 throws StructureException { 111 112 if (structure == null) { 113 throw new IllegalArgumentException("The structure is null"); 114 } 115 116 if (first == null) { 117 throw new IllegalArgumentException("The sequence is null"); 118 } 119 120 // we need to find the ungapped CP site 121 int gappedCpShift = 0; 122 int ungappedCpShift = 0; 123 while (ungappedCpShift < Math.abs(cpSite)) { 124 char c; 125 try { 126 if (cpSite <= 0) { 127 c = second.getSequenceAsString().charAt(gappedCpShift); 128 } else { 129 c = second.getSequenceAsString().charAt(first.getLength()-1 - gappedCpShift); 130 } 131 } catch (StringIndexOutOfBoundsException e) { 132 throw new IllegalArgumentException("CP site of " + cpSite + " is wrong"); 133 } 134 if (c != '-') { 135 ungappedCpShift++; 136 } 137 gappedCpShift++; 138 } 139 140 Atom[] ca1 = StructureTools.getRepresentativeAtomArray(structure); 141 Atom[] ca2 = StructureTools.getRepresentativeAtomArray(structure); // can't use cloneCAArray because it doesn't set parent group.chain.structure 142 143 ProteinSequence antipermuted = null; 144 try { 145 antipermuted = new ProteinSequence(SequenceTools.permuteCyclic(second.getSequenceAsString(), gappedCpShift)); 146 } catch (CompoundNotFoundException e) { 147 // this can't happen, the original sequence comes from a ProteinSequence 148 logger.error("Unexpected error while creating protein sequence: {}. This is most likely a bug.",e.getMessage() ); 149 } 150 151 ResidueNumber[] residues = StructureSequenceMatcher.matchSequenceToStructure(first, structure); 152 ResidueNumber[] antipermutedResidues = StructureSequenceMatcher.matchSequenceToStructure(antipermuted, structure); 153 154 ResidueNumber[] nonpermutedResidues = new ResidueNumber[antipermutedResidues.length]; 155 SequenceTools.permuteCyclic(antipermutedResidues, nonpermutedResidues, -gappedCpShift); 156 157 // nullify ResidueNumbers that have a lowercase sequence character 158 if (first.getUserCollection() != null) { 159 CasePreservingProteinSequenceCreator.setLowercaseToNull(first, residues); 160 } 161 if (second.getUserCollection() != null) { 162 CasePreservingProteinSequenceCreator.setLowercaseToNull(second, nonpermutedResidues); 163 } 164 165// for (int i = 0; i < residues.length; i++) { 166// if (residues[i] == null) { 167// System.out.print("="); 168// } else { 169// System.out.print(sequence.getSequenceAsString().charAt(i)); 170// } 171// } 172// System.out.println(); 173// for (int i = 0; i < residues.length; i++) { 174// if (nonpermutedResidues[i] == null) { 175// System.out.print("="); 176// } else { 177// System.out.print(second.getSequenceAsString().charAt(i)); 178// } 179// } 180// System.out.println(); 181 182 return buildAlignment(ca1, ca2, residues, nonpermutedResidues); 183 184 } 185 186 /** 187 * Reads the file {@code fastaFile}, expecting exactly two sequences which give a pairwise alignment. Uses this and two structures to create an AFPChain corresponding to the alignment. Uses a 188 * {@link CasePreservingProteinSequenceCreator} and assumes that a residue is aligned if and only if it is given by an uppercase letter. 189 * 190 * @see #fastaToAfpChain(ProteinSequence, ProteinSequence, Structure, Structure) 191 * @throws IOException 192 * @throws StructureException 193 */ 194 public static AFPChain fastaFileToAfpChain(File fastaFile, Structure structure1, Structure structure2) 195 throws IOException, StructureException { 196 InputStream inStream = new FileInputStream(fastaFile); 197 SequenceCreatorInterface<AminoAcidCompound> creator = new CasePreservingProteinSequenceCreator( 198 AminoAcidCompoundSet.getAminoAcidCompoundSet()); 199 SequenceHeaderParserInterface<ProteinSequence, AminoAcidCompound> headerParser = new GenericFastaHeaderParser<ProteinSequence, AminoAcidCompound>(); 200 FastaReader<ProteinSequence, AminoAcidCompound> fastaReader = new FastaReader<ProteinSequence, AminoAcidCompound>( 201 inStream, headerParser, creator); 202 LinkedHashMap<String, ProteinSequence> sequences = fastaReader.process(); 203 inStream.close(); 204 return fastaToAfpChain(sequences, structure1, structure2); 205 } 206 207 /** 208 * Returns an AFPChain corresponding to the alignment between {@code structure1} and {@code structure2}, which is given by the gapped protein sequences {@code sequence1} and {@code sequence2}. The 209 * sequences need not correspond to the entire structures, since local alignment is performed to match the sequences to structures. 210 * @throws StructureException 211 * @throws CompoundNotFoundException 212 */ 213 public static AFPChain fastaStringToAfpChain(String sequence1, String sequence2, Structure structure1, 214 Structure structure2) throws StructureException, CompoundNotFoundException { 215 ProteinSequence seq1 = new ProteinSequence(sequence1); 216 ProteinSequence seq2 = new ProteinSequence(sequence2); 217 return fastaToAfpChain(seq1, seq2, structure1, structure2); 218 } 219 220 /** 221 * Uses two sequences each with a corresponding structure to create an AFPChain corresponding to the alignment. Provided only for convenience since FastaReaders return such maps. 222 * 223 * @param sequences 224 * A Map containing exactly two entries from sequence names as Strings to gapped ProteinSequences; the name is ignored 225 * @see #fastaToAfpChain(ProteinSequence, ProteinSequence, Structure, Structure) 226 * @throws StructureException 227 */ 228 public static AFPChain fastaToAfpChain(Map<String, ProteinSequence> sequences, Structure structure1, 229 Structure structure2) throws StructureException { 230 231 if (sequences.size() != 2) { 232 throw new IllegalArgumentException("There must be exactly 2 sequences, but there were " + sequences.size()); 233 } 234 235 if (structure1 == null || structure2 == null) { 236 throw new IllegalArgumentException("A structure is null"); 237 } 238 239 List<ProteinSequence> seqs = new ArrayList<ProteinSequence>(); 240 List<String> names = new ArrayList<String>(2); 241 for (Map.Entry<String, ProteinSequence> entry : sequences.entrySet()) { 242 seqs.add(entry.getValue()); 243 names.add(entry.getKey()); 244 } 245 246 return fastaToAfpChain(seqs.get(0), seqs.get(1), structure1, structure2); 247 } 248 249 /** 250 * TODO Write comment 251 * @param sequence1 252 * @param sequence2 253 * @param structure1 254 * @param structure2 255 * @return 256 * @throws StructureException 257 * @throws CompoundNotFoundException 258 */ 259 public static AFPChain fastaToAfpChain(String sequence1, String sequence2, Structure structure1, 260 Structure structure2) throws StructureException, CompoundNotFoundException { 261 ProteinSequence s1 = new ProteinSequence(sequence1); 262 s1.setUserCollection(getAlignedUserCollection(sequence1)); 263 ProteinSequence s2 = new ProteinSequence(sequence2); 264 s2.setUserCollection(getAlignedUserCollection(sequence2)); 265 return fastaToAfpChain(s1, s2, structure1, structure2); 266 } 267 268 /** 269 * Returns an AFPChain corresponding to the alignment between {@code structure1} and {@code structure2}, which is given by the gapped protein sequences {@code sequence1} and {@code sequence2}. The 270 * sequences need not correspond to the entire structures, since local alignment is performed to match the sequences to structures. Assumes that a residue is aligned if and only if it is given by 271 * an uppercase letter. 272 * @param sequence1 <em>Must</em> have {@link ProteinSequence#getUserCollection()} set to document upper- and lower-case as aligned and unaligned; see {@link #getAlignedUserCollection(String)} 273 * @throws StructureException 274 */ 275 public static AFPChain fastaToAfpChain(ProteinSequence sequence1, ProteinSequence sequence2, Structure structure1, 276 Structure structure2) throws StructureException { 277 278 if (structure1 == null || structure2 == null) { 279 throw new IllegalArgumentException("A structure is null"); 280 } 281 282 if (sequence1 == null || sequence2 == null) { 283 throw new IllegalArgumentException("A sequence is null"); 284 } 285 286 Atom[] ca1 = StructureTools.getRepresentativeAtomArray(structure1); 287 Atom[] ca2 = StructureTools.getRepresentativeAtomArray(structure2); 288 289 ResidueNumber[] residues1 = StructureSequenceMatcher.matchSequenceToStructure(sequence1, structure1); 290 ResidueNumber[] residues2 = StructureSequenceMatcher.matchSequenceToStructure(sequence2, structure2); 291 292 // nullify ResidueNumbers that have a lowercase sequence character 293 if (sequence1.getUserCollection() != null) { 294 CasePreservingProteinSequenceCreator.setLowercaseToNull(sequence1, residues1); 295 } 296 if (sequence2.getUserCollection() != null) { 297 CasePreservingProteinSequenceCreator.setLowercaseToNull(sequence2, residues2); 298 } 299 300 return buildAlignment(ca1, ca2, residues1, residues2); 301 302 } 303 304 /** 305 * Provided only for convenience. 306 * 307 * @see #fastaToAfpChain(ProteinSequence, ProteinSequence, Structure, Structure) 308 * @throws StructureException 309 */ 310 public static AFPChain fastaToAfpChain(SequencePair<Sequence<AminoAcidCompound>, AminoAcidCompound> alignment, 311 Structure structure1, Structure structure2) throws StructureException { 312 List<AlignedSequence<Sequence<AminoAcidCompound>, AminoAcidCompound>> seqs = alignment.getAlignedSequences(); 313 StringBuilder sb1 = new StringBuilder(); 314 for (AminoAcidCompound a : seqs.get(0)) { 315 sb1.append(a.getBase()); 316 } 317 try { 318 ProteinSequence seq1 = new ProteinSequence(sb1.toString()); 319 StringBuilder sb2 = new StringBuilder(); 320 for (AminoAcidCompound a : seqs.get(1)) { 321 sb1.append(a.getBase()); 322 } 323 ProteinSequence seq2 = new ProteinSequence(sb2.toString()); 324 LinkedHashMap<String, ProteinSequence> map = new LinkedHashMap<String, ProteinSequence>(); 325 map.put(structure1.getName(), seq1); 326 map.put(structure2.getName(), seq2); 327 return fastaToAfpChain(map, structure1, structure2); 328 } catch (CompoundNotFoundException e) { 329 logger.error("Unexpected error while creating protein sequences: {}. This is most likely a bug.",e.getMessage()); 330 return null; 331 } 332 } 333 334 /** 335 * Builds an {@link AFPChain} from already-matched arrays of atoms and residues. 336 * 337 * @param ca1 338 * An array of atoms in the first structure 339 * @param ca2 340 * An array of atoms in the second structure 341 * @param residues1 342 * An array of {@link ResidueNumber ResidueNumbers} in the first structure that are aligned. Only null ResidueNumbers are considered to be unaligned 343 * @param residues2 344 * An array of {@link ResidueNumber ResidueNumbers} in the second structure that are aligned. Only null ResidueNumbers are considered to be unaligned 345 * @throws StructureException 346 */ 347 private static AFPChain buildAlignment(Atom[] ca1, Atom[] ca2, ResidueNumber[] residues1, ResidueNumber[] residues2) 348 throws StructureException { 349 350 // remove any gap 351 // this includes the ones introduced by the nullifying above 352 List<ResidueNumber> alignedResiduesList1 = new ArrayList<ResidueNumber>(); 353 List<ResidueNumber> alignedResiduesList2 = new ArrayList<ResidueNumber>(); 354 for (int i = 0; i < residues1.length; i++) { 355 if (residues1[i] != null && residues2[i] != null) { 356 alignedResiduesList1.add(residues1[i]); 357 alignedResiduesList2.add(residues2[i]); 358 } 359 } 360 361 ResidueNumber[] alignedResidues1 = alignedResiduesList1.toArray(new ResidueNumber[alignedResiduesList1.size()]); 362 ResidueNumber[] alignedResidues2 = alignedResiduesList2.toArray(new ResidueNumber[alignedResiduesList2.size()]); 363 364 AFPChain afpChain = AlignmentTools.createAFPChain(ca1, ca2, alignedResidues1, alignedResidues2); 365 afpChain.setAlgorithmName("unknown"); 366 367 AlignmentTools.updateSuperposition(afpChain, ca1, ca2); 368 369 afpChain.setBlockSize(new int[] {afpChain.getNrEQR()}); 370 afpChain.setBlockRmsd(new double[] {afpChain.getTotalRmsdOpt()}); 371 afpChain.setBlockGap(new int[] {afpChain.getGapLen()}); 372 373 return afpChain; 374 375 } 376 377 /** 378 * Takes a protein sequence string with capital and lowercase letters and sets its {@link ProteinSequence#getUserCollection() user collection} to record which letters are uppercase (aligned) and which are lowercase (unaligned). 379 * @param sequence Make sure <em>not</em> to use {@link ProteinSequence#getSequenceAsString()} for this, as it won't preserve upper- and lower-case 380 */ 381 public static List<Object> getAlignedUserCollection(String sequence) { 382 List<Object> aligned = new ArrayList<Object>(sequence.length()); 383 for (char c : sequence.toCharArray()) { 384 aligned.add(Character.isUpperCase(c)); 385 } 386 return aligned; 387 } 388 389 /** 390 * Prints out the XML representation of an AFPChain from a file containing exactly two FASTA sequences. 391 * 392 * @param args 393 * A String array of fasta-file structure-1-name structure-2-name 394 * @throws StructureException 395 * @throws IOException 396 */ 397 public static void main(String[] args) throws StructureException, IOException { 398 if (args.length != 3) { 399 System.err.println("Usage: FastaAFPChainConverter fasta-file structure-1-name structure-2-name"); 400 return; 401 } 402 File fasta = new File(args[0]); 403 Structure structure1 = StructureTools.getStructure(args[1]); 404 Structure structure2 = StructureTools.getStructure(args[2]); 405 if (structure1 == null) throw new IllegalArgumentException("No structure for " + args[1] + " was found"); 406 if (structure2 == null) throw new IllegalArgumentException("No structure for " + args[2] + " was found"); 407 AFPChain afpChain = fastaFileToAfpChain(fasta, structure1, structure2); 408 String xml = AFPChainXMLConverter.toXML(afpChain); 409 System.out.println(xml); 410 } 411 412}