#require("preprocessCore") require(tRanslatome) quantile_normalisation <- function(df){ df_rank <- apply(df,2,rank,ties.method="min") df_sorted <- data.frame(apply(df, 2, sort)) df_mean <- apply(df_sorted, 1, mean) index_to_mean <- function(my_index, my_mean){ return(my_mean[my_index]) } df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean) rownames(df_final) <- rownames(df) return(df_final) } myTranslatome=function(rna,poly,id){ poly1=poly[(rowSums(poly[,1:2])>=5)|(rowSums(poly[,3:4])>=5),]; poly2 = as.data.frame(quantile_normalisation(poly1[,1:4])); poly2$transcriptID=rownames(poly1) rna1=rna[rowSums(rna[,1:4])>=2,]; rna2 = as.data.frame(quantile_normalisation(rna1[,1:4])); rna2$transcriptID=rownames(rna1) ex=merge(rna2,poly2,by="transcriptID"); ex2=merge(rna1,poly1,by="transcriptID"); expressionMatrix0=cbind(ex[,2],ex[,3],ex[,4],ex[,5],ex[,6],ex[,7],ex[,8],ex[,9]) rownames(expressionMatrix0)=ex[,1] colnames(expressionMatrix0)=colnames(ex)[2:9] colnames(expressionMatrix0)=c("tot.wt.a", "tot.wt.b","tot.ko.a","tot.ko.b","pol.wt.a", "pol.wt.b","pol.ko.a", "pol.ko.b") #ex0=expressionMatrix0[(rowSums(expressionMatrix0[,1:4])>0)&(rowSums(expressionMatrix0[,5:8])>0),] ex0=expressionMatrix0 # translatome.analysis <- newTranslatomeDataset(quantile_normalisation(ex0), translatome.analysis <- newTranslatomeDataset( ex0, c("tot.wt.a", "tot.wt.b"), c("tot.ko.a", "tot.ko.b"), c("pol.wt.a", "pol.wt.b"), c("pol.ko.a", "pol.ko.b"), label.level= c("transcriptome", "translatome"), label.condition=c("wt", "ko")) limma.DEGs <- computeDEGs(translatome.analysis, method= "limma", mult.cor=TRUE,significance.threshold=0.05) png(paste("translatome-scatter-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) Scatterplot(limma.DEGs, outputformat="on screen",track="") dev.off() png(paste("translatome-histogram-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) Histogram(limma.DEGs,plottype ="detailed", outputformat="on screen") dev.off() ##enrichment analysis of the selected DEGs # CCEnrichment <- GOEnrichment(limma.DEGs,ontology="CC", classOfDEGs="up",test.method="elim", test.threshold = 0.05) #png(paste("translatome-GOradar-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) # Radar(CCEnrichment) #dev.off() #png(paste("translatome-GOheatmap-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) # Heatmap(CCEnrichment) #dev.off() ##performing a comparison of the biological themes enriched ##in the two levels of gene expression # CCComparison <- GOComparison(CCEnrichment) #png(paste("translatome-GOsimilarity-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) #SimilarityPlot(CCComparison) #dev.off() limma.DEGs <- computeDEGs(translatome.analysis, method= "limma", mult.cor=F,significance.threshold=0.05) png(paste("translatome-scatter-pval-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) Scatterplot(limma.DEGs, outputformat="on screen",track="") dev.off() png(paste("translatome-histogram-pval-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) Histogram(limma.DEGs,plottype ="detailed", outputformat="on screen") dev.off() ##enrichment analysis of the selected DEGs # CCEnrichment <- GOEnrichment(limma.DEGs,ontology="CC", classOfDEGs="up",test.method="elim", test.threshold = 0.05) #png(paste("translatome-GOradar-pval-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) # Radar(CCEnrichment) #dev.off() #png(paste("translatome-GOheatmap-pval-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) # Heatmap(CCEnrichment) #dev.off() ##performing a comparison of the biological themes enriched ##in the two levels of gene expression # CCComparison <- GOComparison(CCEnrichment) #png(paste("translatome-GOsimilarity-pval-",id,".png",sep=""),width=1024,height=1024,pointsize = 22) #SimilarityPlot(CCComparison) #dev.off() limma.DEGs } seqwd("/data/images/proton/DKlab/mr/Polysomes/translatome/") polyMatrix0= read.table("/data/images/proton/DKlab/mr/Polysomes/reads/stringtie/transcript_count_matrix.csv",header=T,sep=",") #conditions = factor(c("KO2","KO2","KOIFN","KOIFN","Ko0","Ko0","KoIL4","KoIL4","Ko6","Ko6","wt0","wt0IL4","wt0","wt0IL4","wt2","wt2","wt6","wt6","wtIFN","wtIFN")) # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 #,KO_0h_replicate_I,KO_0h_replicate_II,KO_2h_replicate_I,KO_2h_replicate_II,KO_6h_replicate_I,KO_6h_replicate_II,KO_IFN-gamma_replicate_II,KO_IFN_replicate_I,KO_IL4_replicate_I,KO_IL4_replicate_II,WT_0h_replicate_I,WT_0h_replicate_II,WT_2h_replicate_I,WT_2h_replicate_II,WT_6h_replicate_I,WT_6h_replicate_II,WT_IFN-gamma_replicate_I,WT_IFN-gamma_replicate_II,WT_IL4_replicate_I,WT_IL4_replicate_II rnaseqMatrix0= read.table("/data/images/proton/DKlab/mr/rnaseq/stringtie/transcript_count_matrix.csv",header=T,sep=",") # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 #conditions = factor(c("KO2","KO2","KOIFN","KOIFN","Ko0","Ko0","KoIL4","KoIL4","Ko6","Ko6","wt0","wt0IL4","wt0","wt0IL4","wt2","wt2","wt6","wt6","wtIFN","wtIFN")) col_ordering = c(11,12,1,2)#0h poly =polyMatrix0[,col_ordering]; poly$transcriptID=rownames(poly) col_ordering = c(11,13,5,6)#0h rna =rnaseqMatrix0[,col_ordering]; rna$transcriptID=rownames(rna) degs0h=myTranslatome(rna,poly,"0h") col_ordering = c(13,14,3,4)#2h poly =polyMatrix0[,col_ordering]; poly$transcriptID=rownames(poly) col_ordering = c(15,16,1,2)#2h rna =rnaseqMatrix0[,col_ordering]; rna$transcriptID=rownames(rna) degs=myTranslatome(rna,poly,"translatome-scatter-2h.png") col_ordering = c(15,16,5,6)#6h poly =polyMatrix0[,col_ordering]; poly$transcriptID=rownames(poly) col_ordering = c(17,18,9,10)#6h rna =rnaseqMatrix0[,col_ordering]; rna$transcriptID=rownames(rna) degs=myTranslatome(rna,poly,"translatome-scatter-6h.png") col_ordering = c(17,18,7,8)#IFN poly =polyMatrix0[,col_ordering]; poly$transcriptID=rownames(poly) col_ordering = c(19,20,3,4)#IFN rna =rnaseqMatrix0[,col_ordering]; rna$transcriptID=rownames(rna) degs=myTranslatome(rna,poly,"translatome-scatter-IFN.png") col_ordering = c(19,20,9,10)#IL4 poly =polyMatrix0[,col_ordering]; poly$transcriptID=rownames(poly) col_ordering = c(12,14,7,8)#IL4 rna =rnaseqMatrix0[,col_ordering]; rna$transcriptID=rownames(rna) degs=myTranslatome(rna,poly,"translatome-scatter-IL4.png") ### R code from vignette source 'tRanslatome_package.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: tRanslatome_package.Rnw:55-75 ################################################### ##loading the tRanslatome package library(tRanslatome) ##loading the training data set data(tRanslatomeSampleData) str(expressionMatrix) num [1:1000, 1:12] 5.89 5.31 4.53 6.51 7.59 ... - attr(*, "dimnames")=List of 2 ..$ : chr [1:1000] "A1CF" "A4GALT" "A4GNT" "AAAS" ... ..$ : chr [1:12] "tot.undiff.a" "tot.undiff.b" "tot.undiff.c" "tot.diff.a" ... > expressionMatrix[1,] tot.undiff.a tot.undiff.b tot.undiff.c tot.diff.a tot.diff.b tot.diff.c 5.894621 5.625169 6.367371 7.851970 7.549682 8.039144 pol.undiff.a pol.undiff.b pol.undiff.c pol.diff.a pol.diff.b pol.diff.c 6.971785 7.203912 5.797122 9.164181 8.996630 8.600587 > head(rownames(expressionMatrix)) [1] "A1CF" "A4GALT" "A4GNT" "AAAS" "AACS" "AADAC" [7] "AAGAB" "AAMP" "AANAT" "AARS" "AARSD1" "AASDHPPT" translatome.analysis <- newTranslatomeDataset(expressionMatrix, c("tot.wt.a", "tot.wt.b"), c("tot.ko.a", "tot.ko.b"), c("pol.wt.a", "pol.wt.b"), c("pol.ko.a", "pol.ko.b"), label.level= c("transcriptome", "translatome"), label.condition=c("wt", "ko")) ##identification of DEGs with the use of the limma statistical method limma.DEGs <- computeDEGs(translatome.analysis, method= "limma", mult.cor=TRUE) Scatterplot(limma.DEGs, outputformat="on screen",track="")