#/////////////////////////////////////////////////////////////////////////////
# LIBRARIES AND CUSTOM FUNCTIONS
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
library(edgeR)
## Loading required package: limma
library(car)
library(gplots)
## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009
##
## Attaching package: 'gplots'
##
## The following object is masked from 'package:stats':
##
## lowess
library(reshape)
library(stringr)
library(limma)
library(ggplot2)
require(biomaRt)
## Loading required package: biomaRt
library(stringr)
library(dplyr)
##
## Attaching package: 'dplyr'
##
## The following objects are masked from 'package:stats':
##
## filter, lag
##
## The following object is masked from '.env':
##
## n
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggdendro)
# Create folder for results files
setwd("~/Projects/2014-04-08_GBP_pig/Analysis/")
#/////////////////////////////////////////////////////////////////////////////
# Read in data ------------------------------------------------------------
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
## Read in the data making the row names the first column
counttable_all <- read.table("../RawData/featurecounts.merge", header=T, row.names=1)
head(counttable_all)
## Alignments.405_duodenum.405_duodenum.featureCounts
## ENSSSCG00000000001 77
## ENSSSCG00000000002 209
## ENSSSCG00000000003 3550
## ENSSSCG00000000004 6
## ENSSSCG00000000005 165
## ENSSSCG00000000006 207
## Alignments.405_ileum.405_ileum.featureCounts
## ENSSSCG00000000001 56
## ENSSSCG00000000002 411
## ENSSSCG00000000003 7000
## ENSSSCG00000000004 25
## ENSSSCG00000000005 379
## ENSSSCG00000000006 352
## Alignments.416_ileum.416_ileum.featureCounts
## ENSSSCG00000000001 81
## ENSSSCG00000000002 896
## ENSSSCG00000000003 3831
## ENSSSCG00000000004 33
## ENSSSCG00000000005 714
## ENSSSCG00000000006 82
## Alignments.41_duodenum.41_duodenum.featureCounts
## ENSSSCG00000000001 86
## ENSSSCG00000000002 213
## ENSSSCG00000000003 3478
## ENSSSCG00000000004 11
## ENSSSCG00000000005 160
## ENSSSCG00000000006 191
## Alignments.576_duodenum.576_duodenum.featureCounts
## ENSSSCG00000000001 24
## ENSSSCG00000000002 97
## ENSSSCG00000000003 1777
## ENSSSCG00000000004 1
## ENSSSCG00000000005 391
## ENSSSCG00000000006 17
## Alignments.576_ileum.576_ileum.featureCounts
## ENSSSCG00000000001 29
## ENSSSCG00000000002 376
## ENSSSCG00000000003 2881
## ENSSSCG00000000004 10
## ENSSSCG00000000005 372
## ENSSSCG00000000006 70
## Alignments.692_duodenum.692_duodenum.featureCounts
## ENSSSCG00000000001 44
## ENSSSCG00000000002 572
## ENSSSCG00000000003 2719
## ENSSSCG00000000004 6
## ENSSSCG00000000005 354
## ENSSSCG00000000006 45
## Alignments.692_ileum.692_ileum.featureCounts
## ENSSSCG00000000001 62
## ENSSSCG00000000002 600
## ENSSSCG00000000003 3426
## ENSSSCG00000000004 9
## ENSSSCG00000000005 442
## ENSSSCG00000000006 54
## Alignments.785_duodenum.785_duodenum.featureCounts
## ENSSSCG00000000001 86
## ENSSSCG00000000002 415
## ENSSSCG00000000003 3548
## ENSSSCG00000000004 27
## ENSSSCG00000000005 308
## ENSSSCG00000000006 497
## Alignments.785_ileum.785_ileum.featureCounts
## ENSSSCG00000000001 48
## ENSSSCG00000000002 408
## ENSSSCG00000000003 6608
## ENSSSCG00000000004 16
## ENSSSCG00000000005 508
## ENSSSCG00000000006 47
## Alignments.798_duodenum.798_duodenum.featureCounts
## ENSSSCG00000000001 104
## ENSSSCG00000000002 321
## ENSSSCG00000000003 1987
## ENSSSCG00000000004 14
## ENSSSCG00000000005 168
## ENSSSCG00000000006 342
## Alignments.798_ileum.798_ileum.featureCounts
## ENSSSCG00000000001 68
## ENSSSCG00000000002 305
## ENSSSCG00000000003 1700
## ENSSSCG00000000004 16
## ENSSSCG00000000005 165
## ENSSSCG00000000006 319
## Alignments.821_duodenum.821_duodenum.featureCounts
## ENSSSCG00000000001 42
## ENSSSCG00000000002 287
## ENSSSCG00000000003 2422
## ENSSSCG00000000004 2
## ENSSSCG00000000005 279
## ENSSSCG00000000006 72
## Alignments.821_ileum.821_ileum.featureCounts
## ENSSSCG00000000001 7
## ENSSSCG00000000002 565
## ENSSSCG00000000003 1204
## ENSSSCG00000000004 2
## ENSSSCG00000000005 1011
## ENSSSCG00000000006 2
## Alignments.825_duodenum.825_duodenum.featureCounts
## ENSSSCG00000000001 105
## ENSSSCG00000000002 353
## ENSSSCG00000000003 2521
## ENSSSCG00000000004 17
## ENSSSCG00000000005 217
## ENSSSCG00000000006 217
## Alignments.834_duodenum.834_duodenum.featureCounts
## ENSSSCG00000000001 99
## ENSSSCG00000000002 429
## ENSSSCG00000000003 1994
## ENSSSCG00000000004 13
## ENSSSCG00000000005 219
## ENSSSCG00000000006 393
## Alignments.834_ileum.834_ileum.featureCounts
## ENSSSCG00000000001 99
## ENSSSCG00000000002 431
## ENSSSCG00000000003 2383
## ENSSSCG00000000004 26
## ENSSSCG00000000005 273
## ENSSSCG00000000006 449
## Alignments.948_duodenum.948_duodenum.featureCounts
## ENSSSCG00000000001 77
## ENSSSCG00000000002 367
## ENSSSCG00000000003 1765
## ENSSSCG00000000004 10
## ENSSSCG00000000005 181
## ENSSSCG00000000006 392
## Alignments.948_ileum.948_ileum.featureCounts
## ENSSSCG00000000001 33
## ENSSSCG00000000002 571
## ENSSSCG00000000003 2598
## ENSSSCG00000000004 24
## ENSSSCG00000000005 290
## ENSSSCG00000000006 358
colnames(counttable_all)<-sapply(strsplit(colnames(counttable_all), "\\."), `[`, 3)
FPKM <- read.table("../RawData/cufflinks.merg", header=T, row.names=1)
colnames(FPKM)<-sapply(strsplit(colnames(FPKM), "\\."), `[`, 2)
#/////////////////////////////////////////////////////////////////////////////
# Bring in metadata
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
#Sham-grisar är 785, 948, 834, 821, 825
#GBP-grisar: 405, 416, 576, 692
Treatment <- c(rep("GBP",8),rep("Sham",11))
Tissue<-sapply(strsplit(colnames(counttable_all), "_"), `[`, 2)
SampleID<-sapply(strsplit(colnames(counttable_all), "_"), `[`, 1)
Meta<-data.frame(colnames(counttable_all),Tissue,SampleID,Treatment)
Meta
## colnames.counttable_all. Tissue SampleID Treatment
## 1 405_duodenum duodenum 405 GBP
## 2 405_ileum ileum 405 GBP
## 3 416_ileum ileum 416 GBP
## 4 41_duodenum duodenum 41 GBP
## 5 576_duodenum duodenum 576 GBP
## 6 576_ileum ileum 576 GBP
## 7 692_duodenum duodenum 692 GBP
## 8 692_ileum ileum 692 GBP
## 9 785_duodenum duodenum 785 Sham
## 10 785_ileum ileum 785 Sham
## 11 798_duodenum duodenum 798 Sham
## 12 798_ileum ileum 798 Sham
## 13 821_duodenum duodenum 821 Sham
## 14 821_ileum ileum 821 Sham
## 15 825_duodenum duodenum 825 Sham
## 16 834_duodenum duodenum 834 Sham
## 17 834_ileum ileum 834 Sham
## 18 948_duodenum duodenum 948 Sham
## 19 948_ileum ileum 948 Sham
#/////////////////////////////////////////////////////////////////////////////
# Prepare mart for gene annotation
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
mart <- useMart("ensembl")
mart <- useMart(biomart = "ensembl", dataset = "sscrofa_gene_ensembl")
annotation <- getBM(attributes=c('hgnc_symbol','ensembl_gene_id','description'),filters ='ensembl_gene_id',
values = rownames(counttable_all), mart = mart)
rownames(annotation) <- annotation$ensembl_gene_id
head(annotation)
## hgnc_symbol ensembl_gene_id description
## ENSSSCG00000005588 ENSSSCG00000005588
## ENSSSCG00000005673 ENSSSCG00000005673
## ENSSSCG00000005789 ENSSSCG00000005789
## ENSSSCG00000005793 ENSSSCG00000005793
## ENSSSCG00000005800 ENSSSCG00000005800
## ENSSSCG00000005808 ENSSSCG00000005808
#/////////////////////////////////////////////////////////////////////////////
# QC
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
# Normalize for QC
y <- DGEList(counts=counttable_all, genes=rownames(counttable_all))
isexpr <- rowSums(cpm(y)>1) >= 3
y <- y[isexpr,]
y <- calcNormFactors(y)
v <- voom(y,plot=F)
#plot(colSums(counttable_all),v$targets$lib.size)
par(mar=c(8,4,5,4))
barplot(colSums(counttable_all)/10e5,las=2,ylab="Million counts",col=factor(paste(Meta$Tissue,Meta$Treatment)),legend=F)
#/////////////////////////////////////////////////////////////////////////////
# Overall data structure/clustering
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
# TREE
d <- dist(t(v$E),method="euclidean") # euclidean distances between the rows
fit <- hclust(d, method="ward")
labelColors = c("#CDB380", "#036564", "#EB6841", "#EDC951")
# cut dendrogram in 4 clusters
clusMember = cutree(fit, 4)
clusMember <-as.integer(Meta$Tissue[fit$order])
names(clusMember) <- as.character(Meta[fit$order,1])
# function to get color labels
colLab <- function(n) {
if (is.leaf(n)) {
a <- attributes(n)
labCol <- labelColors[clusMember[which(names(clusMember) == a$label)]]
attr(n, "nodePar") <- c(a$nodePar, lab.col = labCol)
}
n
}
hcd = as.dendrogram(fit)
# using dendrapply
clusDendro = dendrapply(hcd, colLab)
# make plot
plot(clusDendro, type = "triangle")
heatmap.2(cor(v$E,method="spearman"),col=redblue(256),scale="none",key=TRUE, keysize=1.5,
ColSideColors=ifelse(Meta$Treatment=="Sham","red","black"),
RowSideColors=ifelse(Meta$Tissue=="ileum","grey","steelblue2"),
density.info="none", trace="none", cexRow=0.7,cexCol=1.2)
#par(mfrow=c(2,2))
d <- dist(t(v$E),method="euclidean") # euclidean distances between the rows
fit <- cmdscale(d,eig=TRUE, k=3) # k is the number of dim
fit # view results
## $points
## [,1] [,2] [,3]
## 405_duodenum -33.625 -17.536 34.794
## 405_ileum -19.868 2.135 22.200
## 416_ileum 43.374 32.599 -5.063
## 41_duodenum -30.336 -17.730 36.991
## 576_duodenum 82.095 -51.653 75.546
## 576_ileum 50.028 43.791 12.941
## 692_duodenum 55.337 53.461 11.744
## 692_ileum 54.888 54.425 13.415
## 785_duodenum -62.014 -26.524 -6.687
## 785_ileum 51.263 37.731 19.309
## 798_duodenum -71.189 -16.937 -20.229
## 798_ileum -67.757 -15.648 -21.197
## 821_duodenum -8.169 -28.377 20.339
## 821_ileum 233.507 -59.278 -67.830
## 825_duodenum -42.887 1.959 -15.465
## 834_duodenum -67.239 -15.745 -20.371
## 834_ileum -67.008 -15.696 -20.697
## 948_duodenum -68.653 -15.777 -21.875
## 948_ileum -31.746 54.800 -47.866
##
## $eig
## [1] 1.069e+05 2.282e+04 1.943e+04 1.126e+04 5.520e+03 4.944e+03
## [7] 3.609e+03 3.337e+03 2.695e+03 1.899e+03 1.811e+03 1.566e+03
## [13] 1.172e+03 5.913e+02 4.498e+02 3.915e+02 3.253e+02 2.838e+02
## [19] -1.131e-11
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 0.7892 0.7892
x <- fit$points[,1]
y <- fit$points[,2]
z <- fit$points[,3]
Comp <- data.frame(PC1=x,PC2=y,PC3=z)
PropoVarExplained <- fit$eig/sum(fit$eig)*100
qplot(PC1,PC2,size=3,data=Comp,colour=ifelse(Meta$Tissue=="ileum","grey","steelblue2"),shape=Meta$Treatment)+
xlab(paste("PC1 (",prettyNum(PropoVarExplained[1],digits=3)," %)",sep=""))+
ylab(paste("PC2 (",prettyNum(PropoVarExplained[2],digits=3)," %)",sep=""))
qplot(PC2,PC3,size=3,data=Comp,colour=ifelse(Meta$Tissue=="ileum","grey","steelblue2"),shape=Meta$Treatment)+
xlab(paste("PC2 (",prettyNum(PropoVarExplained[2],digits=3)," %)",sep=""))+
ylab(paste("PC3 (",prettyNum(PropoVarExplained[3],digits=3)," %)",sep=""))
qplot(PC1,PC3,size=3,data=Comp,colour=ifelse(Meta$Tissue=="ileum","grey","steelblue2"),shape=Meta$Treatment)+
xlab(paste("PC1 (",prettyNum(PropoVarExplained[1],digits=3)," %)",sep=""))+
ylab(paste("PC3 (",prettyNum(PropoVarExplained[2],digits=3)," %)",sep=""))
#/////////////////////////////////////////////////////////////////////////////
# Look at specific genes
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
geneList <- c("ENSSSCG00000012019","ENSSSCG00000009131","ENSSSCG00000021255")
par(mfrow=c(2,2))
for(genes in geneList) {
barplot(2^v$E[genes,],col=ifelse(Meta$Tissue=="ileum","grey","steelblue2"),las=2,ylab=paste("log2 CPM\n",genes))
}
#/////////////////////////////////////////////////////////////////////////////
# Inference
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
#/////////////////////////////////////////////////////////////////////////////
# Duodenum versus ileum
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
design <- model.matrix(~Tissue)
colnames(design)
## [1] "(Intercept)" "Tissueileum"
y <- DGEList(counts=counttable_all, genes=rownames(counttable_all))
isexpr <- rowSums(cpm(y)>1) >= 3
y <- y[isexpr,]
y <- calcNormFactors(y)
v <- voom(y,design,plot=F)
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
topTable_tissue <-topTable(fit2,coef="Tissueileum",p.value=1,n=200000,sort="p")
table(significant=topTable_tissue$adj.P.Val < .25,Upregulated=topTable_tissue$logFC > 0)
## Upregulated
## significant FALSE TRUE
## FALSE 7411 7071
# Annotate and write to file
is<-intersect(rownames(topTable_tissue),annotation$ensembl_gene_id)
topTable_tissue_annotated <- cbind(topTable_tissue[is,-1],annotation[is,])
head(topTable_tissue_annotated,n=20)
## logFC AveExpr t P.Value adj.P.Val B
## ENSSSCG00000004572 -1.4703 2.2846 -5.014 6.487e-05 0.4591 0.81757
## ENSSSCG00000013767 -1.4564 1.4717 -4.518 2.055e-04 0.4591 -0.16989
## ENSSSCG00000005249 -0.8455 3.8699 -4.478 2.257e-04 0.4591 0.39004
## ENSSSCG00000000991 -2.5193 0.6224 -4.328 3.208e-04 0.4591 -0.76019
## ENSSSCG00000013266 -1.0300 3.4909 -4.212 4.219e-04 0.4591 -0.16513
## ENSSSCG00000017219 -1.0597 6.2826 -4.184 4.497e-04 0.4591 0.01818
## ENSSSCG00000017251 -0.8144 5.4472 -4.151 4.867e-04 0.4591 -0.06078
## ENSSSCG00000029057 -1.1951 7.5081 -4.073 5.848e-04 0.4591 -0.22769
## ENSSSCG00000017368 -3.0773 -2.8453 -4.052 6.139e-04 0.4591 -1.96933
## ENSSSCG00000010639 -1.7680 4.0367 -4.033 6.422e-04 0.4591 -0.42216
## ENSSSCG00000008205 -1.4234 8.4891 -4.011 6.754e-04 0.4591 -0.40022
## ENSSSCG00000030927 -1.4579 9.4183 -3.992 7.064e-04 0.4591 -0.50850
## ENSSSCG00000028822 -1.5603 0.5602 -3.982 7.236e-04 0.4591 -1.19215
## ENSSSCG00000026605 -1.8690 10.4105 -3.949 7.824e-04 0.4591 -0.68307
## ENSSSCG00000030938 -1.6822 4.5090 -3.927 8.237e-04 0.4591 -0.56629
## ENSSSCG00000007983 -2.5405 1.5782 -3.923 8.300e-04 0.4591 -1.06960
## ENSSSCG00000013263 -1.3325 6.9623 -3.902 8.727e-04 0.4591 -0.54918
## ENSSSCG00000006560 -1.2261 5.3581 -3.876 9.284e-04 0.4591 -0.60953
## ENSSSCG00000006717 -1.7574 5.8278 -3.866 9.507e-04 0.4591 -0.62419
## ENSSSCG00000006243 -1.8549 4.5157 -3.865 9.518e-04 0.4591 -0.68852
## hgnc_symbol ensembl_gene_id
## ENSSSCG00000004572 ENSSSCG00000004572
## ENSSSCG00000013767 PALM3 ENSSSCG00000013767
## ENSSSCG00000005249 FAM189A2 ENSSSCG00000005249
## ENSSSCG00000000991 FOXQ1 ENSSSCG00000000991
## ENSSSCG00000013266 GYLTL1B ENSSSCG00000013266
## ENSSSCG00000017219 HID1 ENSSSCG00000017219
## ENSSSCG00000017251 ENSSSCG00000017251
## ENSSSCG00000029057 ENSSSCG00000029057
## ENSSSCG00000017368 ENSSSCG00000017368
## ENSSSCG00000010639 ENSSSCG00000010639
## ENSSSCG00000008205 ENSSSCG00000008205
## ENSSSCG00000030927 ENSSSCG00000030927
## ENSSSCG00000028822 ENSSSCG00000028822
## ENSSSCG00000026605 ENSSSCG00000026605
## ENSSSCG00000030938 ENSSSCG00000030938
## ENSSSCG00000007983 PDIA2 ENSSSCG00000007983
## ENSSSCG00000013263 CREB3L1 ENSSSCG00000013263
## ENSSSCG00000006560 CREB3L4 ENSSSCG00000006560
## ENSSSCG00000006717 ENSSSCG00000006717
## ENSSSCG00000006243 PENK ENSSSCG00000006243
## description
## ENSSSCG00000004572 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1S077]
## ENSSSCG00000013767 paralemmin 3 [Source:HGNC Symbol;Acc:33274]
## ENSSSCG00000005249 family with sequence similarity 189, member A2 [Source:HGNC Symbol;Acc:24820]
## ENSSSCG00000000991 forkhead box Q1 [Source:HGNC Symbol;Acc:20951]
## ENSSSCG00000013266 glycosyltransferase-like 1B [Source:HGNC Symbol;Acc:16522]
## ENSSSCG00000017219 HID1 domain containing [Source:HGNC Symbol;Acc:15736]
## ENSSSCG00000017251 Sus scrofa SRY (sex determining region Y)-box 9 (SOX9), mRNA. [Source:RefSeq mRNA;Acc:NM_213843]
## ENSSSCG00000029057 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LBK0]
## ENSSSCG00000017368 Pancreatic prohormone precursor Pancreatic hormone Pancreatic icosapeptide [Source:UniProtKB/Swiss-Prot;Acc:P01300]
## ENSSSCG00000010639 hyaluronan-binding protein 2 [Source:RefSeq peptide;Acc:NP_001230619]
## ENSSSCG00000008205 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1STC3]
## ENSSSCG00000030927
## ENSSSCG00000028822 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3L8A2]
## ENSSSCG00000026605 Sus scrofa bactericidal/permeability-increasing protein (BPI), mRNA. [Source:RefSeq mRNA;Acc:NM_001159307]
## ENSSSCG00000030938 Novel protein similar to vertebrate phosphoglycerate dehydrogenase (PHGDH) [Source:UniProtKB/TrEMBL;Acc:A5GFY6]
## ENSSSCG00000007983 protein disulfide isomerase family A, member 2 [Source:HGNC Symbol;Acc:14180]
## ENSSSCG00000013263 cAMP responsive element binding protein 3-like 1 [Source:HGNC Symbol;Acc:18856]
## ENSSSCG00000006560 cAMP responsive element binding protein 3-like 4 [Source:HGNC Symbol;Acc:18854]
## ENSSSCG00000006717 Sus scrofa phosphoglycerate dehydrogenase (PHGDH), mRNA. [Source:RefSeq mRNA;Acc:NM_001123162]
## ENSSSCG00000006243 proenkephalin [Source:HGNC Symbol;Acc:8831]
write.csv(file="topTable_tissue.csv",topTable_tissue_annotated)
par(mfrow=c(3,3))
for (gene in c(topTable_tissue$genes[1:9])) {
boxplot(v$E[gene,]~paste(Meta$Treatment,Meta$Tissue),ylab=gene,col=c("grey","steelblue2"),
main=paste(topTable_tissue_annotated[gene,]$hgnc_symbol, "\nAdj. p-val:",prettyNum(topTable_tissue[gene,]$adj.P.Val)),cex=2)
}
heatmap.2(v$E[topTable_tissue[topTable_tissue$adj.P.Val < .10,]$genes,Meta$Tissue=="tissue"],col=redblue(256),
scale="row",key=TRUE, keysize=1.5,
ColSideColors=ifelse(Meta$Treatment=="Sham","red","black"),
density.info="none", trace="none", cexRow=0.7,cexCol=1.2)
## Error: `x' must have at least 2 rows and 2 columns
par(mfrow=c(1,2))
limma::plotMA(fit2,array=2,main="MA-plot")
abline(h=0)
abline(h=c(-1,1),col="red")
text(topTable_tissue[1:20,]$AveExpr,topTable_tissue[1:20,]$logFC,topTable_tissue$genes[1:20],col="red")
volcanoplot(fit2,coef=2,main="Volcano")
abline(h=-log(.05))
abline(v=c(-1,1),col="red")
#/////////////////////////////////////////////////////////////////////////////
# Ileum+Duodenum
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
par(mfrow=c(1,1))
design <- model.matrix(~Treatment)
colnames(design)
## [1] "(Intercept)" "TreatmentSham"
y <- DGEList(counts=counttable_all, genes=rownames(counttable_all))
isexpr <- rowSums(cpm(y)>1) >= 3
y <- y[isexpr,]
y <- calcNormFactors(y)
v <- voom(y,design,plot=F)
fit <- lmFit(v,design)
fit2 <- eBayes(fit)
topTable_treatment <-topTable(fit2,coef="TreatmentSham",p.value=1,n=200000,sort="p")
table(significant=topTable_treatment$adj.P.Val < .05,Upregulated=topTable_treatment$logFC < 0)
## Upregulated
## significant FALSE TRUE
## FALSE 6471 7026
## TRUE 342 643
# Annotate and write to file
is<-intersect(rownames(topTable_treatment),annotation$ensembl_gene_id)
topTable_treatment_annotated <- cbind(topTable_treatment[is,-1],annotation[is,])
head(topTable_treatment_annotated,n=20)
## logFC AveExpr t P.Value adj.P.Val B
## ENSSSCG00000024018 -1.7322 4.884 -8.178 8.562e-08 0.0008048 8.053
## ENSSSCG00000002524 -1.8229 4.220 -8.041 1.111e-07 0.0008048 7.739
## ENSSSCG00000015120 -1.7817 3.539 -7.499 3.207e-07 0.0015481 6.687
## ENSSSCG00000017983 -1.0648 6.238 -7.253 5.254e-07 0.0015601 6.406
## ENSSSCG00000005916 -1.0791 3.255 -7.241 5.386e-07 0.0015601 6.171
## ENSSSCG00000007707 -1.3602 1.613 -7.126 6.810e-07 0.0016437 5.621
## ENSSSCG00000027903 -1.6148 3.891 -6.631 1.900e-06 0.0036557 5.123
## ENSSSCG00000003083 -1.3045 5.356 -6.554 2.236e-06 0.0036557 5.030
## ENSSSCG00000012257 -1.0047 8.756 -6.519 2.412e-06 0.0036557 4.915
## ENSSSCG00000015567 1.3022 3.780 6.400 3.108e-06 0.0036557 4.621
## ENSSSCG00000027197 -0.9475 5.550 -6.317 3.708e-06 0.0036557 4.551
## ENSSSCG00000018007 -1.5778 1.522 -6.291 3.924e-06 0.0036557 4.154
## ENSSSCG00000022884 -1.7603 6.055 -6.290 3.932e-06 0.0036557 4.495
## ENSSSCG00000003374 -1.4888 5.434 -6.286 3.968e-06 0.0036557 4.486
## ENSSSCG00000010009 -1.4261 5.346 -6.270 4.101e-06 0.0036557 4.454
## ENSSSCG00000027946 -0.9883 8.507 -6.252 4.266e-06 0.0036557 4.397
## ENSSSCG00000003302 -3.0709 3.308 -6.249 4.291e-06 0.0036557 4.329
## ENSSSCG00000005394 -2.1310 11.000 -6.182 4.961e-06 0.0037816 4.090
## ENSSSCG00000017495 -1.0485 4.781 -6.182 4.961e-06 0.0037816 4.270
## ENSSSCG00000013599 -1.7710 5.499 -6.116 5.729e-06 0.0038583 4.137
## hgnc_symbol ensembl_gene_id
## ENSSSCG00000024018 SLC16A3 ENSSSCG00000024018
## ENSSSCG00000002524 AMN ENSSSCG00000002524
## ENSSSCG00000015120 ENSSSCG00000015120
## ENSSSCG00000017983 PER1 ENSSSCG00000017983
## ENSSSCG00000005916 KIAA1875 ENSSSCG00000005916
## ENSSSCG00000007707 VPS37D ENSSSCG00000007707
## ENSSSCG00000027903 ENSSSCG00000027903
## ENSSSCG00000003083 ENSSSCG00000003083
## ENSSSCG00000012257 MAOA ENSSSCG00000012257
## ENSSSCG00000015567 FAM129A ENSSSCG00000015567
## ENSSSCG00000027197 DYRK1B ENSSSCG00000027197
## ENSSSCG00000018007 MYH3 ENSSSCG00000018007
## ENSSSCG00000022884 ENSSSCG00000022884
## ENSSSCG00000003374 ENSSSCG00000003374
## ENSSSCG00000010009 ENSSSCG00000010009
## ENSSSCG00000027946 MVP ENSSSCG00000027946
## ENSSSCG00000003302 ENSSSCG00000003302
## ENSSSCG00000005394 ENSSSCG00000005394
## ENSSSCG00000017495 ENSSSCG00000017495
## ENSSSCG00000013599 ENSSSCG00000013599
## description
## ENSSSCG00000024018 solute carrier family 16 (monocarboxylate transporter), member 3 [Source:HGNC Symbol;Acc:10924]
## ENSSSCG00000002524 amnion associated transmembrane protein [Source:HGNC Symbol;Acc:14604]
## ENSSSCG00000015120 Ubiquitin carboxyl-terminal hydrolase [Source:UniProtKB/TrEMBL;Acc:F1SAF6]
## ENSSSCG00000017983 period circadian clock 1 [Source:HGNC Symbol;Acc:8845]
## ENSSSCG00000005916 KIAA1875 [Source:HGNC Symbol;Acc:26959]
## ENSSSCG00000007707 vacuolar protein sorting 37 homolog D (S. cerevisiae) [Source:HGNC Symbol;Acc:18287]
## ENSSSCG00000027903 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LHS9]
## ENSSSCG00000003083 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RMS6]
## ENSSSCG00000012257 monoamine oxidase A [Source:HGNC Symbol;Acc:6833]
## ENSSSCG00000015567 family with sequence similarity 129, member A [Source:HGNC Symbol;Acc:16784]
## ENSSSCG00000027197 dual-specificity tyrosine-(Y)-phosphorylation regulated kinase 1B [Source:HGNC Symbol;Acc:3092]
## ENSSSCG00000018007 myosin, heavy chain 3, skeletal muscle, embryonic [Source:HGNC Symbol;Acc:7573]
## ENSSSCG00000022884 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LT50]
## ENSSSCG00000003374 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RIL6]
## ENSSSCG00000010009 Sus scrofa galactose-3-O-sulfotransferase 1 (GAL3ST1), mRNA. [Source:RefSeq mRNA;Acc:NM_001244429]
## ENSSSCG00000027946 major vault protein [Source:HGNC Symbol;Acc:7531]
## ENSSSCG00000003302 Sus scrofa transmembrane protein 86B (TMEM86B), mRNA. [Source:RefSeq mRNA;Acc:NM_001243276]
## ENSSSCG00000005394 Fructose-bisphosphate aldolase [Source:UniProtKB/TrEMBL;Acc:F1SSB5]
## ENSSSCG00000017495 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RX94]
## ENSSSCG00000013599 Sus scrofa angiopoietin-like 4 (ANGPTL4), mRNA. [Source:RefSeq mRNA;Acc:NM_001038644]
write.csv(file="topTable_treatment.csv",topTable_treatment_annotated)
par(mfrow=c(1,1))
for (gene in c(topTable_treatment$genes[1:9])) {
boxplot(v$E[gene,]~Meta$Treatment,ylab=gene,col=c("grey","steelblue2"),
main=paste(topTable_treatment_annotated[gene,]$hgnc_symbol, "\nAdj. p-val:",prettyNum(topTable_treatment[gene,]$adj.P.Val)),cex=2)
}
heatmap.2(v$E[topTable_treatment[topTable_treatment$adj.P.Val < .01,]$genes,],col=redblue(256),
scale="row",key=TRUE, keysize=1.5,
ColSideColors=ifelse(Meta$Treatment=="Sham","red","black"),
density.info="none", trace="none", cexRow=0.7,cexCol=1.2)
par(mfrow=c(1,2))
limma::plotMA(fit2,array=2,main="MA-plot")
abline(h=0)
abline(h=c(-1,1),col="red")
text(topTable_treatment[1:20,]$AveExpr,topTable_treatment[1:20,]$logFC,topTable_treatment$genes[1:20],col="red")
volcanoplot(fit2,coef=2,main="Volcano")
abline(h=-log(.05))
abline(v=c(-1,1),col="red")
#/////////////////////////////////////////////////////////////////////////////
# Separate
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
par(mfrow=c(1,1))
Treatment<-Meta$Treatment
Tissue <- Meta$Tissue
TreatTissue <- paste(Treatment,Tissue,sep="_")
design <- model.matrix(~0+TreatTissue)
colnames(design) <- gsub("TreatTissue","",colnames(design))
y <- DGEList(counts=counttable_all, genes=rownames(counttable_all))
isexpr <- rowSums(cpm(y)>1) >= 3
y <- y[isexpr,]
y <- calcNormFactors(y)
v <- voom(y,design,plot=T)
fit <- lmFit(v,design)
cont.matrix <- makeContrasts(
Treat_ileum=GBP_ileum-Sham_ileum,
Treat_duodenum=GBP_duodenum-Sham_duodenum,
Diff=(Sham_ileum-GBP_ileum)-(Sham_duodenum-GBP_duodenum),
levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
results <- decideTests(fit2,adjust.method="fdr")
vennDiagram(results,include="both")
#/////////////////////////////////////////////////////////////////////////////
# Ileum
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
topTable_ileum <-topTable(fit2,coef="Treat_ileum",p.value=1,n=200000,sort="p")
table(significant=topTable_ileum$adj.P.Val < .25,Upregulated=topTable_ileum$logFC > 0)
## Upregulated
## significant FALSE TRUE
## FALSE 6422 7951
## TRUE 20 89
# Annotate and write to file
rownames(annotation) <- annotation$ensembl_gene_id
is<-intersect(rownames(topTable_ileum),annotation$ensembl_gene_id)
topTable_ileum_annotated <- cbind(topTable_ileum[is,-1],annotation[is,])
head(topTable_ileum_annotated,n=20)
## logFC AveExpr t P.Value adj.P.Val B
## ENSSSCG00000013393 -0.9826 5.559 -5.923 1.295e-05 0.09864 3.2786
## ENSSSCG00000023038 1.8903 3.456 5.708 2.021e-05 0.09864 2.6423
## ENSSSCG00000030198 1.7435 3.306 5.543 2.859e-05 0.09864 2.3333
## ENSSSCG00000008886 1.5372 3.091 5.449 3.488e-05 0.09864 2.1141
## ENSSSCG00000002524 1.8098 4.220 5.332 4.483e-05 0.09864 2.1081
## ENSSSCG00000024018 1.6959 4.884 5.329 4.507e-05 0.09864 2.1568
## ENSSSCG00000020695 2.1722 1.576 5.259 5.234e-05 0.09864 1.4399
## ENSSSCG00000013366 1.2922 9.062 5.204 5.889e-05 0.09864 1.7849
## ENSSSCG00000017296 2.8639 5.710 5.097 7.426e-05 0.09864 1.7091
## ENSSSCG00000028814 2.2127 3.277 5.016 8.849e-05 0.09864 1.4508
## ENSSSCG00000017068 4.1442 2.513 5.011 8.935e-05 0.09864 1.2892
## ENSSSCG00000003374 1.5963 5.434 4.992 9.322e-05 0.09864 1.5262
## ENSSSCG00000005484 -1.3340 2.227 -4.991 9.341e-05 0.09864 1.0772
## ENSSSCG00000028873 -0.7727 5.712 -4.924 1.079e-04 0.09864 1.3990
## ENSSSCG00000016328 1.5869 4.270 4.899 1.139e-04 0.09864 1.3155
## ENSSSCG00000000916 -0.7433 7.304 -4.883 1.181e-04 0.09864 1.3141
## ENSSSCG00000017163 6.5191 1.757 4.865 1.229e-04 0.09864 0.9274
## ENSSSCG00000013017 0.7598 6.574 4.848 1.274e-04 0.09864 1.2492
## ENSSSCG00000027903 1.6738 3.891 4.813 1.375e-04 0.09864 1.1212
## ENSSSCG00000028753 1.9381 4.295 4.787 1.456e-04 0.09864 1.1029
## hgnc_symbol ensembl_gene_id
## ENSSSCG00000013393 ENSSSCG00000013393
## ENSSSCG00000023038 ENSSSCG00000023038
## ENSSSCG00000030198 ENSSSCG00000030198
## ENSSSCG00000008886 ENSSSCG00000008886
## ENSSSCG00000002524 AMN ENSSSCG00000002524
## ENSSSCG00000024018 SLC16A3 ENSSSCG00000024018
## ENSSSCG00000020695 ENSSSCG00000020695
## ENSSSCG00000013366 ENSSSCG00000013366
## ENSSSCG00000017296 ENSSSCG00000017296
## ENSSSCG00000028814 SOD3 ENSSSCG00000028814
## ENSSSCG00000017068 FAXDC2 ENSSSCG00000017068
## ENSSSCG00000003374 ENSSSCG00000003374
## ENSSSCG00000005484 ZNF618 ENSSSCG00000005484
## ENSSSCG00000028873 ENSSSCG00000028873
## ENSSSCG00000016328 ENSSSCG00000016328
## ENSSSCG00000000916 LUM ENSSSCG00000000916
## ENSSSCG00000017163 ENPP7 ENSSSCG00000017163
## ENSSSCG00000013017 EHD1 ENSSSCG00000013017
## ENSSSCG00000027903 ENSSSCG00000027903
## ENSSSCG00000028753 ENSSSCG00000028753
## description
## ENSSSCG00000013393 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1S981]
## ENSSSCG00000023038 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LU87]
## ENSSSCG00000030198 Putative envelope polyprotein; Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:B3VQ66]
## ENSSSCG00000008886 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LJ75]
## ENSSSCG00000002524 amnion associated transmembrane protein [Source:HGNC Symbol;Acc:14604]
## ENSSSCG00000024018 solute carrier family 16 (monocarboxylate transporter), member 3 [Source:HGNC Symbol;Acc:10924]
## ENSSSCG00000020695 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LM85]
## ENSSSCG00000013366 Sus scrofa lactate dehydrogenase A (LDHA), mRNA. [Source:RefSeq mRNA;Acc:NM_001172363]
## ENSSSCG00000017296 Sus scrofa angiotensin I converting enzyme (peptidyl-dipeptidase A) 1 (ACE), transcript variant 2, mRNA. [Source:RefSeq mRNA;Acc:NM_001083941]
## ENSSSCG00000028814 superoxide dismutase 3, extracellular [Source:HGNC Symbol;Acc:11181]
## ENSSSCG00000017068 fatty acid hydroxylase domain containing 2 [Source:HGNC Symbol;Acc:1334]
## ENSSSCG00000003374 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RIL6]
## ENSSSCG00000005484 zinc finger protein 618 [Source:HGNC Symbol;Acc:29416]
## ENSSSCG00000028873 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LV47]
## ENSSSCG00000016328 Sus scrofa RAB17, member RAS oncogene family (RAB17), mRNA. [Source:RefSeq mRNA;Acc:NM_001243950]
## ENSSSCG00000000916 lumican [Source:HGNC Symbol;Acc:6724]
## ENSSSCG00000017163 ectonucleotide pyrophosphatase/phosphodiesterase 7 [Source:HGNC Symbol;Acc:23764]
## ENSSSCG00000013017 EH-domain containing 1 [Source:HGNC Symbol;Acc:3242]
## ENSSSCG00000027903 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:I3LHS9]
## ENSSSCG00000028753 Sus scrofa cytochrome P450, family 4, subfamily F, polypeptide 55 (CYP4F55), mRNA. [Source:RefSeq mRNA;Acc:NM_001244636]
write.csv(file="topTable_ileum.csv",topTable_ileum_annotated)
par(mfrow=c(3,3))
for (gene in c(topTable_ileum$genes[1:9])) {
boxplot(v$E[gene,]~paste(Meta$Treatment,Meta$Tissue),ylab=gene,col=c("grey","steelblue2"),
main=paste(topTable_ileum_annotated[gene,]$hgnc_symbol, "\nAdj. p-val:",prettyNum(topTable_ileum[gene,]$adj.P.Val)),cex=2)
}
heatmap.2(v$E[topTable_ileum[topTable_ileum$adj.P.Val < .25,]$genes,Meta$Tissue=="ileum"],col=redblue(256),
scale="row",key=TRUE, keysize=1.5,
ColSideColors=ifelse(Meta[Meta$Tissue=="ileum",]$Treatment=="Sham","red","black"),
density.info="none", trace="none", cexRow=0.7,cexCol=1.2)
par(mfrow=c(1,2))
limma::plotMA(fit2,array=2,main="MA-plot")
abline(h=0)
abline(h=c(-1,1),col="red")
text(topTable_ileum[1:20,]$AveExpr,topTable_ileum[1:20,]$logFC,topTable_ileum$genes[1:20],col="red")
volcanoplot(fit2,coef=2,main="Volcano")
abline(h=-log(.05))
abline(v=c(-1,1),col="red")
#/////////////////////////////////////////////////////////////////////////////
# Duodenum
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
par(mfrow=c(1,1))
topTable_duodenum <-topTable(fit2,coef="Treat_duodenum",p.value=1,n=200000,sort="p")
table(significant=topTable_duodenum$adj.P.Val < .1,Upregulated=topTable_duodenum$logFC > 0)
## Upregulated
## significant FALSE TRUE
## FALSE 7176 7278
## TRUE 1 27
# Annotate and write to file
is<-intersect(rownames(topTable_duodenum),annotation$ensembl_gene_id)
topTable_duodenum_annotated <- cbind(topTable_duodenum[is,-1],annotation[is,])
head(topTable_duodenum_annotated,n=20)
## logFC AveExpr t P.Value adj.P.Val B
## ENSSSCG00000017983 1.2698 6.238 6.311 5.868e-06 0.04520 4.024
## ENSSSCG00000015120 2.0460 3.539 6.091 9.166e-06 0.04520 3.267
## ENSSSCG00000010258 1.1249 4.965 6.080 9.362e-06 0.04520 3.550
## ENSSSCG00000005916 1.2420 3.255 5.747 1.863e-05 0.06744 2.653
## ENSSSCG00000024018 1.7537 4.884 5.616 2.452e-05 0.07101 2.689
## ENSSSCG00000002524 1.8188 4.220 5.431 3.630e-05 0.08383 2.280
## ENSSSCG00000009399 -1.0341 2.354 -5.379 4.052e-05 0.08383 1.715
## ENSSSCG00000027197 1.0710 5.550 5.238 5.477e-05 0.09548 2.010
## ENSSSCG00000003083 1.4904 5.356 5.201 5.934e-05 0.09548 1.933
## ENSSSCG00000013599 2.0644 5.499 5.091 7.515e-05 0.09570 1.717
## ENSSSCG00000027946 1.0899 8.507 5.020 8.771e-05 0.09570 1.545
## ENSSSCG00000021536 1.4725 7.016 4.981 9.543e-05 0.09570 1.514
## ENSSSCG00000012852 1.5248 8.774 4.940 1.044e-04 0.09570 1.379
## ENSSSCG00000013242 1.0811 6.713 4.893 1.155e-04 0.09570 1.344
## ENSSSCG00000011831 2.2145 3.775 4.873 1.207e-04 0.09570 1.255
## ENSSSCG00000012257 1.0409 8.756 4.847 1.278e-04 0.09570 1.202
## ENSSSCG00000000061 1.7357 4.808 4.832 1.319e-04 0.09570 1.195
## ENSSSCG00000002403 0.8313 5.456 4.813 1.375e-04 0.09570 1.180
## ENSSSCG00000009836 1.2824 4.037 4.784 1.465e-04 0.09570 1.054
## ENSSSCG00000007710 1.3234 3.778 4.780 1.478e-04 0.09570 1.027
## hgnc_symbol ensembl_gene_id
## ENSSSCG00000017983 PER1 ENSSSCG00000017983
## ENSSSCG00000015120 ENSSSCG00000015120
## ENSSSCG00000010258 AIFM2 ENSSSCG00000010258
## ENSSSCG00000005916 KIAA1875 ENSSSCG00000005916
## ENSSSCG00000024018 SLC16A3 ENSSSCG00000024018
## ENSSSCG00000002524 AMN ENSSSCG00000002524
## ENSSSCG00000009399 CYSLTR2 ENSSSCG00000009399
## ENSSSCG00000027197 DYRK1B ENSSSCG00000027197
## ENSSSCG00000003083 ENSSSCG00000003083
## ENSSSCG00000013599 ENSSSCG00000013599
## ENSSSCG00000027946 MVP ENSSSCG00000027946
## ENSSSCG00000021536 ENSSSCG00000021536
## ENSSSCG00000012852 CDHR5 ENSSSCG00000012852
## ENSSSCG00000013242 ENSSSCG00000013242
## ENSSSCG00000011831 ENSSSCG00000011831
## ENSSSCG00000012257 MAOA ENSSSCG00000012257
## ENSSSCG00000000061 PMM1 ENSSSCG00000000061
## ENSSSCG00000002403 VIPAS39 ENSSSCG00000002403
## ENSSSCG00000009836 ENSSSCG00000009836
## ENSSSCG00000007710 ENSSSCG00000007710
## description
## ENSSSCG00000017983 period circadian clock 1 [Source:HGNC Symbol;Acc:8845]
## ENSSSCG00000015120 Ubiquitin carboxyl-terminal hydrolase [Source:UniProtKB/TrEMBL;Acc:F1SAF6]
## ENSSSCG00000010258 apoptosis-inducing factor, mitochondrion-associated, 2 [Source:HGNC Symbol;Acc:21411]
## ENSSSCG00000005916 KIAA1875 [Source:HGNC Symbol;Acc:26959]
## ENSSSCG00000024018 solute carrier family 16 (monocarboxylate transporter), member 3 [Source:HGNC Symbol;Acc:10924]
## ENSSSCG00000002524 amnion associated transmembrane protein [Source:HGNC Symbol;Acc:14604]
## ENSSSCG00000009399 cysteinyl leukotriene receptor 2 [Source:HGNC Symbol;Acc:18274]
## ENSSSCG00000027197 dual-specificity tyrosine-(Y)-phosphorylation regulated kinase 1B [Source:HGNC Symbol;Acc:3092]
## ENSSSCG00000003083 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RMS6]
## ENSSSCG00000013599 Sus scrofa angiopoietin-like 4 (ANGPTL4), mRNA. [Source:RefSeq mRNA;Acc:NM_001038644]
## ENSSSCG00000027946 major vault protein [Source:HGNC Symbol;Acc:7531]
## ENSSSCG00000021536 Sus scrofa claudin 3 (CLDN3), mRNA. [Source:RefSeq mRNA;Acc:NM_001160075]
## ENSSSCG00000012852 cadherin-related family member 5 [Source:HGNC Symbol;Acc:7521]
## ENSSSCG00000013242 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1SIC2]
## ENSSSCG00000011831 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1SQX9]
## ENSSSCG00000012257 monoamine oxidase A [Source:HGNC Symbol;Acc:6833]
## ENSSSCG00000000061 phosphomannomutase 1 [Source:HGNC Symbol;Acc:9114]
## ENSSSCG00000002403 VPS33B interacting protein, apical-basolateral polarity regulator, spe-39 homolog [Source:HGNC Symbol;Acc:20347]
## ENSSSCG00000009836 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RL48]
## ENSSSCG00000007710 Uncharacterized protein [Source:UniProtKB/TrEMBL;Acc:F1RJN3]
write.csv(file="topTable_duodenum.csv",topTable_duodenum_annotated)
par(mfrow=c(3,3))
for (gene in c(topTable_duodenum$genes[1:9])) {
boxplot(v$E[gene,]~paste(Meta$Treatment,Meta$Tissue),ylab=gene,col=c("grey","steelblue2"),
main=paste(topTable_duodenum_annotated[gene,]$hgnc_symbol, "\nAdj. p-val:",prettyNum(topTable_duodenum[gene,]$adj.P.Val)),cex=2)
}
heatmap.2(v$E[topTable_duodenum[topTable_duodenum$adj.P.Val < .10,]$genes,Meta$Tissue=="duodenum"],col=redblue(256),
scale="row",key=TRUE, keysize=1.5,
ColSideColors=ifelse(Meta[Meta$Tissue=="duodenum",]$Treatment=="Sham","red","black"),
density.info="none", trace="none", cexRow=0.7,cexCol=1.2)
par(mfrow=c(1,2))
limma::plotMA(fit2,array=2,main="MA-plot")
abline(h=0)
abline(h=c(-1,1),col="red")
text(topTable_duodenum[1:20,]$AveExpr,topTable_duodenum[1:20,]$logFC,topTable_duodenum$genes[1:20],col="red")
volcanoplot(fit2,coef=2,main="Volcano")
abline(h=-log(.05))
abline(v=c(-1,1),col="red")
