AhernTBL_annotation_v01_pass.R

df.Tuk2017<-read.csv("Tukiainen2017_Xchr_Escape.csv", row.names = 1)  
rownames(df.Tuk2017)<-t(matrix(unlist(strsplit(rownames(df.Tuk2017), "[.]")),nrow=2))[,1]
df.Tuk2017[which(df.Tuk2017[,length(colnames(df.Tuk2017))] == "x"),"Reported_XCI_status"]<-"Variable"

df.chhabra<-read.csv("ChhabraWarmflash2021f5.csv", header=T)
colnames(df.chhabra) <- c("gene", "Chhabra")
df.xiang<-read.csv("XiangLi2019s4.csv", header=T)
colnames(df.xiang) <- c("gene", "Xiang")
df.west<-read.csv("WestRoberts2018s1.csv", header=T)
colnames(df.west) <- c("gene", "West", "WestP")

ls.khan<-as.list(read.csv("KhanRoberts2021s3.csv", header=T))
df.zhou<-read.csv("ZhouTang2019s3.csv", header=T)
colnames(df.zhou) <- c("gene", "myAUC", "avg_diff", "power", "pct.1", "pct.2", "cluster")

df.hsiao<-read.csv("HsiaoGilad2020s1.csv", header=T)
colnames(df.hsiao) <- c("ENSEMBL", "gene", "Phase")

names(ls.khan)<-paste0("c",substr(names(ls.khan),9,9), ".K")
ls.khan<-lapply(ls.khan, function(X) {as.vector(unlist(X))[as.vector(unlist(X))!=""]})
mc.khan<-cbind(unlist(ls.khan),rep(names(sapply(ls.khan, length)), sapply(ls.khan, length)))
rownames(mc.khan)<-mc.khan[,1]
df.khan<-as.data.frame(mc.khan)
colnames(df.khan)<-c("gene", "Khan")
df.khan<-df.khan[which(!duplicated(df.khan$gene)),]
df.chhabra<-df.chhabra[which(!duplicated(df.chhabra$gene)),]
df.xiang<-df.xiang[which(!duplicated(df.xiang$gene)),]
df.west<-df.west[which(!duplicated(df.west$gene)),]
df.zhou<-df.zhou[which(!duplicated(df.zhou$gene)),]

colnames(df.hsiao)[c(2:3)]<-c("gene","Hsiao")
colnames(df.zhou)[7]<-"Zhou"
df.hsiao<-df.hsiao[which(!duplicated(df.hsiao$gene)),]
rownames(df.hsiao)<-df.hsiao$gene
rownames(df.chhabra)<-df.chhabra$gene
rownames(df.xiang)<-df.xiang$gene
rownames(df.west)<-df.west$gene
rownames(df.zhou)<-df.zhou$gene

df.west$West<-paste0(df.west$West, ".W")
df.zhou$Zhou<-paste0(df.zhou$Zhou, ".Z")

df.xiang[df.xiang$Xiang=="EVT","Xiang"]<-"EVT.X"
df.xiang[df.xiang$Xiang=="eEVT","Xiang"]<-"EVT.early.X"
df.xiang[df.xiang$Xiang=="postCTB","Xiang"]<-"CTB.post.X"
df.xiang[df.xiang$Xiang=="preCTB","Xiang"]<-"CTB.early.X"
df.xiang[df.xiang$Xiang=="CTB","Xiang"]<-"CTB.X"
df.xiang[df.xiang$Xiang=="eSTB","Xiang"]<-"STB.early.X"
df.xiang[df.xiang$Xiang=="STB","Xiang"]<-"STB.X"

df.chhabra[df.chhabra$Chhabra=="E-AM","Chhabra"]<-"E-AM.C"
df.chhabra[df.chhabra$Chhabra=="EPI","Chhabra"]<-"EPI.C"
df.chhabra[df.chhabra$Chhabra=="TE","Chhabra"]<-"TE.C"

df.hsiao[df.hsiao$Hsiao=="G1/S","Hsiao"]<-".G1/S"
df.hsiao[df.hsiao$Hsiao=="G2","Hsiao"]<-".G2"
df.hsiao[df.hsiao$Hsiao=="M","Hsiao"]<-".M"
df.hsiao[df.hsiao$Hsiao=="M/G1","Hsiao"]<-".M/G1"
df.hsiao[df.hsiao$Hsiao=="S","Hsiao"]<-".S"

datExpr0 <- as.data.frame(t(assay(vsd3)))
datExpr<-datExpr0
colnames(datExpr)<-df[rownames(t(datExpr0)),"gene_name"]

df.khan$Khan<-as.factor(df.khan$Khan)
df.west$West<-as.factor(df.west$West)

df.xiang$Xiang<-as.factor(df.xiang$Xiang)
df.chhabra$Chhabra<-as.factor(df.chhabra$Chhabra)
df.hsiao$Hsiao<-as.factor(df.hsiao$Hsiao)

require(WGCNA)
df2<-df[,c(1:(LastMetaColumn+1))]
df3<-cbind(df2, 
           binarizeCategoricalVariable(df.chhabra[df2$gene_name,"Chhabra"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = ""),
           binarizeCategoricalVariable(df.xiang[df2$gene_name,"Xiang"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = ""),
           binarizeCategoricalVariable(df.west[df2$gene_name,"West"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = ""),
           binarizeCategoricalVariable(df.khan[df2$gene_name,"Khan"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = ""),
           binarizeCategoricalVariable(df.zhou[df2$gene_name,"Zhou"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = ""),
           binarizeCategoricalVariable(df.hsiao[df2$gene_name,"Hsiao"], namePrefix = "", includeLevelVsAll = T, dropUninformative = T, includePairwise = F, nameForAll = "")
)

df3.marker<-df3[,c((LastMetaColumn+2):length(df3))]

df3.marker<-df3.marker[,order(colnames(df3.marker))] 
df.Markers.abs<-as.data.frame(apply(df3.marker, 2, function(X) {colMedians(assay(vsd3)[rownames(df3.marker[which(X==1),]),])}))
rownames(df.Markers.abs)<-colnames(assay(vsd3))

require(matrixStats)
mn.mednorm<-assay(vsd3)-rowMedians(assay(vsd3))
rownames(mn.mednorm)<-df[rownames(mn.mednorm),"gene_name"]
df.Markers.med<-as.data.frame(apply(df3.marker, 2, function(X) {colMeans(mn.mednorm[df[rownames(df3.marker[which(X==1),]),"gene_name"],])}))