OK - so the following will redirect the console output to a file, then look for
a matching chain from the sequence in test.fasta - after that I can grep
testing.log for the appropriate pattern and assign it to a best_chain_id
variable of some sort:
import subprocess
import sys
so = open("testing.log", 'w', 0)
sys.stdout = os.fdopen(sys.stdout.fileno(), 'w', 0)
os.dup2(so.fileno(), sys.stdout.fileno())
align_to_closest_chain("%s"%(subprocess.check_output(["tail","+2","test.fasta"])).upper(),
0.9)
On Dec 30, 2013, at 5:45 PM, Oliver Clarke <[email protected]> wrote:
> Great! Next scripting question that I’m puzzling over - how to take a set of
> sequences (present in a single fasta file), in arbitrary order, match and
> align them to their corresponding homologous subunits in a PDB file, and use
> align_and_mutate to correct the sequence of each subunit (and delete the
> gaps).
>
> Right now, I’m stuck on the first step - how to take a single sequence from a
> fasta file and output the matching chain ID.
>
> I can use:
>
> import subprocess
> align_to_closest_chain("%s"%(subprocess.check_output(["tail","+2","test.fasta"])).upper(),
> 0.95)
>
> to check for a matching chain, but the output of align_to_closest_chain() is
> either 0 or 1.
>
> The matching chain id is present in the console output, but I haven’t yet
> figured out how to parse that and do anything useful with it - the console
> output does not seem be easily accessible via either stdout or stderr (was
> hoping I could redirect the console output to a file and grep it for the
> chain id - ugly but I don’t know a better way).
>
> Also it is not immediately clear to me how the match_fraction parameter
> relates to sequence identity or alignment score for a given match (The
> alignment also seems to be case sensitive, which doesn’t seem right, hence
> the .upper() in the above).
>
> What am i doing wrong - thoughts/suggestions? Would it be better to do the
> sequence handling steps externally?
>
> Oliver.
>
>
>
>
> On Dec 30, 2013, at 2:43 PM, Paul Emsley <[email protected]> wrote:
>
>> On 30/12/13 17:42, Oliver Clarke wrote:
>>> And here is a more compact/elegant version courtesy of Paul that
>>> accomplishes a similar function in half the code - morphing with a radius
>>> decreasing from 12 to 6 Å over an arbitrary number of cycles, using a map
>>> initially blurred with a B-factor of 100, decreasing incrementally to 0
>>> (the original map) over the course of the procedure:
>>>
>>
>> I did some experimenting and I have been won over to the idea of an initial
>> RBR (into a blurred map), so adding that and fixing the typos, reso_morph
>> becomes:
>>
>>
>> def reso_morph(imol, imol_map, n_rounds):
>>
>> turn_off_backup(imol)
>> set_refinement_immediate_replacement(1)
>>
>> sharpen(imol_map, 150)
>> for ch_id in chain_ids(imol):
>> if is_polymer_chain(imol, ch_id):
>> rigid_body_refine_by_atom_selection(imol, '//'+ch_id)
>> accept_regularizement()
>>
>> for round in range(n_rounds):
>> f = float(round)/float(n_rounds)
>> for ch_id in chain_ids(imol):
>> if is_polymer_chain(imol, ch_id):
>>
>> # play with these numbers
>> radius = 8 * (2 - f) - 3
>> sf = 120 * (1 - f)
>>
>> sharpen(imol_map, sf)
>> morph_fit_chain(imol, ch_id, radius)
>>
>>
>
