membrane abinitio on heteromultimeric complexes

[Anonymous]

Hi,

so far, I have successfully applied fold&dock to built symmetric homotetramers of the transmembrane helix of a single-anchored membrane protein. However, now I would like to model a hetero(di/tri)meric complex with each subunit being a single transmembrane a-helix. If at least one of the helices had a different membrane topology, I could simply fuse them into a single polypeptide chain, which Rosetta could fold. The problem is that all helices have the same membrane topology. I tried to introduce a transmembrane helix of opposite topology to connect the other two helices. This linker helix would of course interact with the other helices, too, so I wanted to keep this helix at a distance by using constraints. Unfortunately, distance constraints do not work with the membrane abinitio protocol. I modified the protocol according to https://www.rosettacommons.org/content/membrane-ab-initio-modeling-application-constraints, but it didn’t work for me.

Is there any way I can built these heteromultimers with Rosetta’s abinitio approach? Or can – alternatively – the same degree of flexibility of the helical structure in the abinitio protocol be achieved by a docking approach?

Thanks a lot for your help!

Oliver

[Anonymous]

Membrane code in Rosetta isn’t all that well maintained, so I’m no surprised that it’s rather brittle to starting conditions. I might recommend trying regular abinitio with your conditions, and see if that gives reasonable results. Membrane mode may not actually be adding all that much to the modeling.

Passed that, when you say “didn’t work”, what sort of error message are you getting? Depending on what the problem s, we might be able to knock together a workaround (no guarantees, though).

[Anonymous]

Looking at the file rosetta_source/src/protocols/topology_broker/AsymFoldandDockClaimer.cc, it looks like there is a bug in the setup that’s been fixed in the trunk version. You need to add the following to the AsymFoldandDockClaimer::initialize_dofs() function after the line “std::string chainID(“A_B”);”:


using namespace kinematics;
// If we don't have a docking jump we have to create it before calling setup_foldtree (strangely enough)
if ( pose.fold_tree().num_jump() == 0 ) {
FoldTree f(pose.fold_tree());
f.clear();
f.simple_tree( pose.total_residue() );
f.new_jump( 1, pose.total_residue() , chain_break_res_ );
pose.fold_tree( f );
}


Hopefully that should fix the issue.

[Anonymous]

Hi,
I have a similar issue – a trimer consisting of transmembrane monomers of >500 residues.
I want to remodel only a small part of one monomer but in the presence of the rest of the complex.
I try to mimic the above approach but I got stuck with the fragments file.
As I understood I should adjust the numbering to fit the second and third monomer and then concatenate those fragments files.
So in principle I should change in the fragments file the second field in the following line:

position: 1 neighbors: 200

and the number following the ‘P’ field (it’s a 15th column) in :

1dgw X 72 L L -112.472 129.829 179.474 -0.439 4.922 26.075 3 0.000 P 1 F 12

However, the format of the fragments file allows only a 3-digit number in the 15th column. So I cannot put there a number >999.

Does anybody has an idea how to do it?

The structure of the complex has been prepared by homology modeling and it’s quite certain so I wouldn’t like to change the positions of monomers like in the Fold-dock algorithm. I’d rather use the membrane ab-initio only for the 100-residues part which is not present in a template.

I use the following setup for the topology broker protocol:

USE_INPUT_POSE

CLAIMER MembraneTopologyClaimer
END_CLAIMER

CLAIMER RigidChunkClaimer
PDB starting_complex.pdbA
REGION_FILE complex_part.rigid
END_CLAIMER

The rigid part is 400 residues of the first monomers and the whole second and third monomer of the complex.

Thanks a lot for any suggestion!

[Anonymous]

Thanks,
actually I did something different. I joined monomers with 50-residues long, extended dummy chains and did not change the fragments files at all. As long as I model only a fragment in one monomer it seems that it works.
What bothers me slightly is the Rosetta energy. Do you know any way to discard those dummy chains from the energy calculations in the topology broker protocol? Or maybe it’s already discarded since it’s kept rigid?

[Anonymous]

Hi,
Unfortunately something is still wrong.
Even when I joined fragments files for 615-residue and 881-residue sequences which both corresponding to the trimer sequence the folding simulation failed. The right fragments are not being found.
The models are generated but the fragment which should be folded is just a straight line.

Below lines overflow the log file:

protocols.simple_moves.FragmentMover: couldn’t find fragment to insert!!
core.scoring.MembranePotential: MembraneCenter 561.580 354.646 24.300
core.scoring.MembranePotential: MembraneNormal -0.421 0.895 0.150
core.scoring.MembranePotential: ATOM 9999 X MEM A 999 561.580 354.646 24.300
core.scoring.MembranePotential: ATOM 9999 Y MEM A 999 555.271 368.068 26.545
core.scoring.MembranePotential: ATOM 9999 Z MEM A 999 567.889 341.223 22.055
protocols.simple_moves.FragmentMover: couldn’t find fragment to insert!!

I’ll gladly send you files via email, if needed.

Thanks,

D.

[Anonymous]

Just to add: the fragments file works if I use only a monomer structure. The missing fragment folds perfectly well then.

[Anonymous]

Also, I cut the sequence to 1000 residues, then to 800 residues and finally to 699 residues.
And only in case of 699 residues that overflowing line “protocols.simple_moves.FragmentMover: couldn’t find fragment to insert!!” disappeared.
So maybe there is some array size in the topology broker code that limits protein size?

[Anonymous]

From the example runs you emailed me, I tracked down several problems:

First off was the fragment file. Because you had concatenated two fragment files, you were missing 2 & 8 residues, respectively, from the 3 & 9 length fragment files where you spliced them together. (Because you don’t have partial-length fragments, you don’t get fragments for the n-1 residues in a protein.) It looks like Rosetta is getting confused about the gap in the fragment files, which causes issues. You need to paste in fragments for the bridging region. (Since it’s in the rigid region, the exact fragments you use probably won’t matter.)

Secondly, there’s an issue with your lipophilicity file. For some reason you’re showing values for residue 1497 — which is a problem, as the protein you’re using only has 1496 residues.

The big issue you’re seeing with the “couldn’t find fragment to insert!!” error has to do with a bug in the end bias fragment calculation.

If you run a fragment insertion with equal probability of trying an insertion at each position, what you see is that you don’t get out even insertion of fragments across the protein. Instead, due to lever arm effects and the like, you get bias in the frequency of acceptance at the beginning of the protein versus in the middle, versus at the end. For this reason, fragment insertion often has an “end bias” term turned on. This reweights the probability of inserting a fragment based on where the fragment is in the chain, and the reweighting takes the form of automatically rejecting a fragment insertion some fraction of the time (based on position of the insertion). This way you get a more even acceptance of fragment insertion across the protein.

The problem you ran into is that the weights computation was originally made based on the assumption that the full protein would allow fragment insertions. Your protein is mostly rigid, with comparatively few residues being allowed to have fragment insertions. There’s some attempt to compensate for rigid regions, but for your case the end bias check is thrown off completely, making it such that you’ll almost never get a fragment insertion. That’s why cutting down the length of the protein fixes the issue – the number of rigid residue versus the number of flexible residues drops down to a level where the end bias calculation starts spitting back somewhat reasonable numbers.

The way around this is to effectively turn off the end bias calculation. Unfortunately, I didn’t see a simple on/off flag for this. What you can do instead is put a really large value into the end bias constant, which will effectively turn off the end bias check. The default is 30, but something like “-abinitio:end_bias 10000” should address the issue – at least to the point where you should get fragment insertions.

[Anonymous]

Thanks a lot!
I’ll adjust the files and let know if I succeed.

[Anonymous]

It works fine!
Thanks a lot!
But there is just one more issue with the log file size.
Is there any way to silent lines like below to make the log files less heavy?

core.scoring.MembranePotential: MembraneCenter -10.892 27.640 75.111
core.scoring.MembranePotential: MembraneNormal -0.348 -0.934 -0.081
core.scoring.MembranePotential: ATOM 9999 X MEM A 999 -10.892 27.640 75.111
core.scoring.MembranePotential: ATOM 9999 Y MEM A 999 -16.111 13.629 73.900
core.scoring.MembranePotential: ATOM 9999 Z MEM A 999 -5.673 41.650 76.323

[Anonymous]

Great!
Thanks a lot for a detailed comment!

[Anonymous]

@membrane abinitio:
Is a regular abinitio approach really an alternative for a membrane protein? I find it hard to image that without an implicit membrane the transmembrane helices can be modelled with similar accuracy. Or do you mean instead of using the membrane abinitio protocol, that I should just use regular abinitio with a score weights file reflecting a membrane environment?

In the mean time I have thought about a strategy for assembling the complex. 1) model all helices separately by membrane abinitio; 2) choose top scorers&cluster them; 3) pairwise dockings of representatives of clusters of different complex subunits; 4) assembly complex models by selecting suitable top scoring dockings from 3.
Of course, step 3 is computationally very demanding (read “insane”) and step 4 won’t be fun either. But it all depends on steps 1&2 – what do you think, Rocco, would these first two steps enrich native-like conformations of the helices? Or is the structure of a single transmembrane helix diffusing freely in a membrane not necessarily very similar to the conformation in its bound state within a complex?

@membrane abinitio and distance constraints:
In the comment referred to by the link in my first posting, the user ytao reported that after altering the source code the constraints were listed in the score lists printed by Rosetta jobs for each stage. In my case Rosetta did not print any warning or error message. With the minirosetta application Rosetta was at least so kind to ackknowledge having read a constraint file, but still the score lists did not give any indication on whether the constraints were actually applied. The resulting models looked as if the constraints had not been applied. Anyway, this might be a dead end and not worth pursuing.

[Anonymous]

The membrane code and the membrane scorefunction is not as good as it conceivably could be – there just haven’t been enough work and benchmarking on it. Therefore, I’d anticipate that it’s possible that you might get equivalent results with and without the implicit membrane. I wouldn’t necessarily immediately launch a big job, but if you’re having issues with the membrane ab initio, I’d try a small-scale test of the regular ab initio, and see if the results you’re getting are reasonable. Trying it with a weights file that reduces the solvation penalty for exposed hydrophobics would be worth trying.

—

Unfortunately, from what I understand of what you’re doing, the best way to do it would be with a two-component symmetric system. While this is under current active development, the machinery to do it hasn’t been released yet, and I don’t believe it will be includes in the upcoming 3.5 release. If you want to use the symmetry machinery, you’ll need to model things as a single component symmetry.

Probably the best way to do that is to use the (undocumented) asymmetric fold and dock to join the two helicies, and then use symmetric docking (the non-folding kind) to find the symmetric assembly of the monomers.

To use the asymmetric fold&dock, you just change the broker setup file to something like:

CLAIMER AsymFoldandDockClaimer
loop_file flexible.loop
chain_break_assym_fnd 95
END_CLAIMER

Where 95 is the residue number of the end of the first chain. The loop_file is obligatory, and encompasses the entirety of the two proteins (e.g. all 95+53 residues, with an arbitrary cut point at residue 49)

LOOP 1 148 49 0 0

The other thing you need to do is to create a merged fasta of the two proteins (one sequence with both chains), and concatenate the fragment files, with the fragments for the second chain being shifted in position numbering to match that of the merged fasta files. Remove the flags for the symmetry from the command line, and otherwise it should be more-or-less like symmetric fold and dock.

[Anonymous]

Yes, asymmetric fold&dock is exactly what I was hoping for! Thanks a lot for the hint, Rocco! I followed your instructions and set up a test case. Unfortunately, Rosetta terminates with the error message:

ERROR: Unable to set up interface foldtree because there are no movable jumps
ERROR:: Exit from: src/protocols/docking/util.cc line: 289


Maybe I am just missing a setting or using a wrong one. I hope, you can help me with this.

broker setup file:

CLAIMER AsymFoldandDockClaimer
loop_file 2zw3_18-47_67-98.loop
chain_break_assym_fnd 30
END_CLAIMER

loop file:

LOOP 1 62 40 0 0


Here is my job script with all settings; in addition I have attached the log from the test run.

setenv rosetta3Dir /share/apps/rosetta3.4
setenv rosetta3 ${rosetta3Dir}/rosetta_source/bin/minirosetta.linuxgccrelease
setenv rosetta3DB ${rosetta3Dir}/rosetta_database
time ${rosetta3}
-database ${rosetta3DB}
-in:file:frag3 18-47_67-98.200.3mers
-in:file:frag9 18-47_67-98.200.9mers
-in:file:fasta 2zw3_18-47_67-98.fasta
-run:seed_offset $SGE_TASK_ID
-run:protocol broker
-broker:setup 2zw3_18-47_67-98.tpb
-abinitio:relax
-use_filters true
-psipred_ss2 2zw3_18-47_67-98.psipred_ss2
-abinitio:increase_cycles 1
-abinitio:rg_reweight 0.5
-abinitio:rsd_wt_helix 0.5
-abinitio:rsd_wt_loop 0.5
-user_tag j000$SGE_TASK_ID
-out:nstruct 1
-out:file:fullatom
-out:file:silent_struct_type binary
-out:file:silent default_$SGE_TASK_ID.out
-out:file:scorefile score_$SGE_TASK_ID.fsc


[Anonymous]

I might have discovered a problem. When I joined the fragment files for both chains I didn’t notice that for the last 2 (3mers) and 8 (9mers) residues of the first sequence no fragments are collected. Do I have to elongate the first sequence with the second sequence for fragment making? Is the error message (see above) a direct consequence of missing fragments? I thought the chain break introduced by the asym fold&dock protocol cuts the sequence and thus fragments for the last 2 and 8 residues of the first subsequence are not required. I will add the missing fragments and try running the asym fold&dock protocol again.

edit
I have added the missing fragments, but still the same error message is printed.
/edit

[Anonymous]

Yes, it’s working now!

Thank you very much for taking the time to fix this issue, Rocco!

[Anonymous]

Actually, as I read the code, the current version of the fragment file readers don’t actually look at anything above the 74th column. Some other utility may use the “3 0.000 P 1 F 12” data, but it doesn’t look like to actually use the fragments in Rosetta you need it.

[Anonymous]

If you’re not touching the dummy chains at all during the protocol (especially if you’re holding the things directly on either side rigid), you probably don’t even need the dummy linkers. Typically Rosetta doesn’t care if parts of a chain is missing. (Usually this is due to missing density in crystal structures and the like.) As long as you’re not modeling the rigid body/backbone degrees of freedom through the missing regions, you can just omit the residue. (Note that this will change the pose numbering of the residues downstream of the deleted region, so if fragments span the gap or are on the C-terminal side of the gap they’ll need to be adjusted.

Omitted regions won’t show up in the scorefunction, of course. (The geometry of the gap isn’t penalized unless you somehow set things up to consider the gap a chainbreak, which normally *isn’t* the default.)

[Anonymous]

Oh, certainly. Each of the tracers that print to the log can be controlled individually.

If you want to turn them off completely, you can use the -mute command line option, e.g. “-mute core.scoring.MembranePotential”. You can also selectively enable certain channels, by muting everything and unmuting the channels you want: “-mute all -unmute core.init”

That turns off those tracers completely, though. Perhaps slightly better is to adjust the “volume” of tracer by adjusting the output level. For example “-out:levels core.scoring.MembranePotential:warning” will turn off the regular “info” level messages from the core.scoring.MembranePotential channel, but will allow any “warning” or “error” messages to be printed. The full hierarchy is “error” “warning” “info” “debug” “trace”, where a particular level setting will let that level and higher of messages be printed.

For production runs, what I typically do is turn the default level up to warning (to make sure I catch any potential issues), and then selectively lower the levels on the tracers I’m interested in. (e.g. “-out:levels all:warning core.init:info protocols.jd2:info protocols.simple_filters.SidechainRmsdFilter:info” — which ones you’d want would be dependent on your run conditions). Note that the tracers should be hierarchical, so setting protocols.jd2 would also apply to protocols.jd2.JobDistributor, unless you also explicitly set things for protocols.jd2.JobDistributor too.

[Author]

Posts