-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathThree_species_GITHUB.R
944 lines (786 loc) · 53.3 KB
/
Three_species_GITHUB.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
###############################################################################################
library(readr)
library(DESeq2)
library(HTSFilter)
library(data.table)
d = "F:/LabWork/Three-Species/"
setwd(d)
three_species_annot = fread("Trinotate_annotation/AMEX_ThreeSpecies_associations_wTrinotate_allSPECIES.txt")
amex.eData = read.table("Ax_RSEM.gene.counts.matrix", row.names = 1, header = TRUE, sep='\t')
amex.eData = amex.eData[which(rownames(amex.eData)%in%three_species_annot$AxREFv2),]
#amex.eData = amex.eData[,c(-3,-8)]
amex.pData = data.frame("samp.names" = c("AMEX_D0_1","AMEX_D0_2","AMEX_D0_3",
"AMEX_D1_1","AMEX_D1_2","AMEX_D1_3",
"AMEX_D2_1","AMEX_D2_2","AMEX_D2_3",
"AMEX_D3_1","AMEX_D3_2","AMEX_D3_3",
"AMEX_D4_1","AMEX_D4_2","AMEX_D4_3",
"AMEX_D5_1","AMEX_D5_2","AMEX_D5_3"),
"dpa" = as.factor(c(rep(0,3),rep(1,3),rep(2,3),rep(3,3),rep(4,3),rep(5,3))))
rownames(amex.pData) = colnames(amex.eData)
#amex.dds <- amex.dds[rowSums(fpm(amex.dds)>0)>=18]
register(bpstart(SnowParam(6)))
##### AMEX - generate DESeq set #######
amex.dds <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1,2,4,6,7,9,10,12,13,14,17,18)]), colData = amex.pData[c(1,2,4,6,7,9,10,12,13,14,17,18),], design = ~dpa)
amex_dds2 <- DESeq(amex.dds, test=c("LRT"), reduced =~1, parallel = TRUE)
#amex_filtered <- HTSFilter(amex_dds2, s.len=100, normalization = "DESeq", parallel = TRUE)$filteredData
#amex_LRT = results(amex_dds2, independentFiltering=TRUE)
#hist(amex_LRT$pvalue)
#summary(amex_LRT)
###
amex.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1,2,4,6)]), colData = amex.pData[c(1,2,4,6),], design = ~dpa)
amex.dds_1v0 <- amex.dds_1v0[rowSums(fpm(amex.dds_1v0)>0)>=4]
amex_dds2_D1vs0 <- DESeq(amex.dds_1v0, test=c("Wald"), parallel = TRUE)
amex_D1vD0 = results(amex_dds2_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(amex_D1vD0$pvalue)
summary(amex_D1vD0, alpha = 0.05)
amex.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(c(1,2,7,9))]), colData = amex.pData[c(c(1,2,7,9)),], design = ~dpa)
amex.dds_2v0 <- amex.dds_2v0[rowSums(fpm(amex.dds_2v0)>0)>=4]
amex_dds2_D2vs0 <- DESeq(amex.dds_2v0, test=c("Wald"), parallel = TRUE)
amex_D2vD0 = results(amex_dds2_D2vs0, contrast = c("dpa","2","0"), independentFiltering=TRUE)
hist(amex_D2vD0$pvalue)
summary(amex_D2vD0, alpha = 0.05)
amex.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1:2,10,12)]), colData = amex.pData[c(1:2,10,12),], design = ~dpa)
amex.dds_3v0 <- amex.dds_3v0[rowSums(fpm(amex.dds_3v0)>0)>=4]
amex_dds2_D3vs0 <- DESeq(amex.dds_3v0, test=c("Wald"), parallel = TRUE)
amex_D3vD0 = results(amex_dds2_D3vs0, contrast = c("dpa","3","0"), independentFiltering=TRUE)
hist(amex_D3vD0$pvalue)
summary(amex_D3vD0, alpha = 0.05)
amex.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1:2,13:14)]), colData = amex.pData[c(1:2,13:14),], design = ~dpa)
amex.dds_4v0 <- amex.dds_4v0[rowSums(fpm(amex.dds_4v0)>0)>=4]
amex_dds2_D4vs0 <- DESeq(amex.dds_4v0, test=c("Wald"), parallel = TRUE)
amex_D4vD0 = results(amex_dds2_D4vs0, contrast = c("dpa","4","0"), independentFiltering=TRUE)
hist(amex_D4vD0$pvalue)
summary(amex_D4vD0, alpha = 0.05)
amex.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1:2,17:18)]), colData = amex.pData[c(1:2,17:18),], design = ~dpa)
amex.dds_5v0 <- amex.dds_5v0[rowSums(fpm(amex.dds_5v0)>0)>=4]
amex_dds2_D5vs0 <- DESeq(amex.dds_5v0, test=c("Wald"), parallel = TRUE)
amex_D5vD0 = results(amex_dds2_D5vs0, contrast = c("dpa","5","0"), independentFiltering=TRUE)
hist(amex_D5vD0$pvalue)
summary(amex_D5vD0, alpha = 0.05)
###
amex_D1vD0_df = data.frame(amex_D1vD0)
amex_D1vD0_df$txid = rownames(amex_D1vD0_df)
amex_D2vD0_df = data.frame(amex_D2vD0)
amex_D2vD0_df$txid = rownames(amex_D2vD0_df)
amex_D3vD0_df = data.frame(amex_D3vD0)
amex_D3vD0_df$txid = rownames(amex_D3vD0_df)
amex_D4vD0_df = data.frame(amex_D4vD0)
amex_D4vD0_df$txid = rownames(amex_D4vD0_df)
amex_D5vD0_df = data.frame(amex_D5vD0)
amex_D5vD0_df$txid = rownames(amex_D5vD0_df)
###
amex_independent_tests = union(union(union(union(union(rownames(amex_D1vD0[which(amex_D1vD0$padj < 0.05),]),
rownames(amex_D2vD0[which(amex_D2vD0$padj < 0.05),])),
rownames(amex_D3vD0[which(amex_D3vD0$padj < 0.05),])),
rownames(amex_D4vD0[which(amex_D4vD0$padj < 0.05),])),
rownames(amex_D4vD0[which(amex_D4vD0$padj < 0.05),])),
rownames(amex_D5vD0[which(amex_D5vD0$padj < 0.05),]))
amex_DEG = data.frame(trans = amex_independent_tests)
rownames(amex_DEG) = amex_DEG[,1]
amex_DEG$D1v0_FC = amex_D1vD0_df[match(rownames(amex_DEG),rownames(amex_D1vD0_df)),c(2)]
amex_DEG$D2v0_FC = amex_D2vD0_df[match(rownames(amex_DEG),rownames(amex_D2vD0_df)),c(2)]
amex_DEG$D3v0_FC = amex_D3vD0_df[match(rownames(amex_DEG),rownames(amex_D3vD0_df)),c(2)]
amex_DEG$D4v0_FC = amex_D4vD0_df[match(rownames(amex_DEG),rownames(amex_D4vD0_df)),c(2)]
amex_DEG$D5v0_FC = amex_D5vD0_df[match(rownames(amex_DEG),rownames(amex_D5vD0_df)),c(2)]
amex_DEG$D1v0_FDR = amex_D1vD0_df[match(rownames(amex_DEG),rownames(amex_D1vD0_df)),c(6)]
amex_DEG$D2v0_FDR = amex_D2vD0_df[match(rownames(amex_DEG),rownames(amex_D2vD0_df)),c(6)]
amex_DEG$D3v0_FDR = amex_D3vD0_df[match(rownames(amex_DEG),rownames(amex_D3vD0_df)),c(6)]
amex_DEG$D4v0_FDR = amex_D4vD0_df[match(rownames(amex_DEG),rownames(amex_D4vD0_df)),c(6)]
amex_DEG$D5v0_FDR = amex_D5vD0_df[match(rownames(amex_DEG),rownames(amex_D5vD0_df)),c(6)]
write.table(amex_DEG, "DEG_list/Amex_DEG_fulllist.txt",row.names = TRUE, col.names = TRUE, sep='\t')
pearson_microarray<-cor(fpm(amex_dds2)[rownames(counts(amex_dds2, norm=TRUE))%in%amex_independent_tests,],method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
amex_fpm_norm = data.frame(fpm(amex_dds2)[rownames(counts(amex_dds2, norm=TRUE))%in%amex_independent_tests,])
amex_fpm_norm$D0_avg = rowMeans(amex_fpm_norm[,c(1:2)])
amex_fpm_norm$D1_avg = rowMeans(amex_fpm_norm[,c(3:4)])
amex_fpm_norm$D2_avg = rowMeans(amex_fpm_norm[,c(5:6)])
amex_fpm_norm$D3_avg = rowMeans(amex_fpm_norm[,c(7:8)])
amex_fpm_norm$D4_avg = rowMeans(amex_fpm_norm[,c(9:10)])
amex_fpm_norm$D5_avg = rowMeans(amex_fpm_norm[,c(11:12)])
library(pvclust)
y = pvclust(data=amex_fpm_norm[,c(13:18)],method.hclust="complete",method.dist="correlation",nboot=1000, parallel = TRUE)
tiff(filename = "figures/AMEX_TimepointCluster_Bootstrap.tiff", width = 6, height = 4, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
#plotMDS(fpm(amex_dds2)[rownames(counts(amex_dds2, norm=TRUE))%in%amex_independent_tests,])
nrow(amex_D1vD0[which(amex_D1vD0$padj < 0.05 & amex_D1vD0$log2FoldChange > 0),])
nrow(amex_D2vD0[which(amex_D2vD0$padj < 0.05& amex_D2vD0$log2FoldChange > 0),])
nrow(amex_D3vD0[which(amex_D3vD0$padj < 0.05& amex_D3vD0$log2FoldChange > 0),])
nrow(amex_D4vD0[which(amex_D4vD0$padj < 0.05& amex_D4vD0$log2FoldChange > 0),])
nrow(amex_D5vD0[which(amex_D5vD0$padj < 0.05& amex_D5vD0$log2FoldChange > 0),])
nrow(amex_D1vD0[which(amex_D1vD0$padj < 0.05& amex_D1vD0$log2FoldChange < 0),])
nrow(amex_D2vD0[which(amex_D2vD0$padj < 0.05& amex_D2vD0$log2FoldChange < 0),])
nrow(amex_D3vD0[which(amex_D3vD0$padj < 0.05& amex_D3vD0$log2FoldChange < 0),])
nrow(amex_D4vD0[which(amex_D4vD0$padj < 0.05& amex_D4vD0$log2FoldChange < 0),])
nrow(amex_D5vD0[which(amex_D5vD0$padj < 0.05& amex_D5vD0$log2FoldChange < 0),])
write.table(amex_D1vD0_df, "amex_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D2vD0_df, "amex_D2vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D3vD0_df, "amex_D3vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D4vD0_df, "amex_D4vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D5vD0_df, "amex_D5vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D1vD0_df[which(amex_D1vD0_df$padj<0.05),], "amex_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D2vD0_df[which(amex_D2vD0_df$padj<0.05),], "amex_D2vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D3vD0_df[which(amex_D3vD0_df$padj<0.05),], "amex_D3vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D4vD0_df[which(amex_D4vD0_df$padj<0.05),], "amex_D4vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D5vD0_df[which(amex_D5vD0_df$padj<0.05),], "amex_D5vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
###########
am1 = rownames(amex_D1vD0_df[which(amex_D1vD0_df$padj<0.05),])
am2 = unique(c(rownames(amex_D2vD0_df[which(amex_D2vD0_df$padj<0.05),]),am1))
am3 = unique(c(rownames(amex_D3vD0_df[which(amex_D3vD0_df$padj<0.05),]),am2))
am4 = unique(c(rownames(amex_D4vD0_df[which(amex_D4vD0_df$padj<0.05),]),am3))
am5 = unique(c(rownames(amex_D5vD0_df[which(amex_D5vD0_df$padj<0.05),]),am4))
c(length(am1), length(am2), length(am3), length(am4), length(am5))
#save.image(file = "ThreeSpecies.RData")
########################################
########## AAND ########################
########################################
library(readr)
library(DESeq2)
library(HTSFilter)
d = "F:/LabWork/Three-Species/"
setwd(d)
and.eData = read.table("And_RSEM.gene.counts.matrix", row.names = 1, header = TRUE, sep='\t')
and.eData = and.eData[which(rownames(and.eData)%in%three_species_annot$AndREFv2),]
and.pData = data.frame("samp.names" = c("And_D0_1","And_D0_2","And_D0_3",
"And_D1_1","And_D1_2","And_D1_3",
"And_D2_1","And_D2_2","And_D2_3",
"And_D3_1","And_D3_2","And_D3_3",
"And_D4_1","And_D4_2","And_D4_3",
"And_D5_1","And_D5_2","And_D5_3"),
"dpa" = as.factor(c(rep(0,3),rep(1,3),rep(2,3),rep(3,3),rep(4,3),rep(5,3))))
rownames(and.pData) = colnames(and.eData)
and.dds <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,4,6,7:9,10,12,13,15,16,17)]), colData = and.pData[c(1:3,4,6,7:9,10,12,13,15,16,17),], design = ~dpa)
#and.dds <- and.dds[rowSums(fpm(and.dds)>0)>=9]
register(bpstart(SnowParam(6)))
and_dds <- DESeq(and.dds, test=c("LRT"), reduced =~1, parallel = TRUE)
##### AAND - run comparisons #######
#and_LRT = results(and_dds, independentFiltering=TRUE)
#hist(and_LRT$pvalue)
#summary(and_LRT, alpha = 0.05)
###
and.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,4,6)]), colData = and.pData[c(1:3,4,6),], design = ~dpa)
and.dds_1v0 <- and.dds_1v0[rowSums(fpm(and.dds_1v0)>0)>=5]
and_dds_D1vs0 <- DESeq(and.dds_1v0, test=c("Wald"), parallel = TRUE)
and_D1vD0 = results(and_dds_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(and_D1vD0$pvalue)
summary(and_D1vD0, alpha = 0.05)
and.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,7:9)]), colData = and.pData[c(1:3,7:9),], design = ~dpa)
and.dds_2v0 <- and.dds_2v0[rowSums(fpm(and.dds_2v0)>0)>=6]
and_dds_D2vs0 <- DESeq(and.dds_2v0, test=c("Wald"), parallel = TRUE)
and_D2vD0 = results(and_dds_D2vs0, contrast = c("dpa","2","0"), independentFiltering=TRUE)
hist(and_D2vD0$pvalue)
summary(and_D2vD0, alpha = 0.05)
and.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,10,12)]), colData = and.pData[c(1:3,10,12),], design = ~dpa)
and.dds_3v0 <- and.dds_3v0[rowSums(fpm(and.dds_3v0)>0)>=5]
and_dds_D3vs0 <- DESeq(and.dds_3v0, test=c("Wald"), parallel = TRUE)
and_D3vD0 = results(and_dds_D3vs0, contrast = c("dpa","3","0"), independentFiltering=TRUE)
hist(and_D3vD0$pvalue)
summary(and_D3vD0, alpha = 0.05)
and.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,13,15)]), colData = and.pData[c(1:3,13,15),], design = ~dpa)
and.dds_4v0 <- and.dds_4v0[rowSums(fpm(and.dds_4v0)>0)>=5]
and_dds_D4vs0 <- DESeq(and.dds_4v0, test=c("Wald"), parallel = TRUE)
and_D4vD0 = results(and_dds_D4vs0, contrast = c("dpa","4","0"), independentFiltering=TRUE)
hist(and_D4vD0$pvalue)
summary(and_D4vD0, alpha = 0.05)
and.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,16,17)]), colData = and.pData[c(1:3,16,17),], design = ~dpa)
and.dds_5v0 <- and.dds_5v0[rowSums(fpm(and.dds_5v0)>0)>=5]
and_dds_D5vs0 <- DESeq(and.dds_5v0, test=c("Wald"), parallel = TRUE)
and_D5vD0 = results(and_dds_D5vs0, contrast = c("dpa","5","0"), independentFiltering=TRUE)
hist(and_D5vD0$pvalue)
summary(and_D5vD0, alpha = 0.05)
and_independent_tests = union(union(union(union(union(rownames(and_D1vD0[which(and_D1vD0$padj < 0.05),]),
rownames(and_D2vD0[which(and_D2vD0$padj < 0.05),])),
rownames(and_D3vD0[which(and_D3vD0$padj < 0.05),])),
rownames(and_D4vD0[which(and_D4vD0$padj < 0.05),])),
rownames(and_D4vD0[which(and_D4vD0$padj < 0.05),])),
rownames(and_D5vD0[which(and_D5vD0$padj < 0.05),]))
and_D1vD0_df = data.frame(and_D1vD0)
and_D1vD0_df$txid = rownames(and_D1vD0_df)
and_D2vD0_df = data.frame(and_D2vD0)
and_D2vD0_df$txid = rownames(and_D2vD0_df)
and_D3vD0_df = data.frame(and_D3vD0)
and_D3vD0_df$txid = rownames(and_D3vD0_df)
and_D4vD0_df = data.frame(and_D4vD0)
and_D4vD0_df$txid = rownames(and_D4vD0_df)
and_D5vD0_df = data.frame(and_D5vD0)
and_D5vD0_df$txid = rownames(and_D5vD0_df)
and_DEG = data.frame(trans = and_independent_tests)
rownames(and_DEG) = and_DEG[,1]
and_DEG$D1v0_FC = and_D1vD0_df[match(rownames(and_DEG),rownames(and_D1vD0_df)),c(2)]
and_DEG$D2v0_FC = and_D2vD0_df[match(rownames(and_DEG),rownames(and_D2vD0_df)),c(2)]
and_DEG$D3v0_FC = and_D3vD0_df[match(rownames(and_DEG),rownames(and_D3vD0_df)),c(2)]
and_DEG$D4v0_FC = and_D4vD0_df[match(rownames(and_DEG),rownames(and_D4vD0_df)),c(2)]
and_DEG$D5v0_FC = and_D5vD0_df[match(rownames(and_DEG),rownames(and_D5vD0_df)),c(2)]
and_DEG$D1v0_FDR = and_D1vD0_df[match(rownames(and_DEG),rownames(and_D1vD0_df)),c(6)]
and_DEG$D2v0_FDR = and_D2vD0_df[match(rownames(and_DEG),rownames(and_D2vD0_df)),c(6)]
and_DEG$D3v0_FDR = and_D3vD0_df[match(rownames(and_DEG),rownames(and_D3vD0_df)),c(6)]
and_DEG$D4v0_FDR = and_D4vD0_df[match(rownames(and_DEG),rownames(and_D4vD0_df)),c(6)]
and_DEG$D5v0_FDR = and_D5vD0_df[match(rownames(and_DEG),rownames(and_D5vD0_df)),c(6)]
nrow(and_D1vD0[which(and_D1vD0$padj < 0.05 & and_D1vD0$log2FoldChange > 0),])
nrow(and_D2vD0[which(and_D2vD0$padj < 0.05& and_D2vD0$log2FoldChange > 0),])
nrow(and_D3vD0[which(and_D3vD0$padj < 0.05& and_D3vD0$log2FoldChange > 0),])
nrow(and_D4vD0[which(and_D4vD0$padj < 0.05& and_D4vD0$log2FoldChange > 0),])
nrow(and_D5vD0[which(and_D5vD0$padj < 0.05& and_D5vD0$log2FoldChange > 0),])
nrow(and_D1vD0[which(and_D1vD0$padj < 0.05& and_D1vD0$log2FoldChange < 0),])
nrow(and_D2vD0[which(and_D2vD0$padj < 0.05& and_D2vD0$log2FoldChange < 0),])
nrow(and_D3vD0[which(and_D3vD0$padj < 0.05& and_D3vD0$log2FoldChange < 0),])
nrow(and_D4vD0[which(and_D4vD0$padj < 0.05& and_D4vD0$log2FoldChange < 0),])
nrow(and_D5vD0[which(and_D5vD0$padj < 0.05& and_D5vD0$log2FoldChange < 0),])
write.table(and_D1vD0_df, "and_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D2vD0_df, "and_D2vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D3vD0_df, "and_D3vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D4vD0_df, "and_D4vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D5vD0_df, "and_D5vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D1vD0_df[which(and_D1vD0_df$padj<0.05),], "and_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D2vD0_df[which(and_D2vD0_df$padj<0.05),], "and_D2vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D3vD0_df[which(and_D3vD0_df$padj<0.05),], "and_D3vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D4vD0_df[which(and_D4vD0_df$padj<0.05),], "and_D4vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D5vD0_df[which(and_D5vD0_df$padj<0.05),], "and_D5vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
ad1 = rownames(and_D1vD0_df[which(and_D1vD0_df$padj<0.05),])
ad2 = unique(c(rownames(and_D2vD0_df[which(and_D2vD0_df$padj<0.05),]),ad1))
ad3 = unique(c(rownames(and_D3vD0_df[which(and_D3vD0_df$padj<0.05),]),ad2))
ad4 = unique(c(rownames(and_D4vD0_df[which(and_D4vD0_df$padj<0.05),]),ad3))
ad5 = unique(c(rownames(and_D5vD0_df[which(and_D5vD0_df$padj<0.05),]),ad4))
c(length(ad1), length(ad2), length(ad3), length(ad4), length(ad5))
write.table(and_DEG, "DEG_list/and_DEG_fulllist.txt",row.names = TRUE, col.names = TRUE, sep='\t')
pearson_microarray<-cor(fpm(and_dds)[rownames(counts(and_dds, norm=TRUE))%in%and_independent_tests,],method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
and_fpm_norm = data.frame(fpm(and_dds)[rownames(counts(and_dds, norm=TRUE))%in%and_independent_tests,])
and_fpm_norm$D0_avg = rowMeans(and_fpm_norm[,c(1:3)])
and_fpm_norm$D1_avg = rowMeans(and_fpm_norm[,c(4:5)])
and_fpm_norm$D2_avg = rowMeans(and_fpm_norm[,c(6:8)])
and_fpm_norm$D3_avg = rowMeans(and_fpm_norm[,c(9:10)])
and_fpm_norm$D4_avg = rowMeans(and_fpm_norm[,c(11:12)])
and_fpm_norm$D5_avg = rowMeans(and_fpm_norm[,c(13:14)])
library(pvclust)
y = pvclust(data=and_fpm_norm[,c(15:20)],method.hclust="complete",method.dist="correlation",nboot=1000, parallel = TRUE)
tiff(filename = "figures/AAND_TimepointCluster_Bootstrap.tiff", width = 6, height = 4, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
########################################
########## AMAC ########################
########################################
#load("ThreeSpecies.RData")
mac.eData = read.table("Mac_RSEM.gene.counts.matrix", row.names = 1, header = TRUE, sep='\t')
mac.eData = mac.eData[which(rownames(mac.eData)%in%three_species_annot$MacREFv2),]
#mac.eData = mac.eData[,c(-6)]
mac.pData = data.frame("samp.names" = c("Mac_D0_1","Mac_D0_2","Mac_D0_3",
"Mac_D1_1","Mac_D1_2","Mac_D1_3",
"Mac_D2_1","Mac_D2_2","Mac_D2_3",
"Mac_D3_1","Mac_D3_2","Mac_D3_3",
"Mac_D4_1","Mac_D4_2","Mac_D4_3",
"Mac_D5_1","Mac_D5_2","Mac_D5_3"),
"dpa" = as.factor(c(rep(0,3),rep(1,3),rep(2,3),rep(3,3),rep(4,3),rep(5,3))))
rownames(mac.pData) = colnames(mac.eData)
mac.dds <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,4,5,7:8,10:11,13:14,16:17)]), colData = mac.pData[c(1:3,4,5,7:8,10:11,13:14,16:17),], design = ~dpa)
#mac.dds <- mac.dds[rowSums(fpm(mac.dds)>0)>=9]
register(bpstart(SnowParam(6)))
mac_dds <- DESeq(mac.dds, test=c("LRT"), reduced =~1, parallel = TRUE )
##### AMAC - generate DESeq set #######
mac_LRT = results(mac_dds, independentFiltering=TRUE)
hist(mac_LRT$pvalue)
summary(mac_LRT)
### comparison
mac.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:5)]), colData = mac.pData[c(1:5),], design = ~dpa)
mac.dds_1v0 <- mac.dds_1v0[rowSums(fpm(mac.dds_1v0)>0)>=5]
mac_dds_D1vs0 <- DESeq(mac.dds_1v0, test=c("Wald"), parallel = TRUE)
mac_D1vD0 = results(mac_dds_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(mac_D1vD0$pvalue)
summary(mac_D1vD0, alpha = 0.05)
mac.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,7:8)]), colData = mac.pData[c(1:3,7:8),], design = ~dpa)
mac.dds_2v0 <- mac.dds_2v0[rowSums(fpm(mac.dds_2v0)>0)>=5]
mac_dds_D2vs0 <- DESeq(mac.dds_2v0, test=c("Wald"), parallel = TRUE)
mac_D2vD0 = results(mac_dds_D2vs0, contrast = c("dpa","2","0"), independentFiltering=TRUE)
hist(mac_D2vD0$pvalue)
summary(mac_D2vD0, alpha = 0.05)
mac.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,10:11)]), colData = mac.pData[c(1:3,10:11),], design = ~dpa)
mac.dds_3v0 <- mac.dds_3v0[rowSums(fpm(mac.dds_3v0)>0)>=5]
mac_dds_D3vs0 <- DESeq(mac.dds_3v0, test=c("Wald"), parallel = TRUE)
mac_D3vD0 = results(mac_dds_D3vs0, contrast = c("dpa","3","0"), independentFiltering=TRUE)
hist(mac_D3vD0$pvalue)
summary(mac_D3vD0, alpha = 0.05)
mac.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,13:14)]), colData = mac.pData[c(1:3,13:14),], design = ~dpa)
mac.dds_4v0 <- mac.dds_4v0[rowSums(fpm(mac.dds_4v0)>0)>=5]
mac_dds_D4vs0 <- DESeq(mac.dds_4v0, test=c("Wald"), parallel = TRUE)
mac_D4vD0 = results(mac_dds_D4vs0, contrast = c("dpa","4","0"), independentFiltering=TRUE)
hist(mac_D4vD0$pvalue)
summary(mac_D4vD0, alpha = 0.05)
mac.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,16:17)]), colData = mac.pData[c(1:3,16:17),], design = ~dpa)
mac.dds_5v0 <- mac.dds_5v0[rowSums(fpm(mac.dds_5v0)>0)>=5]
mac_dds_D5vs0 <- DESeq(mac.dds_5v0, test=c("Wald"), parallel = TRUE)
mac_D5vD0 = results(mac_dds_D5vs0, contrast = c("dpa","5","0"), independentFiltering=TRUE)
hist(mac_D5vD0$pvalue)
summary(mac_D5vD0, alpha = 0.05)
mac_independent_tests = union(union(union(union(union(rownames(mac_D1vD0[which(mac_D1vD0$padj < 0.05),]),
rownames(mac_D2vD0[which(mac_D2vD0$padj < 0.05),])),
rownames(mac_D3vD0[which(mac_D3vD0$padj < 0.05),])),
rownames(mac_D4vD0[which(mac_D4vD0$padj < 0.05),])),
rownames(mac_D4vD0[which(mac_D4vD0$padj < 0.05),])),
rownames(mac_D5vD0[which(mac_D5vD0$padj < 0.05),]))
mac_D1vD0_df = data.frame(mac_D1vD0)
mac_D1vD0_df$txid = rownames(mac_D1vD0_df)
mac_D2vD0_df = data.frame(mac_D2vD0)
mac_D2vD0_df$txid = rownames(mac_D2vD0_df)
mac_D3vD0_df = data.frame(mac_D3vD0)
mac_D3vD0_df$txid = rownames(mac_D3vD0_df)
mac_D4vD0_df = data.frame(mac_D4vD0)
mac_D4vD0_df$txid = rownames(mac_D4vD0_df)
mac_D5vD0_df = data.frame(mac_D5vD0)
mac_D5vD0_df$txid = rownames(mac_D5vD0_df)
mac_DEG = data.frame(trans = mac_independent_tests)
rownames(mac_DEG) = mac_DEG[,1]
mac_DEG$D1v0_FC = mac_D1vD0_df[match(rownames(mac_DEG),rownames(mac_D1vD0_df)),c(2)]
mac_DEG$D2v0_FC = mac_D2vD0_df[match(rownames(mac_DEG),rownames(mac_D2vD0_df)),c(2)]
mac_DEG$D3v0_FC = mac_D3vD0_df[match(rownames(mac_DEG),rownames(mac_D3vD0_df)),c(2)]
mac_DEG$D4v0_FC = mac_D4vD0_df[match(rownames(mac_DEG),rownames(mac_D4vD0_df)),c(2)]
mac_DEG$D5v0_FC = mac_D5vD0_df[match(rownames(mac_DEG),rownames(mac_D5vD0_df)),c(2)]
mac_DEG$D1v0_FDR = mac_D1vD0_df[match(rownames(mac_DEG),rownames(mac_D1vD0_df)),c(6)]
mac_DEG$D2v0_FDR = mac_D2vD0_df[match(rownames(mac_DEG),rownames(mac_D2vD0_df)),c(6)]
mac_DEG$D3v0_FDR = mac_D3vD0_df[match(rownames(mac_DEG),rownames(mac_D3vD0_df)),c(6)]
mac_DEG$D4v0_FDR = mac_D4vD0_df[match(rownames(mac_DEG),rownames(mac_D4vD0_df)),c(6)]
mac_DEG$D5v0_FDR = mac_D5vD0_df[match(rownames(mac_DEG),rownames(mac_D5vD0_df)),c(6)]
write.table(mac_DEG, "DEG_list/mac_DEG_fulllist.txt",row.names = TRUE, col.names = TRUE, sep='\t')
pearson_microarray<-cor(fpm(mac_dds)[rownames(counts(mac_dds, norm=TRUE))%in%mac_independent_tests,],method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
mac_fpm_norm = data.frame(fpm(mac_dds)[rownames(counts(mac_dds, norm=TRUE))%in%mac_independent_tests,])
mac_fpm_norm$D0_avg = rowMeans(mac_fpm_norm[,c(1:3)])
mac_fpm_norm$D1_avg = rowMeans(mac_fpm_norm[,c(4:5)])
mac_fpm_norm$D2_avg = rowMeans(mac_fpm_norm[,c(6:7)])
mac_fpm_norm$D3_avg = rowMeans(mac_fpm_norm[,c(8:9)])
mac_fpm_norm$D4_avg = rowMeans(mac_fpm_norm[,c(10:11)])
mac_fpm_norm$D5_avg = rowMeans(mac_fpm_norm[,c(12:13)])
library(pvclust)
y = pvclust(data=mac_fpm_norm[,c(14:19)],method.hclust="complete",method.dist="correlation",nboot=1000, parallel = TRUE)
tiff(filename = "figures/AMAC_TimepointCluster_Bootstrap.tiff", width = 6, height = 4, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
nrow(mac_D1vD0[which(mac_D1vD0$padj < 0.05 & mac_D1vD0$log2FoldChange > 0),])
nrow(mac_D2vD0[which(mac_D2vD0$padj < 0.05& mac_D2vD0$log2FoldChange > 0),])
nrow(mac_D3vD0[which(mac_D3vD0$padj < 0.05& mac_D3vD0$log2FoldChange > 0),])
nrow(mac_D4vD0[which(mac_D4vD0$padj < 0.05& mac_D4vD0$log2FoldChange > 0),])
nrow(mac_D5vD0[which(mac_D5vD0$padj < 0.05& mac_D5vD0$log2FoldChange > 0),])
nrow(mac_D1vD0[which(mac_D1vD0$padj < 0.05 & mac_D1vD0$log2FoldChange < 0),])
nrow(mac_D2vD0[which(mac_D2vD0$padj < 0.05& mac_D2vD0$log2FoldChange < 0),])
nrow(mac_D3vD0[which(mac_D3vD0$padj < 0.05& mac_D3vD0$log2FoldChange < 0),])
nrow(mac_D4vD0[which(mac_D4vD0$padj < 0.05& mac_D4vD0$log2FoldChange < 0),])
nrow(mac_D5vD0[which(mac_D5vD0$padj < 0.05& mac_D5vD0$log2FoldChange < 0),])
write.table(mac_D1vD0_df, "mac_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D2vD0_df, "mac_D2vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D3vD0_df, "mac_D3vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D4vD0_df, "mac_D4vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D5vD0_df, "mac_D5vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D1vD0_df[which(mac_D1vD0_df$padj<0.05),], "mac_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D2vD0_df[which(mac_D2vD0_df$padj<0.05),], "mac_D2vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D3vD0_df[which(mac_D3vD0_df$padj<0.05),], "mac_D3vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D4vD0_df[which(mac_D4vD0_df$padj<0.05),], "mac_D4vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D5vD0_df[which(mac_D5vD0_df$padj<0.05),], "mac_D5vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
mac1 = rownames(mac_D1vD0_df[which(mac_D1vD0_df$padj<0.05),])
mac2 = unique(c(rownames(mac_D2vD0_df[which(mac_D2vD0_df$padj<0.05),]),mac1))
mac3 = unique(c(rownames(mac_D3vD0_df[which(mac_D3vD0_df$padj<0.05),]),mac2))
mac4 = unique(c(rownames(mac_D4vD0_df[which(mac_D4vD0_df$padj<0.05),]),mac3))
mac5 = unique(c(rownames(mac_D5vD0_df[which(mac_D5vD0_df$padj<0.05),]),mac4))
c(length(mac1), length(mac2), length(mac3), length(mac4), length(mac5))
#################################################################################################
amex.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(1,2,4,6)]), colData = amex.pData[c(1,2,4,6),], design = ~dpa)
amex.dds_1v0 <- amex.dds_1v0[rowSums(fpm(amex.dds_1v0)>0)>=4]
amex_dds2_D1vs0 <- DESeq(amex.dds_1v0, test=c("Wald"), parallel = TRUE)
amex_D1vD0 = results(amex_dds2_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(amex_D1vD0$pvalue)
summary(amex_D1vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(amex_dds2_D1vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="ward.D")
plot(pearson_micro_clust)
amex.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(4,6,7,9)]), colData = amex.pData[c(4,6,7,9),], design = ~dpa)
amex.dds_2v0 <- amex.dds_2v0[rowSums(fpm(amex.dds_2v0)>0)>=4]
amex_dds2_D2vs0 <- DESeq(amex.dds_2v0, test=c("Wald"), parallel = TRUE)
amex_D2vD0 = results(amex_dds2_D2vs0, contrast = c("dpa","2","1"), independentFiltering=TRUE)
hist(amex_D2vD0$pvalue)
summary(amex_D2vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(amex_dds2_D2vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
amex.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(7,9,10,12)]), colData = amex.pData[c(7,9,10,12),], design = ~dpa)
amex.dds_3v0 <- amex.dds_3v0[rowSums(fpm(amex.dds_3v0)>0)>=4]
amex_dds2_D3vs0 <- DESeq(amex.dds_3v0, test=c("Wald"), parallel = TRUE)
amex_D3vD0 = results(amex_dds2_D3vs0, contrast = c("dpa","3","2"), independentFiltering=TRUE)
hist(amex_D3vD0$pvalue)
summary(amex_D3vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(amex_dds2_D3vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="ward.D")
plot(pearson_micro_clust)
amex.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(10,12,13:14)]), colData = amex.pData[c(10,12,13:14),], design = ~dpa)
amex.dds_4v0 <- amex.dds_4v0[rowSums(fpm(amex.dds_4v0)>0)>=4]
amex_dds2_D4vs0 <- DESeq(amex.dds_4v0, test=c("Wald"), parallel = TRUE)
amex_D4vD0 = results(amex_dds2_D4vs0, contrast = c("dpa","4","3"), independentFiltering=TRUE)
hist(amex_D4vD0$pvalue)
summary(amex_D4vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(amex_dds2_D4vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
amex.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(amex.eData[,c(13:14,17:18)]), colData = amex.pData[c(13:14,17:18),], design = ~dpa)
amex.dds_5v0 <- amex.dds_5v0[rowSums(fpm(amex.dds_5v0)>0)>=4]
amex_dds2_D5vs0 <- DESeq(amex.dds_5v0, test=c("Wald"), parallel = TRUE)
amex_D5vD0 = results(amex_dds2_D5vs0, contrast = c("dpa","5","4"), independentFiltering=TRUE)
hist(amex_D5vD0$pvalue)
summary(amex_D5vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(amex_dds2_D5vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
###
amex_D1vD0_df = data.frame(amex_D1vD0)
amex_D1vD0_df$txid = rownames(amex_D1vD0_df)
amex_D1vD0_df$geneID = amex_D1vD0_df[which(rownames(amex_D1vD0_df)%in%three_species_annot$AxREFv2),]
amex_D2vD0_df = data.frame(amex_D2vD0)
amex_D2vD0_df$txid = rownames(amex_D2vD0_df)
amex_D2vD0_df$geneID = amex_D2vD0_df[which(rownames(amex_D2vD0_df)%in%three_species_annot$AxREFv2),]
amex_D3vD0_df = data.frame(amex_D3vD0)
amex_D3vD0_df$txid = rownames(amex_D3vD0_df)
amex_D3vD0_df$geneID = amex_D3vD0_df[which(rownames(amex_D3vD0_df)%in%three_species_annot$AxREFv2),]
amex_D4vD0_df = data.frame(amex_D4vD0)
amex_D4vD0_df$txid = rownames(amex_D4vD0_df)
amex_D4vD0_df$geneID = amex_D4vD0_df[which(rownames(amex_D4vD0_df)%in%three_species_annot$AxREFv2),]
amex_D5vD0_df = data.frame(amex_D5vD0)
amex_D5vD0_df$txid = rownames(amex_D5vD0_df)
amex_D5vD0_df$geneID = amex_D5vD0_df[which(rownames(amex_D5vD0_df)%in%three_species_annot$AxREFv2),]
###
amex_independent_tests = union(union(union(union(union(rownames(amex_D1vD0[which(amex_D1vD0$padj < 0.05),]),
rownames(amex_D2vD0[which(amex_D2vD0$padj < 0.05),])),
rownames(amex_D3vD0[which(amex_D3vD0$padj < 0.05),])),
rownames(amex_D4vD0[which(amex_D4vD0$padj < 0.05),])),
rownames(amex_D4vD0[which(amex_D4vD0$padj < 0.05),])),
rownames(amex_D5vD0[which(amex_D5vD0$padj < 0.05),]))
amex_DEG = data.frame(trans = amex_independent_tests)
rownames(amex_DEG) = amex_DEG[,1]
amex_DEG$D1v0_FC = amex_D1vD0_df[match(rownames(amex_DEG),rownames(amex_D1vD0_df)),c(2)]
amex_DEG$D2v0_FC = amex_D2vD0_df[match(rownames(amex_DEG),rownames(amex_D2vD0_df)),c(2)]
amex_DEG$D3v0_FC = amex_D3vD0_df[match(rownames(amex_DEG),rownames(amex_D3vD0_df)),c(2)]
amex_DEG$D4v0_FC = amex_D4vD0_df[match(rownames(amex_DEG),rownames(amex_D4vD0_df)),c(2)]
amex_DEG$D5v0_FC = amex_D5vD0_df[match(rownames(amex_DEG),rownames(amex_D5vD0_df)),c(2)]
amex_DEG$D1v0_FDR = amex_D1vD0_df[match(rownames(amex_DEG),rownames(amex_D1vD0_df)),c(6)]
amex_DEG$D2v0_FDR = amex_D2vD0_df[match(rownames(amex_DEG),rownames(amex_D2vD0_df)),c(6)]
amex_DEG$D3v0_FDR = amex_D3vD0_df[match(rownames(amex_DEG),rownames(amex_D3vD0_df)),c(6)]
amex_DEG$D4v0_FDR = amex_D4vD0_df[match(rownames(amex_DEG),rownames(amex_D4vD0_df)),c(6)]
amex_DEG$D5v0_FDR = amex_D5vD0_df[match(rownames(amex_DEG),rownames(amex_D5vD0_df)),c(6)]
#amex_DEG$geneID = amex_DEG[which(rownames(amex_DEG)%in%three_species_annot$AxREFv2),]
write.table(amex_DEG, "DEG_list/Amex_DEG_fulllist_adjacent.txt",row.names = TRUE, col.names = TRUE, sep='\t')
pearson_microarray<-cor(fpm(amex_dds2)[rownames(counts(amex_dds2, norm=TRUE))%in%amex_independent_tests,],method="spearman")
#pearson_microarray<-cor(fpm(amex_dds2),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
amex_fpm_norm = data.frame(fpm(amex_dds2)[rownames(counts(amex_dds2, norm=TRUE))%in%amex_independent_tests,])
amex_fpm_norm$D0_avg = rowMeans(amex_fpm_norm[,c(1:2)])
amex_fpm_norm$D1_avg = rowMeans(amex_fpm_norm[,c(3:4)])
amex_fpm_norm$D2_avg = rowMeans(amex_fpm_norm[,c(5:6)])
amex_fpm_norm$D3_avg = rowMeans(amex_fpm_norm[,c(7:8)])
amex_fpm_norm$D4_avg = rowMeans(amex_fpm_norm[,c(9:10)])
amex_fpm_norm$D5_avg = rowMeans(amex_fpm_norm[,c(11:12)])
library(pvclust)
y = pvclust(data=amex_fpm_norm[,c(13:18)],method.hclust="complete",method.dist="correlation",nboot=1000, parallel = TRUE)
tiff(filename = "figures/AMEX_TimepointCluster_adjacent_Bootstrap.tiff", width = 6, height = 4, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
nrow(amex_D1vD0[which(amex_D1vD0$padj < 0.05 & amex_D1vD0$log2FoldChange > 0),])
nrow(amex_D2vD0[which(amex_D2vD0$padj < 0.05 & amex_D2vD0$log2FoldChange > 0),])
nrow(amex_D3vD0[which(amex_D3vD0$padj < 0.05 & amex_D3vD0$log2FoldChange > 0),])
nrow(amex_D4vD0[which(amex_D4vD0$padj < 0.05 & amex_D4vD0$log2FoldChange > 0),])
nrow(amex_D5vD0[which(amex_D5vD0$padj < 0.05 & amex_D5vD0$log2FoldChange > 0),])
nrow(amex_D1vD0[which(amex_D1vD0$padj < 0.05 & amex_D1vD0$log2FoldChange < 0),])
nrow(amex_D2vD0[which(amex_D2vD0$padj < 0.05 & amex_D2vD0$log2FoldChange < 0),])
nrow(amex_D3vD0[which(amex_D3vD0$padj < 0.05 & amex_D3vD0$log2FoldChange < 0),])
nrow(amex_D4vD0[which(amex_D4vD0$padj < 0.05 & amex_D4vD0$log2FoldChange < 0),])
nrow(amex_D5vD0[which(amex_D5vD0$padj < 0.05 & amex_D5vD0$log2FoldChange < 0),])
write.table(amex_D1vD0_df, "amex_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D2vD0_df, "amex_D2vD1_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D3vD0_df, "amex_D3vD2_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D4vD0_df, "amex_D4vD3_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D5vD0_df, "amex_D5vD4_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D1vD0_df[which(amex_D1vD0_df$padj<0.05),], "amex_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D2vD0_df[which(amex_D2vD0_df$padj<0.05),], "amex_D2vD1_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D3vD0_df[which(amex_D3vD0_df$padj<0.05),], "amex_D3vD2_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D4vD0_df[which(amex_D4vD0_df$padj<0.05),], "amex_D4vD3_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(amex_D5vD0_df[which(amex_D5vD0_df$padj<0.05),], "amex_D5vD4_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
am1 = rownames(amex_D1vD0_df[which(amex_D1vD0_df$padj<0.05),])
am2 = unique(c(rownames(amex_D2vD0_df[which(amex_D2vD0_df$padj<0.05),]),am1))
am3 = unique(c(rownames(amex_D3vD0_df[which(amex_D3vD0_df$padj<0.05),]),am2))
am4 = unique(c(rownames(amex_D4vD0_df[which(amex_D4vD0_df$padj<0.05),]),am3))
am5 = unique(c(rownames(amex_D5vD0_df[which(amex_D5vD0_df$padj<0.05),]),am4))
##### AAND - run comparisons #######
and.dds <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,4,6,7:9,10,12,13,15,16,17)]), colData = and.pData[c(1:3,4,6,7:9,10,12,13,15,16,17),], design = ~dpa)
#and.dds <- and.dds[rowSums(fpm(and.dds)>0)>=15]
and_dds <- DESeq(and.dds, test=c("LRT"), reduced =~1, parallel = TRUE)
#and_LRT = results(and_dds, independentFiltering=TRUE)
#hist(and_LRT$pvalue)
#summary(and_LRT, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
###
and.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(1:3,4,6)]), colData = and.pData[c(1:3,4,6),], design = ~dpa)
and.dds_1v0 <- and.dds_1v0[rowSums(fpm(and.dds_1v0)>0)>=5]
and_dds_D1vs0 <- DESeq(and.dds_1v0, test=c("Wald"), parallel = TRUE)
and_D1vD0 = results(and_dds_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(and_D1vD0$pvalue)
summary(and_D1vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds_D1vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
and.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(4,6,7:9)]), colData = and.pData[c(4,6,7:9),], design = ~dpa)
and.dds_2v0 <- and.dds_2v0[rowSums(fpm(and.dds_2v0)>0)>=5]
and_dds_D2vs0 <- DESeq(and.dds_2v0, test=c("Wald"), parallel = TRUE)
and_D2vD0 = results(and_dds_D2vs0, contrast = c("dpa","2","1"), independentFiltering=TRUE)
hist(and_D2vD0$pvalue)
summary(and_D2vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds_D2vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
and.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(7:9,10,12)]), colData = and.pData[c(7:9,10,12),], design = ~dpa)
and.dds_3v0 <- and.dds_3v0[rowSums(fpm(and.dds_3v0)>0)>=5]
and_dds_D3vs0 <- DESeq(and.dds_3v0, test=c("Wald"), parallel = TRUE)
and_D3vD0 = results(and_dds_D3vs0, contrast = c("dpa","3","2"), independentFiltering=TRUE)
hist(and_D3vD0$pvalue)
summary(and_D3vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds_D3vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="average")
plot(pearson_micro_clust)
and.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(10,12,13,15)]), colData = and.pData[c(10,12,13,15),], design = ~dpa)
and.dds_4v0 <- and.dds_4v0[rowSums(fpm(and.dds_4v0)>0)>=4]
and_dds_D4vs0 <- DESeq(and.dds_4v0, test=c("Wald"), parallel = TRUE)
and_D4vD0 = results(and_dds_D4vs0, contrast = c("dpa","4","3"), independentFiltering=TRUE)
hist(and_D4vD0$pvalue)
summary(and_D4vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds_D4vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
and.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(and.eData[,c(13,15,16:17)]), colData = and.pData[c(13,15,16:17),], design = ~dpa)
and.dds_5v0 <- and.dds_5v0[rowSums(fpm(and.dds_5v0)>0)>=4]
and_dds_D5vs0 <- DESeq(and.dds_5v0, test=c("Wald"), parallel = TRUE)
and_D5vD0 = results(and_dds_D5vs0, contrast = c("dpa","5","4"), independentFiltering=TRUE)
hist(and_D5vD0$pvalue)
summary(and_D5vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(and_dds_D5vs0,norm=TRUE)+1),method="pearson")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
####
and_independent_tests = union(union(union(union(union(rownames(and_D1vD0[which(and_D1vD0$padj < 0.05),]),
rownames(and_D2vD0[which(and_D2vD0$padj < 0.05),])),
rownames(and_D3vD0[which(and_D3vD0$padj < 0.05),])),
rownames(and_D4vD0[which(and_D4vD0$padj < 0.05),])),
rownames(and_D4vD0[which(and_D4vD0$padj < 0.05),])),
rownames(and_D5vD0[which(and_D5vD0$padj < 0.05),]))
and_D1vD0_df = data.frame(and_D1vD0)
and_D1vD0_df$txid = rownames(and_D1vD0_df)
and_D2vD0_df = data.frame(and_D2vD0)
and_D2vD0_df$txid = rownames(and_D2vD0_df)
and_D3vD0_df = data.frame(and_D3vD0)
and_D3vD0_df$txid = rownames(and_D3vD0_df)
and_D4vD0_df = data.frame(and_D4vD0)
and_D4vD0_df$txid = rownames(and_D4vD0_df)
and_D5vD0_df = data.frame(and_D5vD0)
and_D5vD0_df$txid = rownames(and_D5vD0_df)
and_DEG = data.frame(trans = and_independent_tests)
rownames(and_DEG) = and_DEG[,1]
and_DEG$D1v0_FC = and_D1vD0_df[match(rownames(and_DEG),rownames(and_D1vD0_df)),c(2)]
and_DEG$D2v0_FC = and_D2vD0_df[match(rownames(and_DEG),rownames(and_D2vD0_df)),c(2)]
and_DEG$D3v0_FC = and_D3vD0_df[match(rownames(and_DEG),rownames(and_D3vD0_df)),c(2)]
and_DEG$D4v0_FC = and_D4vD0_df[match(rownames(and_DEG),rownames(and_D4vD0_df)),c(2)]
and_DEG$D5v0_FC = and_D5vD0_df[match(rownames(and_DEG),rownames(and_D5vD0_df)),c(2)]
and_DEG$D1v0_FDR = and_D1vD0_df[match(rownames(and_DEG),rownames(and_D1vD0_df)),c(6)]
and_DEG$D2v0_FDR = and_D2vD0_df[match(rownames(and_DEG),rownames(and_D2vD0_df)),c(6)]
and_DEG$D3v0_FDR = and_D3vD0_df[match(rownames(and_DEG),rownames(and_D3vD0_df)),c(6)]
and_DEG$D4v0_FDR = and_D4vD0_df[match(rownames(and_DEG),rownames(and_D4vD0_df)),c(6)]
and_DEG$D5v0_FDR = and_D5vD0_df[match(rownames(and_DEG),rownames(and_D5vD0_df)),c(6)]
write.table(and_DEG, "DEG_list/and_DEG_fulllist_adjacent.txt",row.names = TRUE, col.names = TRUE, sep='\t')
library(pvclust)
y = pvclust(data=amby_DEGs_Rv0[,c(6:20)],method.hclust="complete",method.dist="correlation",nboot=10000, parallel = TRUE)
tiff(filename = "figures/Ambystoma_DEG_Rv0Comparisons.tiff", width = 6, height = 8, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
nrow(and_D1vD0[which(and_D1vD0$padj < 0.05 & and_D1vD0$log2FoldChange > 0),])
nrow(and_D2vD0[which(and_D2vD0$padj < 0.05 & and_D2vD0$log2FoldChange > 0),])
nrow(and_D3vD0[which(and_D3vD0$padj < 0.05& and_D3vD0$log2FoldChange > 0),])
nrow(and_D4vD0[which(and_D4vD0$padj < 0.05& and_D4vD0$log2FoldChange > 0),])
nrow(and_D5vD0[which(and_D5vD0$padj < 0.05& and_D5vD0$log2FoldChange > 0),])
nrow(and_D1vD0[which(and_D1vD0$padj < 0.05& and_D1vD0$log2FoldChange < 0),])
nrow(and_D2vD0[which(and_D2vD0$padj < 0.05& and_D2vD0$log2FoldChange < 0),])
nrow(and_D3vD0[which(and_D3vD0$padj < 0.05& and_D3vD0$log2FoldChange < 0),])
nrow(and_D4vD0[which(and_D4vD0$padj < 0.05& and_D4vD0$log2FoldChange < 0),])
nrow(and_D5vD0[which(and_D5vD0$padj < 0.05& and_D5vD0$log2FoldChange < 0),])
write.table(and_D1vD0_df, "and_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D2vD0_df, "and_D2vD1_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D3vD0_df, "and_D3vD2_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D4vD0_df, "and_D4vD3_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D5vD0_df, "and_D5vD4_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D1vD0_df[which(and_D1vD0_df$padj<0.05),], "and_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D2vD0_df[which(and_D2vD0_df$padj<0.05),], "and_D2vD1_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D3vD0_df[which(and_D3vD0_df$padj<0.05),], "and_D3vD2_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D4vD0_df[which(and_D4vD0_df$padj<0.05),], "and_D4vD3_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(and_D5vD0_df[which(and_D5vD0_df$padj<0.05),], "and_D5vD4_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
and1 = rownames(and_D1vD0_df[which(and_D1vD0_df$padj<0.05),])
and2 = unique(c(rownames(and_D2vD0_df[which(and_D2vD0_df$padj<0.05),]),and1))
and3 = unique(c(rownames(and_D3vD0_df[which(and_D3vD0_df$padj<0.05),]),and2))
and4 = unique(c(rownames(and_D4vD0_df[which(and_D4vD0_df$padj<0.05),]),and3))
and5 = unique(c(rownames(and_D5vD0_df[which(and_D5vD0_df$padj<0.05),]),and4))
##### AMAC - generate DESeq set #######
#mac_LRT = results(mac_dds, independentFiltering=TRUE)
#hist(mac_LRT$pvalue)
#summary(mac_LRT)
### comparison
mac.dds_1v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(1:3,4:5)]), colData = mac.pData[c(1:3,4:5),], design = ~dpa)
mac.dds_1v0 <- mac.dds_1v0[rowSums(fpm(mac.dds_1v0)>0)>=5]
mac_dds_D1vs0 <- DESeq(mac.dds_1v0, test=c("Wald"), parallel = TRUE)
mac_D1vD0 = results(mac_dds_D1vs0, contrast = c("dpa","1","0"), independentFiltering=TRUE)
hist(mac_D1vD0$pvalue)
summary(mac_D1vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(mac_dds_D1vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
mac.dds_2v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(4:5,7:8)]), colData = mac.pData[c(4:5,7:8),], design = ~dpa)
mac.dds_2v0 <- mac.dds_2v0[rowSums(fpm(mac.dds_2v0)>0)>=4]
mac_dds_D2vs0 <- DESeq(mac.dds_2v0, test=c("Wald"), parallel = TRUE)
mac_D2vD0 = results(mac_dds_D2vs0, contrast = c("dpa","2","1"), independentFiltering=TRUE)
hist(mac_D2vD0$pvalue)
summary(mac_D2vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(mac_dds_D2vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
mac.dds_3v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(7:8,10:11)]), colData = mac.pData[c(7:8,10:11),], design = ~dpa)
mac.dds_3v0 <- mac.dds_3v0[rowSums(fpm(mac.dds_3v0)>0)>=4]
mac_dds_D3vs0 <- DESeq(mac.dds_3v0, test=c("Wald"), parallel = TRUE)
mac_D3vD0 = results(mac_dds_D3vs0, contrast = c("dpa","3","2"), independentFiltering=TRUE)
hist(mac_D3vD0$pvalue)
summary(mac_D3vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(mac_dds_D3vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
mac.dds_4v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(10:11,13:14)]), colData = mac.pData[c(10:11,13:14),], design = ~dpa)
mac.dds_4v0 <- mac.dds_4v0[rowSums(fpm(mac.dds_4v0)>0)>=4]
mac_dds_D4vs0 <- DESeq(mac.dds_4v0, test=c("Wald"), parallel = TRUE)
mac_D4vD0 = results(mac_dds_D4vs0, contrast = c("dpa","4","3"), independentFiltering=TRUE)
hist(mac_D4vD0$pvalue)
summary(mac_D4vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(mac_dds_D4vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="ward.D")
plot(pearson_micro_clust)
mac.dds_5v0 <- DESeqDataSetFromMatrix(countData = round(mac.eData[,c(13:14,16:17)]), colData = mac.pData[c(13:14,16:17),], design = ~dpa)
mac.dds_5v0 <- mac.dds_5v0[rowSums(fpm(mac.dds_5v0)>0)>=4]
mac_dds_D5vs0 <- DESeq(mac.dds_5v0, test=c("Wald"), parallel = TRUE)
mac_D5vD0 = results(mac_dds_D5vs0, contrast = c("dpa","5","4"), independentFiltering=TRUE)
hist(mac_D5vD0$pvalue)
summary(mac_D5vD0, alpha = 0.05)
pearson_microarray<-cor(log2(counts(mac_dds_D5vs0,norm=TRUE)+1),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="ward.D")
plot(pearson_micro_clust)
mac_independent_tests = union(union(union(union(union(rownames(mac_D1vD0[which(mac_D1vD0$padj < 0.05),]),
rownames(mac_D2vD0[which(mac_D2vD0$padj < 0.05),])),
rownames(mac_D3vD0[which(mac_D3vD0$padj < 0.05),])),
rownames(mac_D4vD0[which(mac_D4vD0$padj < 0.05),])),
rownames(mac_D4vD0[which(mac_D4vD0$padj < 0.05),])),
rownames(mac_D5vD0[which(mac_D5vD0$padj < 0.05),]))
mac_D1vD0_df = data.frame(mac_D1vD0)
mac_D1vD0_df$txid = rownames(mac_D1vD0_df)
mac_D2vD0_df = data.frame(mac_D2vD0)
mac_D2vD0_df$txid = rownames(mac_D2vD0_df)
mac_D3vD0_df = data.frame(mac_D3vD0)
mac_D3vD0_df$txid = rownames(mac_D3vD0_df)
mac_D4vD0_df = data.frame(mac_D4vD0)
mac_D4vD0_df$txid = rownames(mac_D4vD0_df)
mac_D5vD0_df = data.frame(mac_D5vD0)
mac_D5vD0_df$txid = rownames(mac_D5vD0_df)
mac_DEG = data.frame(trans = mac_independent_tests)
rownames(mac_DEG) = mac_DEG[,1]
mac_DEG$D1v0_FC = mac_D1vD0_df[match(rownames(mac_DEG),rownames(mac_D1vD0_df)),c(2)]
mac_DEG$D2v0_FC = mac_D2vD0_df[match(rownames(mac_DEG),rownames(mac_D2vD0_df)),c(2)]
mac_DEG$D3v0_FC = mac_D3vD0_df[match(rownames(mac_DEG),rownames(mac_D3vD0_df)),c(2)]
mac_DEG$D4v0_FC = mac_D4vD0_df[match(rownames(mac_DEG),rownames(mac_D4vD0_df)),c(2)]
mac_DEG$D5v0_FC = mac_D5vD0_df[match(rownames(mac_DEG),rownames(mac_D5vD0_df)),c(2)]
mac_DEG$D1v0_FDR = mac_D1vD0_df[match(rownames(mac_DEG),rownames(mac_D1vD0_df)),c(6)]
mac_DEG$D2v0_FDR = mac_D2vD0_df[match(rownames(mac_DEG),rownames(mac_D2vD0_df)),c(6)]
mac_DEG$D3v0_FDR = mac_D3vD0_df[match(rownames(mac_DEG),rownames(mac_D3vD0_df)),c(6)]
mac_DEG$D4v0_FDR = mac_D4vD0_df[match(rownames(mac_DEG),rownames(mac_D4vD0_df)),c(6)]
mac_DEG$D5v0_FDR = mac_D5vD0_df[match(rownames(mac_DEG),rownames(mac_D5vD0_df)),c(6)]
write.table(mac_DEG, "DEG_list/mac_DEG_fulllist_adjacent.txt",row.names = TRUE, col.names = TRUE, sep='\t')
pearson_microarray<-cor(fpm(mac_dds)[rownames(counts(mac_dds, norm=TRUE))%in%mac_independent_tests,],method="spearman")
#pearson_microarray<-cor(fpm(amex_dds2),method="spearman")
diff_pearson<-1-pearson_microarray
dist_pearson_micro<-as.dist(diff_pearson)
pearson_micro_clust<-hclust(dist_pearson_micro,method="complete")
plot(pearson_micro_clust)
library(pvclust)
y = pvclust(data=mac_fpm_norm[,c(14:19)],method.hclust="complete",method.dist="correlation",nboot=1000, parallel = TRUE)
tiff(filename = "figures/AMAC_TimepointCluster_adjacent_Bootstrap.tiff", width = 6, height = 4, units = "in", res = 300)
plot(y)
pvrect(y, alpha=0.94)
dev.off()
nrow(mac_D1vD0[which(mac_D1vD0$padj < 0.05 & mac_D1vD0$log2FoldChange > 0),])
nrow(mac_D2vD0[which(mac_D2vD0$padj < 0.05 & mac_D2vD0$log2FoldChange > 0),])
nrow(mac_D3vD0[which(mac_D3vD0$padj < 0.05& mac_D3vD0$log2FoldChange > 0),])
nrow(mac_D4vD0[which(mac_D4vD0$padj < 0.05& mac_D4vD0$log2FoldChange > 0),])
nrow(mac_D5vD0[which(mac_D5vD0$padj < 0.05& mac_D5vD0$log2FoldChange > 0),])
nrow(mac_D1vD0[which(mac_D1vD0$padj < 0.05& mac_D1vD0$log2FoldChange < 0),])
nrow(mac_D2vD0[which(mac_D2vD0$padj < 0.05& mac_D2vD0$log2FoldChange < 0),])
nrow(mac_D3vD0[which(mac_D3vD0$padj < 0.05& mac_D3vD0$log2FoldChange < 0),])
nrow(mac_D4vD0[which(mac_D4vD0$padj < 0.05& mac_D4vD0$log2FoldChange < 0),])
nrow(mac_D5vD0[which(mac_D5vD0$padj < 0.05& mac_D5vD0$log2FoldChange < 0),])
write.table(mac_D1vD0_df, "mac_D1vD0_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D2vD0_df, "mac_D2vD1_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D3vD0_df, "mac_D3vD2_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D4vD0_df, "mac_D4vD3_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D5vD0_df, "mac_D5vD4_all.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D1vD0_df[which(mac_D1vD0_df$padj<0.05),], "mac_D1vD0_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D2vD0_df[which(mac_D2vD0_df$padj<0.05),], "mac_D2vD1_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D3vD0_df[which(mac_D3vD0_df$padj<0.05),], "mac_D3vD2_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D4vD0_df[which(mac_D4vD0_df$padj<0.05),], "mac_D4vD3_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
write.table(mac_D5vD0_df[which(mac_D5vD0_df$padj<0.05),], "mac_D5vD4_DEG.txt",row.names = TRUE, col.names = TRUE, sep='\t')
mac1 = rownames(mac_D1vD0_df[which(mac_D1vD0_df$padj<0.05),])
mac2 = unique(c(rownames(mac_D2vD0_df[which(mac_D2vD0_df$padj<0.05),]),mac1))
mac3 = unique(c(rownames(mac_D3vD0_df[which(mac_D3vD0_df$padj<0.05),]),mac2))
mac4 = unique(c(rownames(mac_D4vD0_df[which(mac_D4vD0_df$padj<0.05),]),mac3))
mac5 = unique(c(rownames(mac_D5vD0_df[which(mac_D5vD0_df$padj<0.05),]),mac4))