<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Utility functions to parse experiment metadata
Step1:  Read the counts table
This is bacterial mRNA-Seq with samples at 4 different temperatures with or without the addition of BCM. Each condition is sequenced in triplicates. And we are interested in 5'UTRs transcription levels.
Step2:  Normalize to UTR length
Step3:  Notation
X 
Step4:  Principal Component Analisys (PCA)
Linear dimensionality reduction using Singular Value Decomposition of the data to project it to a lower dimensional space. The input data is centered but not scaled for each feature before applying the SVD.
The resulting components are ranked by the amount of variance explained in the data.
Step5:  Aside
Step6:  +BCM samples
<ASSISTANT_TASK:>
Python Code:
def parse_barcodes(bcfile, bc_id='BC'):
    res = {}
    with open(bcfile, 'r') as fi:
        for line in fi:
            fields = line.strip().split(',')
            if fields[0].startswith(bc_id):
                res[fields[0]] = fields[1]
    return res
def parse_exp_config(expfile, bc_dict):
    res = []
    fieldnames = ['id', 'sample', 'cond', 'barcode', 'size', 'region', 'Qbit', 'conc', 'dilution']
    with open(expfile) as fi:
        reader = csv.DictReader(fi, fieldnames=fieldnames)
        for rec in reader:
            if rec['id']:
                res.append({
                    'sample': rec['sample'],
                    'bc_id': rec['barcode'],
                    'bc_seq': bc_dict[rec['barcode']],
                    'temp': int(rec['cond'][:2]),
                    'bcm': '+' in rec['cond'],
                })
    return pd.DataFrame.from_records(res)
Explanation: Utility functions to parse experiment metadata
End of explanation
bc_dict = parse_barcodes('../../data/Lexogen_Sense_RNA-Seq.csv')
exp_df = parse_exp_config('../../data/2017-03-09_NextSeq.csv', bc_dict)
agg_utr = pd.read_csv('../../data/utr.counts.csv')
agg_utr
Explanation: Read the counts table
This is bacterial mRNA-Seq with samples at 4 different temperatures with or without the addition of BCM. Each condition is sequenced in triplicates. And we are interested in 5'UTRs transcription levels.
End of explanation
def normalize(df, edf, columns=None):
    '''
    Prepares the UTR dataframe (`df`) for log transformation.
    
    Adds experiment metadata from `edf`.
    Adds pseudocounts to `utr_counts` and `UTR_length`.
    Normalizes counts to UTR length.
    '''
    def pseudo_counts(x):
        return x + 1 if x == 0 else x
    df = df.merge(edf, how='left', on='sample')
    if columns is not None:
        df = df[columns]
    # Add pseudocounts to allow log transform later
    df['utr_counts'] = df['utr_counts'].apply(pseudo_counts)
    df['UTR_length'] = df['UTR_length'].apply(pseudo_counts)
    df['utr_norm'] = df['utr_counts'] / df['UTR_length']
    return df
columns = ['gene', 'TSS', 'start', 'end', 'UTR_length',
               'utr_counts', 'sample', 'bcm', 'temp']
utr = normalize(agg_utr, exp_df, columns)
utr
Explanation: Normalize to UTR length
End of explanation
# build expression matrix
X = pd.DataFrame()
samples = []
for sample in set(utr['sample']):
    mask = (utr['sample']==sample) & (utr['bcm']==False)
    if not utr[mask].empty:
        X[sample] = utr[mask]['utr_norm'].values
        samples.append(sample)
    
# Same as .fit() and then .transform()
X_std = StandardScaler().fit_transform(X.values.T)
X_std
Explanation: Notation
X : data (features)
Y : predicted value
-BCM samples
From scikit-learn docs:
Standardize features by removing the mean and scaling to unit variance
The standard score of a sample x is calculated as:
z = (x - u) / s
where u is the mean of the training samples or zero if with_mean=False, and s is the standard deviation of the training samples or one if with_std=False.
End of explanation
sklearn_pca = sklearnPCA(n_components=10)
Y = sklearn_pca.fit_transform(X_std)
print(Y)
print(sklearn_pca.explained_variance_)
print(sklearn_pca.explained_variance_ratio_)
Explanation: Principal Component Analisys (PCA)
Linear dimensionality reduction using Singular Value Decomposition of the data to project it to a lower dimensional space. The input data is centered but not scaled for each feature before applying the SVD.
The resulting components are ranked by the amount of variance explained in the data.
End of explanation
vdf = pd.DataFrame()
vdf['PC'] = [(i+1) for i,x in enumerate(sklearn_pca.explained_variance_ratio_)]
vdf['var'] = sklearn_pca.explained_variance_ratio_
(ggplot(vdf, aes(x='PC', y='var'))
    + geom_point(size=5, alpha=0.3)
    + ylab('Explained variance')
    + ggtitle('Unfiltered -BCM')
)
pca_df = pd.DataFrame()
pca_df['cond'] = ['%doC' % exp_df[exp_df['sample']==sample]['temp'] for sample in samples]
pca_df['PC1'] = Y[:,0]
pca_df['PC2'] = Y[:,1]
pca_df
(ggplot(pca_df, aes(x='PC1', y='PC2', color='cond'))
    + geom_point(size=5, alpha=0.3)
    + ggtitle('Unfiltered -BCM')
)
Explanation: Aside: ggplot terminology
from: https://beanumber.github.io/sds192/lab-ggplot2.html
Geometric Objects (geom)
Geometric objects or geoms are the actual marks we put on a plot. Examples include:
points (geom_point, for scatter plots, dot plots, etc)
lines (geom_line, for time series, trend lines, etc)
boxplot (geom_boxplot, for, well, boxplots!)
… and many more!
A plot should have at least one geom, but there is no upper limit. You can add a geom to a plot using the + operator.
Aesthetic Mapping (aes)
In ggplot2, aesthetic means “something you can see”. Each aesthetic is a mapping between a visual cue and a variable. Examples include:
position (i.e., on the x and y axes)
color (“outside” color)
fill (“inside” color)
shape (of points)
line type
size
Each type of geom accepts only a subset of all aesthetics—refer to the geom help pages to see what mappings each geom accepts. Aesthetic mappings are set with the aes() function.
End of explanation
# build expression matrix
X = pd.DataFrame()
samples = []
for sample in set(utr['sample']):
    mask = (utr['sample']==sample) & (utr['bcm']==True)
    if not utr[mask].empty:
        X[sample] = utr[mask]['utr_norm'].values
        samples.append(sample)
    
X_std = StandardScaler().fit_transform(X.values.T)
X_std
sklearn_pca = sklearnPCA(n_components=10)
Y = sklearn_pca.fit_transform(X_std)
print(Y)
print(sklearn_pca.explained_variance_)
print(sklearn_pca.explained_variance_ratio_)
vdf = pd.DataFrame()
vdf['PC'] = [(i+1) for i,x in enumerate(sklearn_pca.explained_variance_ratio_)]
vdf['var'] = sklearn_pca.explained_variance_ratio_
(ggplot(vdf, aes(x='PC', y='var'))
    + geom_point(size=5, alpha=0.3)
    + ylab('Explained variance')
    + ggtitle('Unfiltered +BCM')
)
pca_df = pd.DataFrame()
pca_df['cond'] = ['%doC' % exp_df[exp_df['sample']==sample]['temp'] for sample in samples]
pca_df['PC1'] = Y[:,0]
pca_df['PC2'] = Y[:,1]
pca_df
(ggplot(pca_df, aes(x='PC1', y='PC2', color='cond'))
    + geom_point(size=5, alpha=0.3)
    + ggtitle('Unfiltered +BCM')
)
Explanation: +BCM samples
End of explanation
<END_TASK>