Figure1.sh

for f in ../data/pats/*; do cat $f | wc -l; done | Rscript -e 'data=scan("stdin");pdf("../results/countDensity.pdf",width=10,height=6);par(mar=c(5.1,5.1,2.1,2.1));dens=density(log10(data),adjust=1.2);print(10^dens$x[which.max(dens$y)]);plot(dens,xlab="log10(Number of Patients)",main="",cex.lab=1.5,lwd=4,ylim=c(0,.46),xlim=c(0,log10(max(data))),xaxs="i",yaxs="i");abline(v=log10(251),col="red",lwd=4,lty=2);text(2.8,0.4,"1 in 2,000\nprevalence",col="red");dev.off()'

join <(cut -f1,2 ../data/dis.data | sort -k1,1) <(sort -k1,1 ../data/optum.counts) | awk '{print $1"\t"$2/502493*100"\t"$3/38734971*100}' | cut -f2,3 | Rscript -e 'data=read.table("stdin");pdf("../results/compareOptum.pdf",width=10,height=6);par(mar=c(5.1,5.1,2.1,2.1));plot(data[,1],data[,2],log="xy",xlab="Prevalence in UK Biobank (%)",ylab="Prevalence in Optum (%)",cex.lab=1.5);abline(lm(log10(data[,2])~log10(data[,1])),lwd=3,lty=2,col="red");text(.1,2e-3,paste0("r=",round(cor(data[,1],data[,2]),2)),col="red",cex=2);dev.off();print(cor.test(data[,1],data[,2])$p.val)'

for f in ../data/pats/*; do basename $f; done | awk 'NR==FNR{orpha[$1]=1;next} ($1 in orpha)' - ../results/diseases.txt | cut -f4 | sort | uniq -c | sort -nk1,1 | while read n cat; do printf "$n\t$cat\n"; done | awk -F'\t' 'NR==FNR{pats[$3]=$1+$2;next} {print $2"\t"$1"\t"pats[$2]}' ../data/sexGroups.counts - | Rscript -e 'data=read.table("stdin",sep="\t");pdf("../results/catBar.pdf",width=10,height=6);par(mar=c(11.1,5.1,1.1,7.1));barplot(data[,2],col="lightblue",ylab="Number of Diseases",cex.lab=2,cex.axis=1.5);par(new=TRUE);par(lwd=4);barplot(data[,3],col=NA,border="red",axes=FALSE);mtext("Number of Individuals",side=4,col="red",line=5,cex=2);axis(4,col="red",col.axis="red",las=1,cex.axis=1.5);text(seq(.5,19.86,by=1.21),0,srt=60,adj=c(1,1),xpd=TRUE,labels=data[,1],cex=1.5);dev.off()'

cut -f3,4 ../data/dis.data | awk '$1!="NA"' | Rscript -e 'library(scales);library(ggplot2);data=read.table("stdin");data[,1]=data[,1]*100;print(cor.test(data[,1],data[,2]));colnames(data)=c("x","y");pdf("../results/agesex.pdf",width=10,height=6);ggplot(data,aes(x=x,y=y))+geom_hex(bins=20)+scale_fill_continuous(type="viridis",breaks=pretty_breaks(),alpha=0.4)+geom_point(shape=8)+theme_bw()+xlab("Percentage Male (%)")+ylab("Mean Age at Recruitment")+labs(fill="Number of\nDiseases")+theme(text = element_text(size = 20))+geom_vline(xintercept=45.5957,col="red",linetype="dashed")+geom_hline(yintercept=56.5286,col="red",linetype="dashed");dev.off()'

cut -f5- ../data/dis.data | awk '$1!="NA"' | Rscript -e 'data=read.table("stdin");colnames(data)=c("Digestive","Circulatory","Musculoskeletal","Genitourinary","Endocrine/metabolic","Neurological","Other Immune","Respiratory","Injury/Poisoning","Eye/Ear","Neoplasms","Skin/Subcutaneous","Infectious/Parasitic","Pregnancy/Childbirth","Blood","Congenital");pdf("../results/comorbid.pdf",width=18,height=6);par(mar=c(7.1,5.1,2.1,2.1));data=data[,order(-apply(data,2,median))]*100;boxplot(data,outline=FALSE,ylab="Percentage of Patients (%)",xaxt="n",cex.lab=1.5);for(idx in 1:ncol(data))points(jitter(rep(idx,nrow(data)),amount=0.3),data[,idx],pch=16,col=rgb(0,0,0,0.2));text(seq(1,16,length.out=ncol(data)),-7,srt=45,adj=1,xpd=TRUE,labels=colnames(data));dev.off()'