library("tidyr")
library("ggpubr")
library("stringr")
library("plyr")
library("dplyr")Bioinformatics Class Project Day 3
1 Introduction
Last week, we set out to test the hypothesis that TEs are contributing to the evolution of the worker-queen differentiation in eusocial shrimp. We began with raw RNAseq data from three workers and three queens. I have assembled a transcriptome from the RNAseq data and ran Repeat Masker on the transcriptome assembly. We analyzed the raw RNAseq data with Galaxy to find the differential expressed genes and Finally, we loaded and explored two data sets in R:
df- List of transcripts that are SHARED in all individuals with results related to differential expression between queens and workers (mainlylog2FoldChangeandpadj).df_transcript_te- List of transposable elements (TEs) and their locations in the transcripts that has TEs.
Our task today is to combine these data sets into a master data set and to decide what analyses can test whether the data support the hypothesis.
2 Load libraries
3 Make a master data set
How do we combine df and df_transcript_te? Examine these files to see if there’s common column.
Note that df has only transcripts that are mapped across all individuals, while df_transcript_te has only transcripts that have TEs. None of these files have ALL transcripts in the assembled transcriptome. But we do not need that information.
Write some codes to combine df and df_transcript_te and save the it as df_tranTE. Look at the different ways to combine dataframes using dplyr here, decide which one should be used.
# your code
# This should return true if you did this correctly
# nrow(df_tranTE) == nrow(df)4 Discussion
If the hypothesis is true, what pattern would you expect to see in
df_tranTE? These are your predictions/analysis.What are the predictor and response in each analysis?
What statistic tests should be used for each analysis?
5 Predictions
Below, we will run through a few example analysis using all TEs.
You do not need to report Prediction 1 in your project report. Instead, I’ll provide a sample report based on Prediction 1.
For Prediction 2 and Prediction 3, you will need to
(1) report the main analysis here (with all TEs) and
(2) modify the codes to run it for your assigned TE class and report it in the project report.
Name TE_class Akhila DNA Diandra DNA Mohab LTR Rebeca LTR Sanjna LINE Steven LINE Tyla SINE Solomon SINE
Prediction 1
If TEs are contributing to the evolution of the worker-queen differentiation (our hypothesis), there would be more TEs in DE genes than non-DE genes.
What are the predictor & response variables?
Predictor = DE (Y or N, categorical)
Response = TE (Y or N, categorical)
What statistic tests should be used?
- Chi-sq test of independence
To do this, we will need to create two new columns DE and TE in df_tranTE that indicated whether the gene is DE or not, and has TE or not. We can use mutate in combination with case_when.
df_tranTE = mutate(df_tranTE,
DE = case_when(padj < 0.05 ~ 'DE',
padj >= 0.05 ~ 'NS'),
TE = case_when(TE_class %in% c("DNA",
"LTR",
"LINE",
"SINE") ~ 'TE',
.default = 'Not')
)
# In the second case_when(), when the TE_class is not the ones in c(), the cell will be assigned "Not"Then we can use table to make a contingency table that has the count of transcripts in all combinations of the columns DE and TE.
obs_table = table(df_tranTE$TE, df_tranTE$DE)
obs_table
DE NS
Not 224 26108
TE 39 3462# Your codeRemember that for a chi-square test of independence, the null hypothesis is that the two variables (TE and DE) are independent. When p < 0.05, we can reject the null hypothesis. When p > 0.05, we fail to reject the null hypothesis.
What does the test result tell you?
Let’s plot this using the raw data df_tranTE.
ggplot(data = df_tranTE, aes(x = TE, fill = TE)) +
geom_bar() +
facet_grid(DE~., scales = "free_y") +
labs(x = "Transposable Elements", y = "Count")
We can further test whether TE_class affect whether a transcript is DE or not.
chisq.test(with(df_tranTE, table(TE_class, DE)))Warning in chisq.test(with(df_tranTE, table(TE_class, DE))): Chi-squared
approximation may be incorrect
Pearson's Chi-squared test
data: with(df_tranTE, table(TE_class, DE))
X-squared = 1.2042, df = 3, p-value = 0.752What does this result tell you?
Summary Statement
In summary we can say that (note the use of 3 significant figures in reporting stats.):
The proportions of transcripts having TEs or not are not statistically different between transcripts that are and are not differentially expressed between queens and workers (Chi-squared = 2.160, d.f. = 1, p = 0.142). Further, whether a transcript is differentially expressed or not is not dependent on the major TE class (DNA, LTR, LINE, and SINE) (Chi-squared = 1.204, d.f. = 3, p = 0.752). These results do not support our hypothesis that TEs are contributing to the evolution of the worker-queen differentiation.