Skip to content

Instantly share code, notes, and snippets.

@CnrLwlss
Created May 25, 2018 09:59
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save CnrLwlss/7f7262ca1d9ecaf9787f7322ca089567 to your computer and use it in GitHub Desktop.
Save CnrLwlss/7f7262ca1d9ecaf9787f7322ca089567 to your computer and use it in GitHub Desktop.
Generating some exploratory plots of mitocyto data.
# Read data and rename columns
dat = read.delim("mitocyto_merged_results.csv",sep=",",stringsAsFactors=FALSE)
colnames(dat) = c("value","id","channel","patient_id","patient_type")
# Specify which ids correspond to patients, and which to control
dat$patient_type = ifelse(dat$patient_id=="M1105","patient","control")
dat$patient_id = paste(toupper(substr(dat$patient_type,1,1)),dat$patient_id,sep="_")
# Specify some colours for plotting
dat$colour = "black"
dat$colour[dat$patient_type=="patient"]=rgb(1,0,0,0.3)
dat$colour[dat$patient_type=="control"]=rgb(0,0,1,0.3)
# Different ways to group the data
channels = sort(unique(dat$channel))
patient_types = sort(unique(dat$patient_type))
patient_ids = sort(unique(dat$patient_id))
# Stripchart reports, one per channel (compare between individuals)
pdf("ByChannel.pdf",width=16,height=12,pointsize=20)
op = par(mfrow=c(3,4),mar = c(3, 4, 0, 0.3) + 0.1)
for(channel in channels){
d = dat[dat$channel == channel,]
xgrps = levels(factor(d$patient_id))
xcols = sapply(xgrps,function(x) unique(d$col[d$patient_id==x]))
stripchart(d$value~d$patient_id,col=xcols,vertical=TRUE,method="jitter",ylab=channel,pch=16,cex=1)
}
par(op)
dev.off()
# Stripchart reports, one per individual (compare between channels)
pdf("ByPatient.pdf",width=8,height=8,pointsize=14)
op = par(mar = c(8, 4, 2, 2) + 0.1)
for(patient_id in patient_ids){
d = dat[dat$patient_id == patient_id,]
xgrps = levels(factor(d$channels))
xcols = rgb(0,0,0,0.05)
stripchart(d$value~d$channel,col=xcols,vertical=TRUE,
method="jitter",jitter=0.2,ylab=patient_id,pch=16,las=2,
ylim=c(min(dat$value),max(dat$value)),cex=2,
main=paste("N =",length(d$value)))
}
par(op)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment