<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Compare the read depth and number of strains
This data is the average read depth of each metagenome. The table in read_depth.strains.tsv has the SRA ID, the average read depth across the amplicon region, and the number of strains that we recovered for each of the amplicon regions, A, B, and C.
Step1:  Filter the data
Because (0,0) correlates strongly, we filter out any row where the read_depth and the strain count are 0. We have to do this on a per SRR basis.
Step2:  Note that we have reduced our matrix from having 11,054 entries with all the zeros to only having 1,397 entries now!
plot the data.
This just provides a quick overview of the data
Step3:  Note that this plot is skewed by a few outliers. Lets limit it to anything where read_depth < 1000 and redraw the plot
Step4:  When we zoom in, this doesn't look like a strong correlation. Note that there are a lot of data points here compared to the whole data set. In the data set excluding (0,0) we had 1,397 entries, and now we have 1,386 entries, so we only removed 9 values!
Linear regression
What is the correlation between these two data sets. Note that we use the data with all the non-zero's removed (1,397 data points).
As a reminder, the statsmodels OLS uses y ~ x
Step5:  Removing the outliers
Remember the nine outliers above? If we use the data set where we have removed them, we can compare the r<sup>2</sup> value for that dataset.
<ASSISTANT_TASK:>
Python Code:
#instantiate our environment
import os
import sys
%matplotlib inline
import pandas as pd
import statsmodels.api as sm
# read the data into a pandas dataframe
df = pd.read_csv("read_depth.strains.tsv", header=0, delimiter="\t")
print("Shape: {}".format(df.shape))
df.head()
Explanation: Compare the read depth and number of strains
This data is the average read depth of each metagenome. The table in read_depth.strains.tsv has the SRA ID, the average read depth across the amplicon region, and the number of strains that we recovered for each of the amplicon regions, A, B, and C.
End of explanation
dfa = df[(df["A_read_depth"] > 0) & (df["A_strains"] > 0)]
dfb = df[(df["B_read_depth"] > 0) & (df["B_strains"] > 0)]
dfc = df[(df["C_read_depth"] > 0) & (df["C_strains"] > 0)]
print("Shape: {}".format(dfa.shape))
dfa.head()
Explanation: Filter the data
Because (0,0) correlates strongly, we filter out any row where the read_depth and the strain count are 0. We have to do this on a per SRR basis.
End of explanation
ax = dfa.plot('A_read_depth', 'A_strains', kind='scatter')
ax.set(ylabel="# strains", xlabel="read depth")
Explanation: Note that we have reduced our matrix from having 11,054 entries with all the zeros to only having 1,397 entries now!
plot the data.
This just provides a quick overview of the data
End of explanation
dfas = dfa[dfa['A_read_depth'] < 1000]
print("Shape: {}".format(dfas.shape))
ax = dfas.plot('A_read_depth', 'A_strains', kind='scatter')
ax.set(ylabel="# strains", xlabel="read depth")
Explanation: Note that this plot is skewed by a few outliers. Lets limit it to anything where read_depth < 1000 and redraw the plot
End of explanation
model = sm.OLS(dfa['A_strains'], dfa['A_read_depth']).fit()
predictions = model.predict(dfa['A_read_depth'])
model.summary()
Explanation: When we zoom in, this doesn't look like a strong correlation. Note that there are a lot of data points here compared to the whole data set. In the data set excluding (0,0) we had 1,397 entries, and now we have 1,386 entries, so we only removed 9 values!
Linear regression
What is the correlation between these two data sets. Note that we use the data with all the non-zero's removed (1,397 data points).
As a reminder, the statsmodels OLS uses y ~ x
End of explanation
model = sm.OLS(dfas['A_strains'], dfas['A_read_depth']).fit()
predictions = model.predict(dfas['A_read_depth'])
model.summary()
Explanation: Removing the outliers
Remember the nine outliers above? If we use the data set where we have removed them, we can compare the r<sup>2</sup> value for that dataset.
End of explanation
<END_TASK>