Merge pull request biocore#29 from cameronmartino/fix_issues

Re-centering for trajectories and Rank-2 fix.

gemelli/ctf.py

@@ -134,10 +134,10 @@ def ctf_helper(table: biom.Table,
     # <https://github.com/biocore/emperor/commit
     # /a93f029548c421cb0ba365b4294f7a5a6b0209ce>
     if n_components == 2:
-        TF.subjects['PC3'] = [0] * len(TF.subjects.index)
-        TF.features['PC3'] = [0] * len(TF.features.index)
-        TF.proportion_explained.loc['PC3', :] = 0
-        TF.eigvals.loc['PC3', :] = 0
+        TF.subjects.loc[:, 'PC3'] = [0] * len(TF.subjects.index)
+        TF.features.loc[:, 'PC3'] = [0] * len(TF.features.index)
+        TF.proportion_explained['PC3'] = 0
+        TF.eigvals['PC3'] = 0
 
     # save ordination results
     short_method_name = 'CTF_Biplot'
@@ -190,6 +190,10 @@ def ctf_helper(table: biom.Table,
         # ensure index name for q2
         straj.index.name = "#SampleID"
         # save traj.
+        keep_PC_traj = [col for col in straj.columns
+                        if 'PC' in col]
+        straj[keep_PC_traj] -= straj[keep_PC_traj].mean()
+        ftraj[keep_PC_traj] -= ftraj[keep_PC_traj].mean()
         subject_trajectories[condition] = straj
         ftraj.index = ftraj.index.astype(str)
         feature_trajectories[condition] = ftraj

gemelli/factorization.py

@@ -230,19 +230,19 @@ def _fit(self):
             self.biplot_components = np.min(possible_comp)
             X = X - X.mean(axis=0)
             X = X - X.mean(axis=1).reshape(-1, 1)
-            u, s, v = svd(X)
+            u, s, v = svd(X, full_matrices=False)
             u = u[:, :self.biplot_components]
             v = v.T[:, :self.biplot_components]
-            p = s**2 / np.sum(s**2)
+            p = s * (1 / s.sum())
             p = np.array(p[:self.biplot_components])
             s = np.diag(s[:self.biplot_components])
             # save the re-centered biplot
             self.features = v
             self.subjects = u
         else:
             # just make prop-exp
-            p = np.array(np.diag(s)**2 /
-                         np.sum(np.diag(s)**2))
+            p = np.array(np.diag(s) /
+                         1 / np.sum(np.diag(s)))
         # save all eigen values
         self.eigvals = np.diag(s)
         # the proortion explained for n_components
@@ -265,7 +265,7 @@ def _fit(self):
             # temporary list of components in feature trajectory
             feature_temp_trajectory = []
             # for each component in the rank given to TensorFactorization
-            for component in range(self.n_components):
+            for component in range(self.n_components)[::-1]:
                 # component condition-subject trajectory
                 dtmp = np.dot(loads[0][:, [component]],
                               condition[:, [component]].T).flatten()
@@ -277,6 +277,9 @@ def _fit(self):
             # combine all n_components
             subject_temp_trajectory = np.array(subject_temp_trajectory).T
             feature_temp_trajectory = np.array(feature_temp_trajectory).T
+            # double check check centered
+            subject_temp_trajectory -= subject_temp_trajectory.mean(axis=0)
+            feature_temp_trajectory -= feature_temp_trajectory.mean(axis=0)
             # save subject-condition trajectory and distance matrix
             self.subject_trajectory.append(subject_temp_trajectory)
             self.subject_distances.append(

gemelli/q2/tests/test_method.py

@@ -68,6 +68,7 @@ def setUp(self):
         self.q2meta = Metadata(self.meta_table)
         self.subj = 'host_subject_id'
         self.state = 'context'
+        self.out_ = os_path_sep.join(self.in_table.split(os_path_sep)[:-1])
 
     def test_qiime2_ctf(self):
         """Tests that the Q2 and standalone ctf results match.
@@ -90,7 +91,6 @@ def test_qiime2_ctf(self):
         # ...First, though, we need to write the contents of self.q2table to a
         # BIOM file, so gemelli can understand it.
         # Derived from a line in test_standalone_ctf()
-        out_ = os_path_sep.join(self.in_table.split(os_path_sep)[:-1])
         # Run gemelli outside of QIIME 2...
         runner = CliRunner()
         result = runner.invoke(standalone_ctf,
@@ -103,7 +103,7 @@ def test_qiime2_ctf(self):
                                 '--state-column-1',
                                 'context',
                                 '--output-dir',
-                                out_])
+                                self.out_])
         # check exit code was 0 (indicating success)
         self.assertEqual(result.exit_code, 0)
         # ...and read in the resulting output files. This code was derived from
@@ -125,6 +125,18 @@ def test_qiime2_ctf(self):
                         absolute_sort(q2ftraj[comp_col].values),
                         atol=.5)
 
+    def test_ctf_rank2(self):
+        """Tests that ctf with rank < 3
+        """
+        # Run gemelli
+        res = q2gemelli.actions.ctf(table=self.q2table,
+                                    sample_metadata=self.q2meta,
+                                    individual_id_column=self.subj,
+                                    state_column=self.state,
+                                    n_components=2)
+        # check exit code was 0 (indicating success)
+        self.assertEqual(len(res), 5)
+
 
 if __name__ == "__main__":
     unittest.main()

gemelli/scripts/tests/data/expected-context-distance-matrix.tsv

@@ -1,10 +1,10 @@
-	sample_261	sample_466	sample_434	sample_234	sample_67	sample_140	sample_165	sample_55	sample_111
-sample_261	0.0	67.91465121470452	55.487976558660314	54.17097763707955	71.8869053567693	64.03065302962949	63.26084447452229	55.46972695955536	65.37625367108693
-sample_466	67.91465121470452	0.0	16.89836551768841	18.515833230343233	5.792857602618866	24.140045672118926	41.243368307630504	22.430254689969164	23.31408177994854
-sample_434	55.487976558660314	16.89836551768841	0.0	20.853948091617333	21.37930730220944	13.836844131466576	27.641068078011795	7.10793622006	14.00252658225181
+	sample_261	sample_527	sample_434	sample_234	sample_218	sample_140	sample_477	sample_325	sample_201
+sample_261	0.0	67.91465121470452	55.48797655866031	54.17097763707955	71.88690535676929	64.03065302962949	63.26084447452228	55.46972695955536	65.37625367108693
+sample_527	67.91465121470452	0.0	16.89836551768841	18.515833230343233	5.792857602618866	24.140045672118926	41.243368307630504	22.430254689969164	23.31408177994854
+sample_434	55.48797655866031	16.89836551768841	0.0	20.853948091617333	21.379307302209437	13.836844131466576	27.641068078011795	7.10793622006	14.00252658225181
 sample_234	54.17097763707955	18.515833230343233	20.853948091617333	0.0	20.46729361418852	34.211896461946814	47.1955237836692	27.504836133965814	34.027908559260226
-sample_67	71.8869053567693	5.792857602618866	21.37930730220944	20.46729361418852	0.0	28.07049490096316	44.159191560420325	27.11028036731012	26.850820924757976
+sample_218	71.88690535676929	5.792857602618866	21.379307302209437	20.46729361418852	0.0	28.07049490096316	44.159191560420325	27.110280367310125	26.850820924757976
 sample_140	64.03065302962949	24.140045672118926	13.836844131466576	34.211896461946814	28.07049490096316	0.0	19.747932075060884	9.243704807547592	2.493866919694757
-sample_165	63.26084447452229	41.243368307630504	27.641068078011795	47.1955237836692	44.159191560420325	19.747932075060884	0.0	22.717764881650552	19.75348391748244
-sample_55	55.46972695955536	22.430254689969164	7.10793622006	27.504836133965814	27.11028036731012	9.243704807547592	22.717764881650552	0.0	10.463563482783073
-sample_111	65.37625367108693	23.31408177994854	14.00252658225181	34.027908559260226	26.850820924757976	2.493866919694757	19.75348391748244	10.463563482783073	0.0
+sample_477	63.26084447452228	41.243368307630504	27.641068078011795	47.1955237836692	44.159191560420325	19.747932075060884	0.0	22.717764881650552	19.75348391748244
+sample_325	55.46972695955536	22.430254689969164	7.10793622006	27.504836133965814	27.110280367310125	9.243704807547592	22.717764881650552	0.0	10.463563482783073
+sample_201	65.37625367108693	23.31408177994854	14.00252658225181	34.027908559260226	26.850820924757976	2.493866919694757	19.75348391748244	10.463563482783073	0.0