Transcriptional Changes Fade Prior to Long-Term Memory for Sensitization of the <i>Aplysia</i> Siphon-Withdrawal Reflex.

1. #Analysis of LTS microarray data, 5 days after training
2. # v1.1 - Bob Calin-Jageman, 3/3/2017
3. # v1.1 - 3/3/2017 -- initial version for pre-registration, based on prior script for 1day analysis
6. #Notes
7. #This r script analyses differential expression for trained vs. untrained samples of pleural ganglia taken from
8. # Aplysia 5 days after LTS training (4 shocks, 90mA of 0.5s on 0.5s off for 10s, 30min between shocks).
9. # 16 Biological replicates were harvested. A two-color approach was used on a custom 8x60k Agilent microarray
10. # Dye color was counterbalanced across conditions.
12. # This script reads a design file, reads the agilent file, adjusts for background signal using normexp algorithm,
13. # normalizes for dye intensity within each array, filters out control genes, then averages by probeset (for transcripts measured
14. # multiple probes.
16. # Next, this script tests a number of research questions about regulation of expression in this experiment.
18. # Development of script--was based primarily on Limma user guide and this guide for two-color Agilent analysis:
19. # http://matticklab.com/index.php?title=Two_channel_analysis_of_Agilent_microarray_data_with_Limma
20. # with some additional help from the references cited below.
22. # Decisions:
23. # Median or mean signal? Zahurak et al. (2007) report the choice makes no practical difference in their analysis of Agilent
24. # two-color arrays. Therefore, default of median was selected.
25. # Background subtraction? Zahurak et al. (2007) report best results with no correction. However, in a more extensive comparison
26. # Ritchie et al. (2007) report best results with normexp+offset for background correction.
27. # Normexp offset? Offset needs to be empirically determined by inspection of MA plots, but analysis is not super-sensitive to
28. # choice of offset value (k). Although Ritchie et al. used an offset of 50 based on their previous work, inspection of
29. # MA plots of 15, 30, 50, and 100 show basically no from k=30-100. MAPlots for different
30. # k values from this data set are in the MAPlots subfolder. k value is currently set to 30 based on a prior microarray
31. # analysis, but will need to be set based on MAPlots from this 5 day data set, when available
32. # Normalize Within Arrays? Zahurak et al. (2007) found best results for Agilent two-color with loess normalization. This also
33. # seems to have become the standard within-array normalization strategy.
34. # How to assess practical significance? Most statistical testing is done against an LFC of 0, which leads to some results
35. # which are significant but not meaningful in terms of the magnitude of regulation. This can be fixed by a) adding an LFC
36. # threshold after significance testing. Alternatively, McCarthy & Smith, 2009 recommend conducting statistical testing
37. # against a more reasonable null hypothesis, thus testing for practical and statistical significance at the same time.
38. # as using eBayes with an LFC of 0 resulted in many statistically significant but small differences between groups,
39. # the approach of checking against a more rigorous null was used. Specifically, an LFC of 1.5 was selected (0.176).
40. # At the end of the script, comments show how to do the LFC=0 approach for the null.
43. #Setup
44. # will need to setup directory to wherever microarray data is eventually posted
45. # setwd("Z:/projects/2015 - LTS maintenance and decay - paired pleurals/LTS 1 week - microarray - paired pleurals")
46. library(limma)
48. #Read design and data
49. # will need to setup targets file with correct file names and colors
50. targets <- readTargets("targets.txt")
51. RG <- read.maimages(targets, source="agilent.median")
53. #Preprocessing
54. RG <- backgroundCorrect(RG, method="normexp", offset=30)
55. MA <- normalizeWithinArrays(RG, method="loess")
56. write.table(MA, file="./output/GeoTable.txt")
57. MA2 <- MA[MA$genes$ControlType==0,]
58. MA2.avg <- avereps(MA2, ID=MA2$genes$GeneName)
59. write.table(MA2.avg, file="./output/MA2Table_combined.txt")
61. #Make the MAPlots
62. for (i in 1:nrow(targets)) {
63. jpeg(file=paste("./output/", targets$FileName[i]," - MAplot.jpeg", sep=""))
64. plotMA(MA, i, status=MA$genes$ControlType)
65. dev.off()
66. }
68. #Research questions:
69. ## 1) How many of the 24h regulated transcripts are regulated at this time point?
70. ## Transcripts persistently regulated from 24h to 5 day could be related to retention. Transcripts regulated in this pattern will be identified by reading the list of 24h transcripts and then comparing trained/untrained expression at 5 days within just this set of transcripts. The test will be a standard empirical Bayes t-test. Significant differences will be identified by testing against a null of at least 10% change in eaither direction with correction for multiple comparisons (B-H, FDR = 0.05).
71. ## 1b) To what extent is regulation at 1 day similar to regulation at 5 days?
72. ## Will use the genas function. In this case we will do a combined analysis across the 1day and 5day microarray data and then identify the degree to which the pattern of regulation at 1 day is preserved at 5 days. This will be done with the p_union switch to look only at the connection amongst regulated transcripts.
74. ## 2) What transcripts are NEWLY regulated at this time point, 5 days after LTS training?
75. ## Transcripts regulated at 5 day but not 1 day could represent transcripts related to active forgetting processes. This question will be addressed by running the same analysis as above (comparing trained/untrained) across the whole microarray WITHOUT the the list of 1-day regulated. This will identify *newly* regulated transcripts. Again the approach will be a simple comparison of trained to untrained expression, usin the moderated empirical Bayes t-test comparing the trained and untrained side. Significant differences will be identified by testing against a null of at least 10% change in either direction with correction for mulltiple comparisons (B-H, FDR = 0.05)
76. ## 2a) How many missed NEWLY regulated transcripts?
77. ## Will use the propTrueNull function with the convest approach to assess how many false negatives might have occured in the microarray less the 1day transcripts.
79. ## 3) What transcripts might be related to the degree of forgetting?
80. # Will examine transcripts that significantly predict forgetting scores
81. ## Forgetting scores are calculated as log(T_1d/U_1d) - log(T_5d/U_5d)
82. ## The first term represents sharpness of memory at 1 day. The second term represents sharpness of memory at 5 days. The difference represents amount of memory lost, with higher numbers indicating more forgetting. This analyis will add a term for forgetting and examine across the whole set of 5-day microarray data.
84. ## 4) What happens to 5-day regulated transcripts at 7 days?
85. ## Not in this script, but read the 7 day data, restrict down to 5-day regulated transcripts, and test for LFC>1.1
86. ## In other words, repeat question 1 exactly but using 7 day data and the list of 5-day regulated genes
87. ## 4a) How much of the pattern of regulation at 5 days is preserved at 7 days?
88. ## Not in this script, but read 5 and 7 day data, use genas to examine correlation using p.union
89. ## In other words, repeat analysis for question 1a but comparing 5 and 7 day data.
91. ## Setup design matrix
92. design <- modelMatrix(targets, ref="Control")
93. design <- cbind(Dye=1, design)
96. #Research question 1 (and same pattern for 4)
97. # restrict down to just the transcripts that were up at 1h (or 7 days for question 4)
98. restricted = read.csv("refined_24h.txt")
99. MA2_restricted <- MA2.avg[MA2.avg$genes$GeneName %in% restricted$Transcripts, ]
100. fit_restricted <- lmFit(MA2_restricted, design)
101. regulated_restricted <- treat(fit_restricted, lfc=log2(1.1))
102. regulatedTable_restricted <- topTreat(regulated_restricted, coef="Trained", number=nrow(regulated_restricted))
103. write.table(regulatedTable_restricted, file="./output/LTS5day_retention_related_lfc1_10.txt")
106. #Research question 2 - Test for significant FCs > log (1.1) only amongst genes not regulated at 1day
107. MA2_unique <- MA2.avg[!(MA2.avg$genes$GeneName %in% restricted$Transcripts), ]
108. fit_unique <- lmFit(MA2_unique, design)
109. regulated_unique <- treat(fit_unique, lfc=log2(1.1))
110. regulatedTable_unique <- topTreat(regulated_unique, coef="Trained", number=nrow(regulated_unique))
111. write.table(regulatedTable_unique, file="./output/LTS5day_forgetting_related_lfc_1_10.txt")
113. #Research question 2a
114. propTrueNull(regulatedTable_unique$P.Value, method="convest")
117. #Research question 3
118. i <- targets$Cy5=="Control"
119. #THis is the array of forgetting scores for the 16 arrays. Be sure it matches order of arrays.
120. mscore <- c(0.443219007884632, 0.131376770816482, 0.639579778630896, 0.247227348439673, 0.622799920975808, 0.200682793496743, 0.937516492070804, 0.157722366671975, 0.620463296505017, 0.114347963363, 0.462221812522715, 0.163534611418199, 0.816323784510695, 0.281895301846058, 0.630446328710265, 0.380411562929875)
121. mscore[i] = -mscore[i]
122. design_cor <- cbind(design, mscore=mscore)
123. fit_cor <- lmFit(MA2.avg, design_cor)
124. regulated_cor <- treat(fit_cor, lfc=log2(1.1))
125. regulated_cor_condition <- topTreat(regulated_cor, coef="Trained", number=nrow(regulated_cor))
126. regulated_cor_retention <- topTreat(regulated_cor, coef="mscore")
127. write.table(regulated_cor_condition, file="./output/LTS5day_corr_cond_lfc1_10.txt")
128. write.table(regulated_cor_retention, file="./output/LTS5day_corr_retention_lfc1_10.txt")
131. #Research question 1a
132. ## Re-do script to read in not only 5 day but also 1 day data, each as own group, then use geneas function to compare regulation
133. #Read design and data
134. # will need to setup targets file with correct file names and colors
136. ## !! This feels dangerous. Better to split this into a different script--it would be easy to not properly rename
137. ## all the veriables for this analysis, and thus end up smooshing together this analysis with leftover portions
138. ## of the other analyses.
140. targets_comb <- readTargets("targets_compare_1d5d.txt")
141. RG_comb <- read.maimages(targets_comb, source="agilent.median")
143. #Preprocessing
144. RG_comb <- backgroundCorrect(RG_comb, method="normexp", offset=30)
145. MA_comb <- normalizeWithinArrays(RG_comb, method="loess")
146. MA2_comb <- MA_comb[MA_comb$genes$ControlType==0,]
147. MA2_comb.avg <- avereps(MA2_comb, ID=MA2$genes$GeneName)
149. #Test for significant FCs > log (1.1)
150. design_comb <- modelMatrix(targets_comb, ref="Control")
151. design_comb <- cbind(Dye=1, design_comb)
152. fit_comb <- lmFit(MA2_comb.avg, design_comb)
153. regulated_comb <- treat(fit_comb, lfc=log2(1.1))
155. genas(fit_comb, coef=c("h1", "h24"), subset="p.union", plot=TRUE, alpha=0.4)
156. genas(fit_comb, coef=c("d1", "d5"), subset="p.union", plot=TRUE, alpha=0.4)

lts 5day - microarray script - preregistered.R, no license · at the source