基础知识
数据处理案例

Relion_subtomo_averaging_code_review

Cryo-ET subtomo averaging是冷冻电镜解出高分辨率生物结构的重要过程，本文是对relion subtomo部分主要程序功能的代码阅读，对relion的解读以及部分数据示例主要来自源码阅读以及relion教程。

relion subtomo的基本流程

  1. 数据输入：
    • 现有软件（如IMOD、Warp等）做好的ET bin1数据（tilt series）与align参数
    • 现有软件（如CTFFIND等）估计的ctf参数
    • 现有软件（如dynamo等）解出的particles 3D位置与对应欧拉角
  2. Make pseudo-subtomograms：根据每个particle位置使用FBP重建出每个particle的3D volume。第一步先以bin4分辨率重建。
  3. De novo 3D model generation：随机挑选bin4的particle volume average出一个初始的initial model。
  4. Initial 3D refinement：以上一步的model为reference，计算FSC并使用SGD优化每一个particle的坐标与欧拉角直到收敛。
  5. Reconstruct particle：使用优化好的particle的所有2D crop image，使用FBP重建出bin2分辨率的reference map。与第2步不同的是，这一步不是每个particle单独重建自己的volume，而是所有particle一起重建一个reference map。
  6. 重复第2步，用第4步的优化结果以bin2分辨率重建每个particle。
  7. 3D auto refine：以第五步结果为reference map，优化bin2 particle。
  8. 重复5、6、7步，在bin1分辨率上重建与优化，得到最终高分辨率reference map。

接下来会详细介绍上述的具体步骤，每个步骤可大致分为数据输入、数据输出以及代码解析三部分，代码解析中会有用中文写的注释进行说明。
注：由于3D initial model与refine部分的内容非常重要，需要进一步的阅读与考证，因此本文的解析只限于前两个步骤，之后会持续更新。

Main Step 1：数据输入（import tomograms and import coordinates）

在这一步，relion会把ET tomogram与对应particle信息全部输入进来，对于tomogram，需要原始tilt series、align参数文件（imod）、ctf参数（ctffind或ctfplotter）；对于particles，需要每个particle的3D坐标与欧拉旋转角。
源码位置：source/jaz/tomography/apps/import_tomograms.cpp source/jaz/tomography/apps/import_particles.cpp

  • Import Tomograms：输入tilt series与对应信息，.star文件如下。 
      data_

  loop_
  _rlnTomoName
  _rlnTomoTiltSeriesName
  _rlnTomoImportCtfFindFile
  _rlnTomoImportImodDir
  _rlnTomoImportFractionalDose
  _rlnTomoImportOrderList

  TS_01    tomograms/TS_01/01.mrc  tomograms/TS_01/01_output.txt    tomograms/TS_01    3.0    input/order_list.csv 
  TS_03    tomograms/TS_03/03.mrc  tomograms/TS_03/03_output.txt    tomograms/TS_03    3.0    input/order_list.csv 
  TS_43    tomograms/TS_43/43.mrc  tomograms/TS_43/43_output.txt    tomograms/TS_43    3.1    input/order_list.csv 
  TS_45    tomograms/TS_45/45.mrc  tomograms/TS_45/45_output.txt    tomograms/TS_45    3.1    input/order_list.csv 
  TS_54    tomograms/TS_54/54.mrc  tomograms/TS_54/54_output.txt    tomograms/TS_54    3.0    input/order_list.csv

    几乎所有relion的输入都可以被描述成.star文件。

    • rlnTomoName: tomogram的名字
    • rlnTomoTiltSeriesName: tomogram raw stack mrc路径
    • rlnTomoImportCtfFindFile: 对应的ctf参数文件（如果是ctfplotter的结果就是rlnTomoImportCtfPlotterFile）
    • rlnTomoImportImodDir: tomogram对应IMOD文件夹路径（用来获得align参数）
    • rlnTomoImportFractionalDose: 数据拍摄用的电子剂量
    • rlnTomoImportOrderList: 标注数据拍摄顺序，表格里每行有两个数（i，j），意味着第i个被拍的角度是j。
  • Import coordinates：输入particle坐标、欧拉角与对应的tomograms名字。 
      data_particles

  loop_
  _rlnTomoName #1
  _rlnTomoParticleId #2
  _rlnTomoManifoldIndex #3
  _rlnCoordinateX #4
  _rlnCoordinateY #5
  _rlnCoordinateZ #6
  _rlnOriginXAngst #7
  _rlnOriginYAngst #8
  _rlnOriginZAngst #9
  _rlnAngleRot #10
  _rlnAngleTilt #11
  _rlnAnglePsi #12
  _rlnClassNumber #13
  _rlnRandomSubset #14
  TS_01	1	1	644	1164	72	-0.000	-0.000	-0.000	-8.72	167.29	87.57	-1	2
  TS_01	2	1	624	1112	64	-0.000	-0.000	-0.000	-130.65	175.35	84.29	-1	1
  TS_01	3	1	676	1120	68	-0.000	-0.000	-0.000	-117.88	170.53	95.89	-1	2
  TS_01	4	1	696	1168	84	-0.000	-0.000	-0.000	53.64	158.62	86.53	-1	1
  TS_01	5	1	592	1156	72	-0.000	-0.000	-0.000	-78.74	169.66	78.39	-1	2
  TS_01	6	1	659	1071	67	0.026	-0.088	0.012	54.34	178.38	-31.41	1	1
  TS_01	7	1	712	1076	64	-0.000	-0.000	-0.000	61.58	177.85	156.65	-1	2
  TS_01	8	1	660	1200	100	-0.000	-0.000	-0.000	106.21	148.26	74.99	-1	1
  ...

    数据里的3D坐标应该与IMOD文件夹内tilt.com文件对应。

      $tilt -StandardInput
  InputProjections 01_dose-filt.ali
  OutputFile 01_dose-filt_full.rec
  IMAGEBINNED 1
  TILTFILE 01_dose-filt.tlt
  THICKNESS 1360
  RADIAL 0.35 0.05
  XAXISTILT 0.0
  LOG 0.0
  SCALE 0.0 330.0
  PERPENDICULAR 
  MODE 2
  FULLIMAGE 3710 3710
  SUBSETSTART 0 1825
  AdjustOrigin 
  ActionIfGPUFails 1,2
  XTILTFILE 01_dose-filt.xtilt
  OFFSET 0.0
  SHIFT 0.0 300.0
  UseGPU 0
  $if (-e ./savework) ./savework

    THICHNESS与FULLIMAGE就是3D volume的维度。

Main step 2：Make pseudo-subtomograms

数据输入后，relion通过particle坐标投影回2D图片进行FBP重建。
源码位置: source/jaz/tomography/programs/subtomo.cpp

  • 数据输入：
    • tomograms.star:
      data_global

  loop_ 
  _rlnTomoName #1 
  _rlnTomoTiltSeriesName #2 
  _rlnTomoFrameCount #3 
  _rlnTomoSizeX #4 
  _rlnTomoSizeY #5 
  _rlnTomoSizeZ #6 
  _rlnTomoHand #7 
  _rlnOpticsGroupName #8 
  _rlnTomoTiltSeriesPixelSize #9 
  _rlnVoltage #10 
  _rlnSphericalAberration #11 
  _rlnAmplitudeContrast #12 
  _rlnTomoImportFractionalDose #13 
      TS_01 tomograms/TS_01/01.mrc           40         3710         3710         1360     -1.00000 opticsGroup1     1.350000   300.000000     2.700000     0.070000     3.000000 
      TS_03 tomograms/TS_03/03.mrc           40         3710         3710         1120     -1.00000 opticsGroup1     1.350000   300.000000     2.700000     0.070000     3.000000 
      TS_43 tomograms/TS_43/43.mrc           41         3710         3710         1200     -1.00000 opticsGroup1     1.350000   300.000000     2.700000     0.070000     3.100000 
      TS_45 tomograms/TS_45/45.mrc           41         3710         3710         1200     -1.00000 opticsGroup1     1.350000   300.000000     2.700000     0.070000     3.100000 
      TS_54 tomograms/TS_54/54.mrc           41         3710         3710         1200     -1.00000 opticsGroup1     1.350000   300.000000     2.700000     0.070000     3.000000

  # version 30001

  data_TS_01

  loop_ 
  _rlnTomoProjX #1 
  _rlnTomoProjY #2 
  _rlnTomoProjZ #3 
  _rlnTomoProjW #4 
  _rlnDefocusU #5 
  _rlnDefocusV #6 
  _rlnDefocusAngle #7 
  _rlnCtfScalefactor #8 
  _rlnMicrographPreExposure #9 
  [0.0446311671183,-0.996636282201,-0.0687259705396,3671.47512481] [0.542807022999,0.0819463245359,-0.835849516828,1134.92684561] [0.838670567945,0,0.544639035015,-1762.57143309]  [0,0,0,1] 36642.020000 36513.230000    -21.44723     1.000000   114.000000 
  [0.0481668410418,-0.996636282201,-0.0662959691761,3652.3982384] [0.585808108577,0.0819463245359,-0.806295689511,1032.33863359] [0.809016994375,0,0.587785252292,-1723.90164368]  [0,0,0,1] 40659.750000 40045.230000     -0.91484     1.000000   105.000000 
  [0.051570493002,-0.996636282201,-0.0636842551693,3652.69289457] [0.627203534849,0.0819463245359,-0.774531861754,962.935527375] [0.777145961457,0,0.62932039105,-1680.50675969]  [0,0,0,1] 40413.310000 40192.590000     -9.89749     1.000000   102.000000 
  [0.0548327828959,-0.996636183544,-0.0608979749071,3639.0873037] [0.666879773814,0.0819463082018,-0.740645022684,862.386149688] [0.743144825477,0,0.669130606359,-1632.50572348]  [0,0,0,1] 39817.330000 38824.790000    22.175120     1.000000    93.000000 
  [0.0579448017726,-0.996636282201,-0.0579448017726,3637.20509583] [0.704728273521,0.0819463245359,-0.704728273521,794.060306913] [0.707106781187,0,0.707106781187,-1580.03010256]  [0,0,0,1] 39895.570000 39376.050000    -59.23525     1.000000    90.000000 
  [0.0608979870457,-0.996636282201,-0.0548327938256,3623.02400791] [0.740645096001,0.0819463245359,-0.666879839828,697.175028762] [0.669130606359,0,0.743144825477,-1523.22372895]  [0,0,0,1] 40997.750000 40550.260000    41.593690     1.000000    81.000000 
  [0.0636842551693,-0.996636282201,-0.051570493002,3620.16706308] [0.774531861754,0.0819463245359,-0.627203534849,636.106604797] [0.62932039105,0,0.777145961457,-1462.24230499]  [0,0,0,1] 40370.910000 39761.240000     8.227517     1.000000    78.000000 
  [0.0662959691761,-0.996636282201,-0.0481668410418,3607.70796738] [0.806295689511,0.0819463245359,-0.585808108577,543.250808196] [0.587785252292,0,0.809016994375,-1397.2529765]  [0,0,0,1] 39999.480000 39411.810000    70.785380     1.000000    69.000000 
  [0.0687259705396,-0.996636282201,-0.0446311671183,3605.5238664] [0.835849516828,0.0819463245359,-0.542807022999,488.766946972] [0.544639035015,0,0.838670567945,-1328.43387472]  [0,0,0,1] 40137.500000 39760.600000    -58.24971     1.000000    66.000000 
  [0.0709675987948,-0.996636282201,-0.0409731622679,3592.33824821] [0.863112338719,0.0819463245359,-0.4983181411,401.401576476] [0.5,0,0.866025403784,-1255.97362803]  [0,0,0,1] 39349.730000 38793.070000    88.441160     1.000000    57.000000 
  [0.0730147097947,-0.996636282201,-0.0372028528279,3588.68703015] [0.888009429684,0.0819463245359,-0.452463403815,355.726113698] [0.45399049974,0,0.891006524188,-1180.07084493]  [0,0,0,1] 39822.450000 39138.010000    -88.34856     1.000000    54.000000 
  [0.0748616925503,-0.996636282201,-0.0333305729541,3578.55399784] [0.910472548526,0.0819463245359,-0.40536849579,274.52722088] [0.406736643076,0,0.913545457643,-1100.93356967]  [0,0,0,1] 39827.790000 38858.570000    37.508620     1.000000    45.000000 
  [0.0765034846101,-0.996636282201,-0.0293669362967,3580.20836769] [0.9304401254,0.0819463245359,-0.357162500895,241.405650451] [0.358367949545,0,0.933580426497,-1018.77871202]  [0,0,0,1] 39505.590000 38739.020000    72.908280     1.000000    42.000000 
  [0.0779355859363,-0.996636282201,-0.0253228069082,3571.11169502] [0.947857430563,0.0819463245359,-0.307977548411,168.405709197] [0.309016994375,0,0.951056516295,-933.831452739]  [0,0,0,1] 39406.160000 38893.030000    52.051360     1.000000    33.000000 
  [0.0791540712387,-0.996636282201,-0.021209269466,3570.27766288] [0.962676724395,0.0819463245359,-0.257948450874,128.42558249] [0.258819045103,0,0.965925826289,-846.324626351]  [0,0,0,1] 39531.560000 39091.820000    58.772460     1.000000    30.000000 
  [0.0801556007337,-0.996636282201,-0.0170375988906,3564.42952706] [0.974857388239,0.0819463245359,-0.207212334563,83.6588849601] [0.207911690818,0,0.978147600734,-756.49808299]  [0,0,0,1] 39427.630000 39046.290000    59.961690     1.000000    21.000000 
  [0.0809374292987,-0.996636282201,-0.0128192294408,3562.97333088] [0.984366035744,0.0819463245359,-0.155908263646,29.0829757451] [0.15643446504,0,0.987688340595,-664.598030968]  [0,0,0,1] 39532.620000 38599.300000    52.084560     1.000000    18.000000 
  [0.0814974139959,-0.996636282201,-0.00856572337417,3558.53487849] [0.991176604367,0.0819463245359,-0.104176859015,-1.18713694823] [0.104528463268,0,0.994521895368,-570.876361938]  [0,0,0,1] 39014.650000 38968.380000    -63.00058     1.000000     9.000000 
  [0.0818340199461,-0.996636282201,-0.00428873925518,3553.4951174] [0.995270426814,0.0819463245359,-0.0521599128554,-51.6324692773] [0.0523359562429,0,0.998629534755,-475.589960481]  [0,0,0,1] 39227.170000 38871.850000    12.309220     1.000000     6.000000 
  [0.0819463245359,-0.996636282201,0,3552.56073262] [0.996636282201,0.0819463245359,0,-90.8439859756] [0,0,1,-379]  [0,0,0,1] 39068.530000 38951.920000    78.665920     1.000000     0.000000 
  [0.0818340199461,-0.996636282201,0.00428873925518,3546.95381987] [0.995270426814,0.0819463245359,0.0521599128554,-131.676259583] [-0.0523359562429,0,0.998629534755,-281.371226863]  [0,0,0,1] 39057.400000 38974.840000     8.182144     1.000000     3.000000 
  [0.0814974139959,-0.996636282201,0.00856572337417,3544.15688517] [0.991176604367,0.0819463245359,0.104176859015,-160.567451435] [-0.104528463268,0,0.994521895368,-182.971234751]  [0,0,0,1] 39035.180000 38782.230000    32.230680     1.000000    12.000000 
  [0.0809374292987,-0.996636282201,0.0128192294408,3543.84477494] [0.984366035744,0.0819463245359,0.155908263646,-161.275971526] [-0.15643446504,0,0.987688340595,-84.0697312034]  [0,0,0,1] 39043.290000 38865.190000    42.476760     1.000000    15.000000 
  [0.0801556007337,-0.996636282201,0.0170375988906,3542.63972101] [0.974857388239,0.0819463245359,0.207212334563,-182.680225678] [-0.207911690818,0,0.978147600734,15.0622016342]  [0,0,0,1] 38972.400000 38458.930000    54.047440     1.000000    24.000000 
  [0.0791540712387,-0.996636282201,0.021209269466,3540.90314453] [0.962676724395,0.0819463245359,0.257948450874,-168.240747201] [-0.258819045103,0,0.965925826289,114.152850024]  [0,0,0,1] 39261.530000 38814.960000    43.824540     1.000000    27.000000 
  [0.0779355859363,-0.996636282201,0.0253228069082,3539.20343897] [0.947857430563,0.0819463245359,0.307977548411,-182.310795287] [-0.309016994375,0,0.951056516295,212.930613387]  [0,0,0,1] 39145.120000 38548.910000    85.985650     1.000000    36.000000 
  [0.0765034846101,-0.996636282201,0.0293669362967,3540.28693794] [0.9304401254,0.0819463245359,0.357162500895,-166.436328861] [-0.358367949545,0,0.933580426497,311.124748739]  [0,0,0,1] 39068.230000 38945.820000    -20.39378     1.000000    39.000000 
  [0.0748616925503,-0.996636282201,0.0333305729541,3539.16020523] [0.910472548526,0.0819463245359,0.40536849579,-161.481184149] [-0.406736643076,0,0.913545457643,408.466112781]  [0,0,0,1] 38644.770000 38312.490000    28.723170     1.000000    48.000000 
  [0.0730147097947,-0.996636282201,0.0372028528279,3542.74050015] [0.888009429684,0.0819463245359,0.452463403815,-140.14705467] [-0.45399049974,0,0.891006524188,504.687899599]  [0,0,0,1] 38551.200000 38377.930000    -32.28560     1.000000    51.000000 
  [0.0709675987948,-0.996636282201,0.0409731622679,3538.60222138] [0.863112338719,0.0819463245359,0.4983181411,-119.113337126] [-0.5,0,0.866025403784,599.526371966]  [0,0,0,1] 38648.360000 38014.520000    -86.12270     1.000000    60.000000 
  [0.0687259705396,-0.996636282201,0.0446311671183,3546.16627805] [0.835849516828,0.0819463245359,0.542807022999,-74.8817386918] [-0.544639035015,0,0.838670567945,692.721584219]  [0,0,0,1] 39392.190000 39280.810000    79.419320     1.000000    63.000000 
  [0.0662959691761,-0.996636282201,0.0481668410418,3546.93628163] [0.806295689511,0.0819463245359,0.585808108577,-52.3131546671] [-0.587785252292,0,0.809016994375,784.018094761]  [0,0,0,1] 38233.800000 38121.190000    28.456480     1.000000    72.000000 
  [0.0636842551693,-0.996636282201,0.051570493002,3560.33286861] [0.774531861754,0.0819463245359,0.627203534849,2.32667005088] [-0.62932039105,0,0.777145961457,873.165666201]  [0,0,0,1] 38817.150000 37993.650000    65.808050     1.000000    75.000000 
  [0.0608979870457,-0.996636282201,0.0548327938256,3551.50050725] [0.740645096001,0.0819463245359,0.666879839828,31.4474137123] [-0.669130606359,0,0.743144825477,959.919951243]  [0,0,0,1] 38704.830000 38435.350000    43.695990     1.000000    84.000000 
  [0.0579448017726,-0.996636282201,0.0579448017726,3566.87644746] [0.704728273521,0.0819463245359,0.704728273521,98.7503580288] [-0.707106781187,0,0.707106781187,1044.04316242]  [0,0,0,1] 38836.910000 38042.740000    30.870360     1.000000    87.000000 
  [0.0548327828959,-0.996636183544,0.0608979749071,3558.7445638] [0.666879773814,0.0819463082018,0.740645022684,139.21135038] [-0.743144825477,0,0.669130606359,1125.30472386]  [0,0,0,1] 38737.180000 37888.440000    -65.01125     1.000000    96.000000 
  [0.051570493002,-0.996636282201,0.0636842551693,3577.08153097] [0.627203534849,0.0819463245359,0.774531861754,220.985939563] [-0.777145961457,0,0.62932039105,1203.48190328]  [0,0,0,1] 38045.040000 37368.970000     1.267321     1.000000    99.000000 
  [0.0481668410418,-0.996636282201,0.0662959691761,3567.41008641] [0.585808108577,0.0819463245359,0.806295689511,271.129430647] [-0.809016994375,0,0.587785252292,1278.36042244]  [0,0,0,1] 37602.840000 37077.300000    -52.20618     1.000000   108.000000 
  [0.0446311671183,-0.996636282201,0.0687259705396,3586.40901829] [0.542807022999,0.0819463245359,0.835849516828,361.164407973] [-0.838670567945,0,0.544639035015,1349.73504455]  [0,0,0,1] 39523.340000 38642.250000    -20.61330     1.000000   111.000000 
  [0.0409731622679,-0.996636282201,0.0709675987948,3579.85223744] [0.4983181411,0.0819463245359,0.863112338719,418.265980414] [-0.866025403784,0,0.5,1417.41013672]  [0,0,0,1] 36923.730000 36471.380000    84.331470     1.000000   120.000000 
  ...
    • particles.star：
      # version 30001

  data_optics

  loop_ 
  _rlnOpticsGroup #1 
  _rlnOpticsGroupName #2 
  _rlnSphericalAberration #3 
  _rlnVoltage #4 
  _rlnTomoTiltSeriesPixelSize #5 
          1 opticsGroup1     2.700000   300.000000     1.350000 
    

  # version 30001

  data_particles

  loop_ 
  _rlnTomoName #1 
  _rlnTomoParticleId #2 
  _rlnTomoManifoldIndex #3 
  _rlnCoordinateX #4 
  _rlnCoordinateY #5 
  _rlnCoordinateZ #6 
  _rlnOriginXAngst #7 
  _rlnOriginYAngst #8 
  _rlnOriginZAngst #9 
  _rlnAngleRot #10 
  _rlnAngleTilt #11 
  _rlnAnglePsi #12 
  _rlnClassNumber #13 
  _rlnRandomSubset #14 
  _rlnTomoParticleName #15 
  _rlnOpticsGroup #16 
      TS_01            1            1   644.000000  1164.000000    72.000000     0.000000     0.000000     0.000000     -8.72000   167.290000    87.570000           -1            2    TS_01/1            1 
      TS_01            2            1   624.000000  1112.000000    64.000000     0.000000     0.000000     0.000000   -130.65000   175.350000    84.290000           -1            1    TS_01/2            1 
      TS_01            3            1   676.000000  1120.000000    68.000000     0.000000     0.000000     0.000000   -117.88000   170.530000    95.890000           -1            2    TS_01/3            1 
      TS_01            4            1   696.000000  1168.000000    84.000000     0.000000     0.000000     0.000000    53.640000   158.620000    86.530000           -1            1    TS_01/4            1 
      TS_01            5            1   592.000000  1156.000000    72.000000     0.000000     0.000000     0.000000    -78.74000   169.660000    78.390000           -1            2    TS_01/5            1 
      TS_01            6            1   659.000000  1071.000000    67.000000     0.000000     0.000000     0.000000    54.340000   178.380000    -31.41000            1            1    TS_01/6            1 
      TS_01            7            1   712.000000  1076.000000    64.000000     0.000000     0.000000     0.000000    61.580000   177.850000   156.650000           -1            2    TS_01/7            1 
      TS_01            8            1   660.000000  1200.000000   100.000000     0.000000     0.000000     0.000000   106.210000   148.260000    74.990000           -1            1    TS_01/8            1 
      TS_01            9            1   612.000000  1204.000000    80.000000     0.000000     0.000000     0.000000   -118.20000   169.440000    98.700000           -1            2    TS_01/9            1 
      TS_01           10            1   570.000000  1106.000000    64.000000     0.000000     0.000000     0.000000    -19.63000   173.260000    16.530000            1            1   TS_01/10            1
      ...

  1. 数据输出：输出每个particle的FBP重建结果以及对应的CTF weight volume。

     |-Subtomograms
     |-1_data.mrc
     |-1_weight.mrc
     |-2_data.mrc
     |-2_weight.mrc
     |-......

  2. Program detail：

     void SubtomoProgram::run()
 {
     //输入参数
     if (cropSize < 0) cropSize = boxSize;
        
     //默认需要ctf
     bool do_ctf = true;

     //s2D是指每个particle对应的2D crop image的size
     const long int s2D = boxSize; 
        
     //s3D是指最后输出的particle volume boxsize
     const long int s3D = cropSize; 
     //sh3D是ctf weight volume的x轴长度，因为relion的fft是rfft，x轴只有正半轴
     const long int sh3D = s3D / 2 + 1; 
        
     //s02D指的是在bin1rawdata上s2D size对应的大小
     const long int s02D = (int)(binning * s2D + 0.5); 
        
     const double relative_box_scale = cropSize / (double) boxSize;
     const double binned_pixel_size = binning * particleSet.getOriginalPixelSize(0);
     ...
 }

    遍历particles，投影到各个view上crop出boxSize大小的2D image。

     void SubtomoProgram::processTomograms(
         const std::vector<int>& tomoIndices,
         const TomogramSet& tomogramSet,
         const ParticleSet& particleSet,
         const std::vector<std::vector<ParticleIndex>>& particles,
         const AberrationsCache& aberrationsCache,
         long int s02D,
         long int s2D,
         long int s3D,
         double relative_box_scale,
         bool do_ctf,
         int verbosity,
         BufferedImage<float>& sum_data,
         BufferedImage<float>& sum_weights )
 {
     ...
     #pragma omp parallel for num_threads(outer_thread_num)
     for (int p = 0; p < pc; p++){
         //particle对应3D坐标
         const std::vector<d3Vector> traj = particleSet.getTrajectoryInPixels(
                         part_id, fc, tomogram.optics.pixelSize, !apply_offsets);

         if (!tomogram.isVisibleAtAll(traj, s2D / 2.0))
         {
             continue;
         }

         //判断particle可不可以被每个view看到，每个view判断一次
         const std::vector<bool> isVisible = tomogram.determineVisiblity(traj, s2D / 2.0);

         std::vector<d4Matrix> projCut(fc), projPart(fc);

         //2D图片数据（用于FBP）
         BufferedImage<fComplex> particleStack = BufferedImage<fComplex>(sh2D,s2D,fc);
         //CTF volume
         BufferedImage<float> weightStack(sh2D,s2D,fc);
         //crop 2D image的函数
         TomoExtraction::extractAt3D_Fourier(
                     tomogram.stack, s02D, binning, tomogram, traj, isVisible,
                     particleStack, projCut, inner_thread_num, do_circle_precrop);
         ...
     }
     ...
 }

    下方函数源码位置:src/jaz/tomography/extraction.h

     template <typename T>
 void TomoExtraction::extractAt3D_Fourier(
         const RawImage<T>& stack, int s, double bin,
         const Tomogram& tomogram,
         const std::vector<gravis::d3Vector>& trajectory,
         const std::vector<bool>& isVisible,
         RawImage<tComplex<T>>& out,
         std::vector<gravis::d4Matrix>& projOut,
         int num_threads,
         bool circle_crop)
 {
     const int fc = tomogram.frameCount;
     std::vector<gravis::d2Vector> centers(fc);
        
     for (int f = 0; f < fc; f++)
     {
         //Apply投影矩阵获得3D坐标到图片上的像素位置
         centers[f] = tomogram.projectPoint(trajectory[f], f);
     }
        
         //根据2Dcenter crop image
     extractAt2D_Fourier(
         stack, s, bin, tomogram.projectionMatrices, centers, isVisible,
         out, projOut, num_threads, circle_crop);
 }

 template <typename T>
 void TomoExtraction::extractAt2D_Fourier(
         const RawImage<T>& stack, int s, double bin,
         const std::vector<gravis::d4Matrix>& projIn,
         const std::vector<gravis::d2Vector>& centers,
         const std::vector<bool>& isVisible,
         RawImage<tComplex<T>>& out,
         std::vector<gravis::d4Matrix>& projOut,
         int num_threads,
         bool circle_crop)
 {
     const int sh = s/2 + 1;
     const int fc = stack.zdim;

     //如果有bin的需求，那么会在这里bin数据
     const int sb = (int)(s / bin + 0.5);
     const int sbh = sb/2 + 1;
        
     BufferedImage<T> smallStack(s,s,fc);
     projOut.resize(fc);
                
     std::vector<gravis::d2Vector> integralShift(fc);
                
     for (int f = 0; f < fc; f++)
     {
         integralShift[f] = gravis::d2Vector(
                 round(centers[f].x) - s/2,
                 round(centers[f].y) - s/2);
     }
        
     //简单的crop,如果 not visible，就全填0
     extractSquares(stack, s, s, integralShift, isVisible, smallStack, false, num_threads);
        
     if (circle_crop) //默认关闭
     {
         //在ctf前提前在原数据上crop circle
         cropCircle(smallStack, 0, EDGE_FALLOFF, num_threads);
     }
        
     std::vector<gravis::d2Vector> posInNewImg(fc);

     //根据crop位置调整每张图的投影矩阵
     for (int f = 0; f < fc; f++)
     {
         projOut[f] = projIn[f];
            
         projOut[f](0,3) += sb/2 - centers[f].x;
         projOut[f](1,3) += sb/2 - centers[f].y;
            
         posInNewImg[f] = (centers[f] - integralShift[f]) / bin;
     }	
        
     BufferedImage<tComplex<T>> smallStackFS(sh,s,fc);

     NewStackHelper::FourierTransformStack_fast(smallStack, smallStackFS, true, num_threads);

     if (bin != 1.0)
     {
         //频域apply binning
         smallStackFS = Resampling::FourierCrop_fftwHalfStack(
                 smallStackFS, bin, num_threads);

         //对应的数据scale
         smallStackFS /= (T) bin;
     }

     //返回频域数据
     NewStackHelper::shiftStack(smallStackFS, posInNewImg, out, true, num_threads);
 }

    回到subtomo.cpp

                 for (int f = 0; f < fc; f++)
             {
                 //看不见的view全是0，不用管
                 if (!isVisible[f]) continue;

                 d3Matrix A;

                 //计算particle自带的旋转（如果需要的话），默认不需要，即A是Identity矩阵
                 if (apply_orientations)//是否要在重建的particle上apply欧拉角旋转,默认关闭
                 {
                     A = particleSet.getMatrix3x3(part_id);
                 }
                 else
                 {
                     A = particleSet.getSubtomogramMatrix(part_id);
                 }
                    
                 projPart[f] = projCut[f] * d4Matrix(A);
                 if (do_ctf)//默认开启
                 {
                     const d3Vector pos = (apply_offsets) ? particleSet.getPosition(part_id) : particleSet.getParticleCoord(part_id);

                     //根据输入的ctf参数计算每个2Dcropimage对应的频域ctf image
                     CTF ctf = tomogram.getCtf(f, pos);
                     BufferedImage<float> ctfImg(sh2D, s2D);
                     ctf.draw(s2D, s2D, binnedPixelSize, gammaOffset, &ctfImg(0,0,0));

                     //默认开启，即把particle从黑的变成白的
                     const float sign = flip_value? -1.f : 1.f;
                     for (int y = 0; y < s2D;  y++)
                     for (int x = 0; x < xRanges(y,f); x++)
                     {
                         const double c = ctfImg(x,y) * doseWeights(x,y,f);

                         particleStack(x,y,f) *= sign * c;
                         weightStack(x,y,f) = c * c;
                     }
                 }
             }
             ...
             //根据boxsize与cropsize做circle crop
             const int boundary = (boxSize - cropSize) / 2;
             if (do_gridding_precorrection || do_circle_crop)
             {
                 BufferedImage<float> particlesRS;

                 //在时域上做
                 particlesRS = NewStackHelper::inverseFourierTransformStack(particleStack);
                 if (do_circle_crop)//默认开启
                 {
                     const double crop_boundary = do_narrow_circle_crop? boundary : 0.0;
                     TomoExtraction::cropCircle(particlesRS, crop_boundary, 5, num_threads);//在boundary边缘5个pixel进行平滑处理，再外面全设为0
                 }

                 if (do_gridding_precorrection)//默认关闭
                 {
                     TomoExtraction::griddingPreCorrect(particlesRS, boundary, num_threads);
                 }
                 particleStack = NewStackHelper::FourierTransformStack(particlesRS);
             }

             //开始做FBP
             BufferedImage<fComplex> dataImgFS(sh3D,s3D,s3D);
             dataImgFS.fill(fComplex(0.0, 0.0));

             BufferedImage<float> ctfImgFS(sh3D,s3D,s3D),
                     dataImgRS(s3D,s3D,s3D), dataImgDivRS(s3D,s3D,s3D),
                     multiImageFS(sh3D,s3D,s3D);

             ctfImgFS.fill(0.0);
             dataImgRS.fill(0.0);
             dataImgDivRS.fill(0.0);

             for (int f = 0; f < fc; f++)
             {
                 if (isVisible[f])
                 {
                     //FBP时会将旋转后超过cropsize的部分扔掉，即在3D频域crop volume
                     FourierBackprojection::backprojectSlice_forward_with_multiplicity(
                         &xRanges(0,f),
                         particleStack.getSliceRef(f),
                         weightStack.getSliceRef(f),
                         projPart[f] * relative_box_scale,//投影矩阵根据cropsize与boxsize做scale，relative_box_scale = crop_size/bos_size.
                         dataImgFS,
                         ctfImgFS,
                         multiImageFS);
                        
                 }
             }

             Centering::shiftInSitu(dataImgFS);

             //因为有频域crop，需要调整数据scale
             if (s3D != s2D)
             {
                 dataImgFS *= (float) sqrt(s2D / (double) s3D);
             }

             FFT::inverseFourierTransform(dataImgFS, dataImgRS, FFT::Both);

             //保存结果
             dataImgRS.write(outData, binnedPixelSize, write_float16);
             if (write_ctf)
             {
                 Centering::fftwHalfToHumanFull(ctfImgFS).write(outCTF, 1.0 / binnedPixelSize, write_float16);
             }
             ...

    重建完后下一步就是3D initial model与autorefine。

Git 教程
GitHub 教程
下一页 →

相关文档

基础知识
数据处理案例