<jupyter_start><jupyter_text>PCA Python vs ROriginally, R was used to calculate PCA using both `princomp` and `prcomp`. However, rpy2 stopped was intorducing some issues on the galaxy server. I decided to switch the calculation over to a pure python solution. scikit-learn has a PCA package which we can used, but it only does SVD and matches the output of `prcomp` with its default values. Here I am testing and figuing out how to output the different values.<jupyter_code>import pandas as pd
import numpy as np
from sklearn.decomposition import PCA<jupyter_output><empty_output><jupyter_text>Import Example Data<jupyter_code>dat = pd.read_table('../example_data/ST000015_log.tsv')
dat.set_index('Name', inplace=True)
dat[:3]<jupyter_output><empty_output><jupyter_text>Use R to calculate PCAMi has looked at this already, but wanted to put the R example here to be complete. Here are the two R methods to output PCA<jupyter_code>%%R -i dat
# First method uses princomp to calulate PCA using eigenvalues and eigenvectors
pr = princomp(dat)
#str(pr)
loadings = pr$loadings
scores = pr$scores
#summary(pr)
%%R -i dat
pr = prcomp(dat)
#str(pr)
loadings = pr$rotation
scores = pr$x
sd = pr$sdev
#summary(pr)<jupyter_output><empty_output><jupyter_text>Use Python to calculate PCAscikit-learn has a PCA package that we will use. It uses the SVD method, so results match the prcomp from R. http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html Generate PCA with default settings<jupyter_code># Initiate PCA class
pca = PCA()

# Fit the model and transform data
scores = pca.fit_transform(dat)

# Get loadings
loadings = pca.components_

# R also outputs the following in their summaries
sd = loadings.std(axis=0)
propVar = pca.explained_variance_ratio_
cumPropVar = propVar.cumsum()<jupyter_output><empty_output><jupyter_text>I compared these results with prcomp and they are identical, note that the python version formats the data in scientific notation. Build output tables that match the original PCA script Build comment blockAt the top of each output file, the original R version includes the standard deviation and the proportion of variance explained. I want to first build this block.<jupyter_code># Labels used for the comment block
labels = np.array(['#Std. deviation', '#Proportion of variance explained', '#Cumulative proportion of variance explained'])

# Stack the data into a matrix
data = np.vstack([sd, propVar, cumPropVar])

# Add the labels to the first position in the matrix
block = np.column_stack([labels, data])
# Create header
header = np.array(['Comp{}'.format(x+1) for x in range(loadings.shape[1])])
compoundIndex = np.hstack([dat.index.name, dat.index])
sampleIndex = np.hstack(['sampleID', dat.columns])

# Create loadings output
loadHead = np.vstack([header, loadings])
loadIndex = np.column_stack([sampleIndex, loadHead])
loadOut = np.vstack([block, loadIndex])

# Create scores output
scoreHead = np.vstack([header, scores])
scoreIndex = np.column_stack([compoundIndex, scoreHead])
scoreOut = np.vstack([block, scoreIndex])
np.savetxt('/home/jfear/tmp/dan.tsv', loadOut, fmt="%s", delimiter='\t')
bob = pd.DataFrame(loadOut)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09541627689429373, 'fraction_numerical': 0.0034299968818210166, 'mean_word_length': 5.840085287846482, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '413295', 'license': 'MIT', 'n_tokens_mistral': 1009, 'n_tokens_neox': 974, 'n_words': 360, 'path': 'notebook/pca_testing.ipynb', 'repo_name': 'GalaxyDream/SECIMTools', 'size': '3207'}