merge hmp and msd16s

merge_hmp_msd16s

load("HMPv35_100nt.biom_MRexperiment.rdata")
library(msd16s)

k = which(pData(obj)$bodysite == "UBERON_stool")
obj = obj[,k]


mergeData = function(o1,o2,norm=FALSE){
	o1names = rownames(o1)
	o2names = rownames(o2)

	unames = unique(c(o1names,o2names))
	ns = ncol(o1) + ncol(o2)
	nf = length(unames)
	mat = array(NA,dim=c(nf,ns))
	rownames(mat) = unames

	colnames(mat) = c(colnames(o1),colnames(o2))
	mat[o1names,1:ncol(o1)] = MRcounts(o1,norm=norm,sl=1)[o1names,]
	mat[o2names,(ncol(o1)+1):ncol(mat)] = MRcounts(o2,norm=norm,sl=1)[o2names,]
	mat
}

obj2 = obj[,-which(duplicated(pData(obj)$ANONYMIZEDNAME))]
obj2 = filterData(obj2,depth=1)
gnames = fData(obj2)[,6]
gnames = gsub("\\[","",gnames)
gnames = gsub("\\]","",gnames)
gnames = gsub("g__","",as.character(gnames))
gnames[which(gnames=="")] = "no_match"
fData(obj2)[,6] = gnames
hmpnf = normFactors(cumNorm(obj2,p=0.5))


library(msd16s)
msd16s = msd16s[,which(pData(msd16s)$Type=="Control")]
msd16s = filterData(msd16s,depth=1)
mnf = normFactors(msd16s)


mm = mergeData(aggTax(msd16s,"genus"),aggTax(obj2,"X6"))
mm[is.na(mm)] = 0

country = as.character(c(rep("USA",ncol(obj2)),as.character(pData(msd16s)$Country)))
country[country=="Gambia"] = "The Gambia"

study=ifelse(sapply(strsplit(colnames(mm),"\\."),length)>1,2,1)

pd = data.frame(samples = colnames(mm),study = study,country=country)
rownames(pd) = colnames(mm)
pd = AnnotatedDataFrame(pd)
merged = newMRexperiment(mm,phenoData = pd,normFactors=c(mnf,hmpnf))


trials = pData(merged)$study
heatmapColColors=brewer.pal(12,"Set3")[as.integer(trials)];
heatmapCols = colorRampPalette(brewer.pal(9, "RdBu"))(50)
plotMRheatmap(obj=merged,n=20,cexRow = 0.75,cexCol = 0.5,trace="none",col = heatmapCols,ColSideColors = heatmapColColors)
legend("left", legend=unique(trials),unique(trials), fill=brewer.pal(12,"Set3")[1:5],box.col=NA)