####  Loading libraries  ####

suppressPackageStartupMessages({
  library(MAST);  library(immunarch);  library(plyr); library(data.table);  library(reshape);  library(scater); library(ComplexHeatmap);
  library(edgeR);  library(ROCR);  library(scRNA.seq.funcs);  library(gridExtra);  library(NMF);  library(RColorBrewer);  library(ggplot2);
  library(dplyr);  library(MASS);  library(pscl);  library(ape);  library(factoextra);  library(NbClust);  library(igraph);  library(stringdist);
  library(tcR);  library(purrr);  library(tidyr);  library(dendextend);  library(doParallel);  library(foreach);  library(dtplyr);  library(tidyverse);
  library(ggpubr);  library(corrplot);  library(e1071);  library(bootstrap);  library(DAAG);  library(RColorBrewer);  library(ppcor);  library(factoextra);
  library(mixOmics);  library(glmnet);  library(psych);  library(caret);  library(randomForest);  library(PRROC);  library(viridis); library(grid);
  library(ggplotify);  library(ROCR);   library(jpeg);   library(lme4);   library(glmmTMB);   library(reshape);   library(protr);  library(msa);  library(Biostrings);  library(ggseqlogo)
})

# RA vs Controls & Whole Blood vs Synovial tissue

reads.sample<-readRDS("Data2analyze/ChainReadsCount.rds")
reads.sample.prop<-readRDS("Data2analyze/ChainReadsProp.rds")

x1<-read.csv("Data2analyze/Metadata.csv", header=T)
x2<-read.table("Data2analyze/Donor2RNA_Codes.txt", header=T)
ra.wk0<-as.character(x1[x1$Week=="wk0", "Sample"]);
ra.wk12<-as.character(x2$RNA2);
ctrl<-as.character(x1[x1$IMID=="CTRL", "Sample"]);
samples2keep<-c(ra.wk0,ra.wk12,ctrl)

tech.var<-read.table("Data2analyze/Technical_Var.txt", header=T)

reads.sample<-reads.sample[(reads.sample$Sample %in% samples2keep),]
reads.sample.prop<-reads.sample.prop[(reads.sample.prop$Sample %in% samples2keep),]

rs.sum<-reads.sample[,-1]
rs.sum$Reads.All<-rowSums(rs.sum)
rs.sum$Sample<-rownames(rs.sum)
colnames(rs.sum)<-c("Reads.TRA", "Reads.TRB", "Reads.TRG", "Reads.TRD", "Reads.IGH", "Reads.IGL", "Reads.IGK", "Reads.All", "Sample")

df.counts.long<-melt(reads.sample);  colnames(df.counts.long)[2]<-"Chain"
df.prop.long<-melt(reads.sample.prop);  colnames(df.prop.long)[2]<-"Chain"
metadata<-read.csv("Data2analyze/Metadata.csv", header=T)
df.counts.long.pheno<-merge(metadata, df.counts.long, by="Sample")
df.prop.long.pheno<-merge(metadata, df.prop.long, by="Sample")

df.sums.long1<-merge(df.prop.long.pheno, rs.sum, by="Sample")
df.sums.long<-merge(df.sums.long1,tech.var,by="Sample")

df.prop.long.pheno.ra<-df.prop.long.pheno[df.prop.long.pheno$IMID=="RA",]
df.prop.long.pheno.ra$ActQT<-as.numeric(as.character(df.prop.long.pheno.ra$ActQT))
df.prop.long.pheno.ctrl<-df.prop.long.pheno[df.prop.long.pheno$IMID=="CTRL",]

df.ra.igl<-df.prop.long.pheno.ra[df.prop.long.pheno.ra$Chain=="IGL",c("Sample","value")]
df.ra.igl$Chain<-rep("IGL",nrow(df.ra.igl))
df.ra.igk<-df.prop.long.pheno.ra[df.prop.long.pheno.ra$Chain=="IGK",c("Sample","value")]
df.ra.igk$Chain<-rep("IGK",nrow(df.ra.igk))
df.ra.merged<-merge(df.ra.igl, df.ra.igk, by="Sample")
df.ra.merged$Ratio<-df.ra.merged$value.y/df.ra.merged$value.x

df.ctrl.igl<-df.prop.long.pheno.ctrl[df.prop.long.pheno.ctrl$Chain=="IGL",c("Sample","value")]
df.ctrl.igl$Chain<-rep("IGL",nrow(df.ctrl.igl))
df.ctrl.igk<-df.prop.long.pheno.ctrl[df.prop.long.pheno.ctrl$Chain=="IGK",c("Sample","value")]
df.ctrl.igk$Chain<-rep("IGK",nrow(df.ctrl.igk))
df.ctrl.merged<-merge(df.ctrl.igl, df.ctrl.igk, by="Sample")
df.ctrl.merged$Ratio<-df.ctrl.merged$value.y/df.ctrl.merged$value.x

chain.usage.out<-data.frame("Chain"=rep("NA",7), "Beta"=rep(10,7), "Pvalue"=rep(10,7), stringsAsFactors=F)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  df.chain<-df.sums.long[(df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  df.chain$IMID<-as.numeric(df.chain$IMID)
  df.chain$IMID[which(df.chain$IMID==1)]<-0
  df.chain$IMID[which(df.chain$IMID==2)]<-1
  df.chain$Gender<-as.numeric(df.chain$Gender)
    
  model.adj.reads.all<-glm(IMID~value+Reads.All+Age+Gender, data=df.chain, family="binomial")
  pval<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  beta<-coefficients(summary(model.adj.reads.all))[,"Estimate"][2]
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-beta
  chain.usage.out[count,3]<-pval
}

p<-ggboxplot(df.sums.long, x="Chain", y="value", color="IMID", add="jitter", fill="IMID") +
  labs(x="", y="\nExpression (% reads)\n") + scale_color_manual(values=c("black","black")) +
  scale_fill_manual(values=c("#7DBB7D","#F69191")) + theme(legend.position="right", text = element_text(size=17)) + ylim(0,max(df.sums.long$value)*1.05) +
  annotate(geom="text", x=1, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[1,3],3)), color="black") +
  annotate(geom="text", x=2, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[2,3],3)), color="black") +
  annotate(geom="text", x=3, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[3,3],3)), color="black") +
  annotate(geom="text", x=4, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[4,3],3)), color="black") +
  annotate(geom="text", x=5, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[5,3],3)), color="black") +
  annotate(geom="text", x=6, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[6,3],4)), color="black") +
  annotate(geom="text", x=7, y=max(df.sums.long$value)*1.04, size=4.5, label=paste0("P= ",round(chain.usage.out[7,3],3)), color="black") +
  

jpeg("Output_Data/ChainUsage_RACTRL.jpeg", width=3700, height=2700, res=300)
plot(p)
dev.off()

df.reads.final<-readRDS("Data2analyze/ReadsMJLewis.rds")
ids<-unique(df.reads.final$Sample)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
df.chainusage<-data.frame("Sample"="test", "Chain"="test", "Usage"=0, stringsAsFactors=F)
count=0
for (sample in ids) {
  for (c in chains) {
    chain.reads<-df.reads.final[(df.reads.final$Sample==sample & df.reads.final$Chain==c),"Reads.Chain"]
    if (length(chain.reads)>0) {
      total.reads<-df.reads.final[(df.reads.final$Sample==sample & df.reads.final$Chain==c),"Reads.All"]
      chain.usage<-as.numeric(chain.reads)/as.numeric(total.reads)*100
    } else {
      chain.usage<-0
    }
    count=count+1
    print(paste0(sample," - ",c," - ",count))
    df.chainusage[count,]<-c(sample,c,chain.usage)
  }
}
df.chainusage$Chain<-as.factor(df.chainusage$Chain)
df.chainusage$Usage<-as.numeric(df.chainusage$Usage)
df.chainusage$key<-paste0(df.chainusage$Sample,"_",df.chainusage$Chain)

df.reads.final$key<-paste0(df.reads.final$Sample,"_",df.reads.final$Chain)
df.data<-merge(df.chainusage, df.reads.final, by="key")
table(df.data$Sample.x==df.data$Sample.y)
df.data<-df.data[,-c(5,6)]

df2plot<-df.data[df.data$Pathotype!="ungraded",]
p1<-ggboxplot(df2plot, x="Chain.x", y="Usage", color="Tissue", add="jitter", fill="Tissue") +
  labs(title="", x="", y=paste0("\nExpression (% reads)\n")) +
  scale_color_manual(values=c("black","black")) + scale_fill_manual(values=c("#E7B800","#00AFBB")) +
  theme(plot.title=element_text(hjust=0.5, size=16, face="bold"), text=element_text(size=17), legend.position="right", axis.text.x = element_text(hjust=1)) +
  stat_compare_means(aes(group=Tissue), method="wilcox.test", label="p.format", size=4.3, label.y=max(df2plot$Usage*1.20), hide.ns=F)

jpeg("Output_Data/ChainUsage_BloodSynovium.jpeg", width=3700, height=2700, res=300)
plot(p1)
dev.off()


## RA phenotypes

pdf("Output_Data/RApheno_chainusage.pdf", width=10)

# RespGM
chain.usage.out<-data.frame("Chain"=rep("NA",7), "Beta"=rep(10,7), "Pvalue"=rep(10,7), stringsAsFactors=F)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  
  df.chain$Resp<-as.numeric(df.chain$Resp)
  df.chain$Resp[which(df.chain$Resp==4)]<-0
  df.chain$Resp[which(df.chain$Resp==2)]<-1
  df.chain$Resp[which(df.chain$Resp==3)]<-1
  df.chain$Gender<-as.numeric(df.chain$Gender)
  
  model.adj.reads.all<-glm(Resp~value+Reads.All+Age+Gender, data=df.chain, family="binomial")
  pval<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  beta<-coefficients(summary(model.adj.reads.all))[,"Estimate"][2]
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-beta
  chain.usage.out[count,3]<-pval
}
df2plot<-df.sums.long[(df.sums.long$IMID=="RA" & !is.na(df.sums.long$Resp)),]
df2plot$Resp<-as.character(df2plot$Resp)
df2plot$Resp[df2plot$Resp=="MOD"]<-"Responders"
df2plot$Resp[df2plot$Resp=="GOOD"]<-"Responders"
df2plot$Resp[df2plot$Resp=="NON_RESP"]<-"Non Responders"
p<-ggboxplot(df2plot, x="Chain", y="value", color="Resp", add="jitter", fill="Resp") +
  labs(x="", y="\nExpression (% reads)\n", title="\nResponse to TNFi therapy\n") + scale_color_manual(values=c("black","black")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) + theme(legend.position="right", plot.title=element_text(hjust=0.5)) + ylim(0,max(df2plot$value)*1.05) +
  annotate(geom="text", x=1, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[1,3],2)), color="black") +
  annotate(geom="text", x=2, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[2,3],2)), color="black") +
  annotate(geom="text", x=3, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[4,3],2)), color="black") +
  annotate(geom="text", x=4, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[3,3],2)), color="black") +
  annotate(geom="text", x=5, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[5,3],2)), color="black") +
  annotate(geom="text", x=6, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[6,3],2)), color="black") +
  annotate(geom="text", x=7, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[7,3],2)), color="black")
plot(p)

# ACPA
chain.usage.out<-data.frame("Chain"=rep("NA",7), "Beta"=rep(10,7), "Pvalue"=rep(10,7), stringsAsFactors=F)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  df.chain$ACPA<-as.numeric(df.chain$ACPA)
  df.chain$ACPA[which(df.chain$ACPA==2)]<-0
  df.chain$ACPA[which(df.chain$ACPA==3)]<-1
  df.chain$Gender<-as.numeric(df.chain$Gender)
  model.adj.reads.all<-glm(ACPA~value+Reads.All+Age+Gender, data=df.chain, family="binomial")
  pval<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  beta<-coefficients(summary(model.adj.reads.all))[,"Estimate"][2]
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-beta
  chain.usage.out[count,3]<-pval
}
df2plot<-df.sums.long[(df.sums.long$IMID=="RA" & !is.na(df.sums.long$ACPA)),]
df2plot$ACPA<-as.character(df2plot$ACPA)
df2plot$ACPA[df2plot$ACPA=="Pos"]<-"Positive"
df2plot$ACPA[df2plot$ACPA=="Neg"]<-"Negative"
p<-ggboxplot(df2plot, x="Chain", y="value", color="ACPA", add="jitter", fill="ACPA") +
  labs(x="", y="\nExpression (% reads)\n", title="\nACPA\n") + scale_color_manual(values=c("black","black")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) + theme(legend.position="right", plot.title=element_text(hjust=0.5)) + ylim(0,max(df2plot$value)*1.05) +
  annotate(geom="text", x=1, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[1,3],2)), color="black") +
  annotate(geom="text", x=2, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[2,3],2)), color="black") +
  annotate(geom="text", x=3, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[4,3],2)), color="black") +
  annotate(geom="text", x=4, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[3,3],2)), color="black") +
  annotate(geom="text", x=5, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[5,3],2)), color="black") +
  annotate(geom="text", x=6, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[6,3],2)), color="black") +
  annotate(geom="text", x=7, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[7,3],2)), color="black")
plot(p)

# RF
chain.usage.out<-data.frame("Chain"=rep("NA",7), "Beta"=rep(10,7), "Pvalue"=rep(10,7), stringsAsFactors=F)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  df.chain$RF<-as.numeric(df.chain$RF)
  df.chain$RF[which(df.chain$RF==2)]<-0
  df.chain$RF[which(df.chain$RF==3)]<-1
  df.chain$Gender<-as.numeric(df.chain$Gender)
  model.adj.reads.all<-glm(RF~value+Reads.All+Age+Gender, data=df.chain, family="binomial")
  pval<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  beta<-coefficients(summary(model.adj.reads.all))[,"Estimate"][2]
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-beta
  chain.usage.out[count,3]<-pval
}
df2plot<-df.sums.long[(df.sums.long$IMID=="RA" & !is.na(df.sums.long$RF)),]
df2plot$RF<-as.character(df2plot$RF)
df2plot$RF[df2plot$RF=="Pos"]<-"Positive"
df2plot$RF[df2plot$RF=="Neg"]<-"Negative"
p<-ggboxplot(df2plot, x="Chain", y="value", color="RF", add="jitter", fill="RF") +
  labs(x="", y="\nExpression (% reads)\n", title="\nRF\n") + scale_color_manual(values=c("black","black")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) + theme(legend.position="right", plot.title=element_text(hjust=0.5)) + ylim(0,max(df2plot$value)*1.05) +
  annotate(geom="text", x=1, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[1,3],2)), color="black") +
  annotate(geom="text", x=2, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[2,3],2)), color="black") +
  annotate(geom="text", x=3, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[4,3],2)), color="black") +
  annotate(geom="text", x=4, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[3,3],2)), color="black") +
  annotate(geom="text", x=5, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[5,3],2)), color="black") +
  annotate(geom="text", x=6, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[6,3],2)), color="black") +
  annotate(geom="text", x=7, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[7,3],2)), color="black")
plot(p)

# ActB
chain.usage.out<-data.frame("Chain"=rep("NA",7), "Beta"=rep(10,7), "Pvalue"=rep(10,7), stringsAsFactors=F)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  df.chain$ActB<-as.numeric(df.chain$ActB)
  df.chain$ActB[which(df.chain$ActB==2)]<-1
  df.chain$ActB[which(df.chain$ActB==3)]<-0
  df.chain$Gender<-as.numeric(df.chain$Gender)
  model.adj.reads.all<-glm(ActB~value+Reads.All+Age+Gender, data=df.chain, family="binomial")
  pval<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  beta<-coefficients(summary(model.adj.reads.all))[,"Estimate"][2]
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-beta
  chain.usage.out[count,3]<-pval
}
df2plot<-df.sums.long[(df.sums.long$IMID=="RA" & !is.na(df.sums.long$ActB)),]
df2plot$ActB<-as.character(df2plot$ActB)
p<-ggboxplot(df2plot, x="Chain", y="value", color="ActB", add="jitter", fill="ActB") +
  labs(x="", y="\nExpression (% reads)\n",title="\nDisease Activity\n") + scale_color_manual(values=c("black","black")) +
  scale_fill_manual(values=c("#E69F00", "#56B4E9")) + theme(legend.position="right", plot.title=element_text(hjust=0.5)) + ylim(0,max(df2plot$value)*1.05) +
  annotate(geom="text", x=1, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[1,3],2)), color="black") +
  annotate(geom="text", x=2, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[2,3],2)), color="black") +
  annotate(geom="text", x=3, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[4,3],2)), color="black") +
  annotate(geom="text", x=4, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[3,3],2)), color="black") +
  annotate(geom="text", x=5, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[5,3],2)), color="black") +
  annotate(geom="text", x=6, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[6,3],2)), color="black") +
  annotate(geom="text", x=7, y=max(df2plot$value)*1.04, size=4, label=paste0("P= ",round(chain.usage.out[7,3],2)), color="black")
plot(p)
dev.off()



# ACPA
chain.usage.out<-data.frame("Chain"=rep("NA",7), "P.unadj"=rep("NA",7), "P.adj.chain.reads"=rep("NA",7), "P.adj.all.reads"=rep("NA",7))
chain.usage.out$Chain<-as.character(chain.usage.out$Chain)
chain.usage.out$P.unadj<-as.numeric(chain.usage.out$P.unadj)
chain.usage.out$P.adj.all.reads<-as.numeric(chain.usage.out$P.adj.all.reads)
chain.usage.out$P.adj.chain.reads<-as.numeric(chain.usage.out$P.adj.chain.reads)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  
  df.chain$ACPA<-as.numeric(df.chain$ACPA)
  df.chain$ACPA[which(df.chain$ACPA==3)]<-1
  df.chain$ACPA[which(df.chain$ACPA==2)]<-0
  
  model.unadj<-glm(ACPA~value, data=df.chain, family="binomial")
  pval.unadj<-coefficients(summary(model.unadj))[,"Pr(>|z|)"][2]
  
  model.adj.reads.all<-glm(ACPA~value+Reads.All, data=df.chain, family="binomial")
  pval.adj.reads.all<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  
  df.chain.filt<-df.chain[,c("ACPA","value",chain.colname)]
  model.adj.reads.chain<-glm(df.chain.filt[,1]~df.chain.filt[,2]+df.chain.filt[,3], data=df.chain.filt, family="binomial")
  pval.adj.reads.chain<-coefficients(summary(model.adj.reads.chain))[,"Pr(>|z|)"][2]  
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-pval.unadj
  chain.usage.out[count,3]<-pval.adj.reads.chain
  chain.usage.out[count,4]<-pval.adj.reads.all
}
saveRDS(chain.usage.out,"Output_Data/ChainUsage_Adj_ACPA.rds")

# RF
chain.usage.out<-data.frame("Chain"=rep("NA",7), "P.unadj"=rep("NA",7), "P.adj.chain.reads"=rep("NA",7), "P.adj.all.reads"=rep("NA",7))
chain.usage.out$Chain<-as.character(chain.usage.out$Chain)
chain.usage.out$P.unadj<-as.numeric(chain.usage.out$P.unadj)
chain.usage.out$P.adj.all.reads<-as.numeric(chain.usage.out$P.adj.all.reads)
chain.usage.out$P.adj.chain.reads<-as.numeric(chain.usage.out$P.adj.chain.reads)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  
  df.chain$RF<-as.numeric(df.chain$RF)
  df.chain$RF[which(df.chain$RF==3)]<-1
  df.chain$RF[which(df.chain$RF==2)]<-0
  
  model.unadj<-glm(RF~value, data=df.chain, family="binomial")
  pval.unadj<-coefficients(summary(model.unadj))[,"Pr(>|z|)"][2]
  
  model.adj.reads.all<-glm(RF~value+Reads.All, data=df.chain, family="binomial")
  pval.adj.reads.all<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  
  df.chain.filt<-df.chain[,c("RF","value",chain.colname)]
  model.adj.reads.chain<-glm(df.chain.filt[,1]~df.chain.filt[,2]+df.chain.filt[,3], data=df.chain.filt, family="binomial")
  pval.adj.reads.chain<-coefficients(summary(model.adj.reads.chain))[,"Pr(>|z|)"][2]  
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-pval.unadj
  chain.usage.out[count,3]<-pval.adj.reads.chain
  chain.usage.out[count,4]<-pval.adj.reads.all
}
saveRDS(chain.usage.out,"Output_Data/ChainUsage_Adj_RF.rds")

# Act
chain.usage.out<-data.frame("Chain"=rep("NA",7), "P.unadj"=rep("NA",7), "P.adj.chain.reads"=rep("NA",7), "P.adj.all.reads"=rep("NA",7))
chain.usage.out$Chain<-as.character(chain.usage.out$Chain)
chain.usage.out$P.unadj<-as.numeric(chain.usage.out$P.unadj)
chain.usage.out$P.adj.all.reads<-as.numeric(chain.usage.out$P.adj.all.reads)
chain.usage.out$P.adj.chain.reads<-as.numeric(chain.usage.out$P.adj.chain.reads)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  
  df.chain$ActQT<-as.numeric(as.character(df.chain$ActQT))
  
  model.unadj<-lm(ActQT~value, data=df.chain)
  pval.unadj<-coefficients(summary(model.unadj))[,"Pr(>|t|)"][2]
  
  model.adj.reads.all<-lm(ActQT~value+Reads.All, data=df.chain)
  pval.adj.reads.all<-coefficients(summary(model.adj.reads.all))[,"Pr(>|t|)"][2]
  
  df.chain.filt<-df.chain[,c("ActQT","value",chain.colname)]
  model.adj.reads.chain<-lm(df.chain.filt[,1]~df.chain.filt[,2]+df.chain.filt[,3], data=df.chain.filt)
  pval.adj.reads.chain<-coefficients(summary(model.adj.reads.chain))[,"Pr(>|t|)"][2]  
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-pval.unadj
  chain.usage.out[count,3]<-pval.adj.reads.chain
  chain.usage.out[count,4]<-pval.adj.reads.all
}
saveRDS(chain.usage.out,"Output_Data/ChainUsage_Adj_Act.rds")
sig.chain<-chain.usage.out[chain.usage.out$P.adj.all.reads<0.05,"Chain"]
for (c in sig.chain) {
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  df.chain$ActQT<-as.numeric(as.character(df.chain$ActQT))
  pval<-chain.usage.out[chain.usage.out$Chain==c,"P.adj.all.reads"]
  p<-ggscatter(df.chain, x="ActQT", y="value", add="reg", conf.int=T, color=cols[count]) +
    labs(title="", x="\nDisease Activity (DAS28)\n", y=paste0("\n",c," Expression (% reads)\n")) + ylim(0,5) +
    theme(plot.title=element_text(hjust=0.5, size=14, face="bold"), legend.title=element_blank(),legend.position = "bottom", axis.text=element_text(size=12)) +
    annotate(geom="text", x=min(df.chain$ActQT, na.rm=T)*1.4, y=max(df.chain$value)*1.1, size=5, label=paste0("P-value = ",round(pval,7)), color="black")
  jpeg("Output_Data/ChainUsage_Act_SigChain.jpeg", width=2300, height=2300, res=300)
  plot(p)
  dev.off()
}

# Resp
chain.usage.out<-data.frame("Chain"=rep("NA",7), "P.unadj"=rep("NA",7), "P.adj.chain.reads"=rep("NA",7), "P.adj.all.reads"=rep("NA",7))
chain.usage.out$Chain<-as.character(chain.usage.out$Chain)
chain.usage.out$P.unadj<-as.numeric(chain.usage.out$P.unadj)
chain.usage.out$P.adj.all.reads<-as.numeric(chain.usage.out$P.adj.all.reads)
chain.usage.out$P.adj.chain.reads<-as.numeric(chain.usage.out$P.adj.chain.reads)
chains<-c("TRA","TRB","TRD","TRG","IGH","IGL","IGK")
count=0
for (c in chains) {
  count=count+1
  
  df.chain<-df.sums.long[(df.sums.long$IMID=="RA" & df.sums.long$Chain==c),]
  chain.colname<-paste0("Reads.",c)
  
  df.chain<-df.chain[(df.chain$Resp=="GOOD" | df.chain$Resp=="NON_RESP"),]
  df.chain$Resp<-as.numeric(df.chain$Resp)
  df.chain$Resp[which(df.chain$Resp==4)]<-0
  df.chain$Resp[which(df.chain$Resp==2)]<-1
  
  model.unadj<-glm(Resp~value, data=df.chain, family="binomial")
  pval.unadj<-coefficients(summary(model.unadj))[,"Pr(>|z|)"][2]
  
  model.adj.reads.all<-glm(Resp~value+Reads.All, data=df.chain, family="binomial")
  pval.adj.reads.all<-coefficients(summary(model.adj.reads.all))[,"Pr(>|z|)"][2]
  
  df.chain.filt<-df.chain[,c("Resp","value",chain.colname)]
  model.adj.reads.chain<-glm(df.chain.filt[,1]~df.chain.filt[,2]+df.chain.filt[,3], data=df.chain.filt, family="binomial")
  pval.adj.reads.chain<-coefficients(summary(model.adj.reads.chain))[,"Pr(>|z|)"][2]  
  
  chain.usage.out[count,1]<-c
  chain.usage.out[count,2]<-pval.unadj
  chain.usage.out[count,3]<-pval.adj.reads.chain
  chain.usage.out[count,4]<-pval.adj.reads.all
}
saveRDS(chain.usage.out,"Output_Data/ChainUsage_Adj_RespGNR.rds")


## ChainUsagePheno_wk0wk12

df.data<-readRDS("Data2analyze/SummaryAIRR.rds");
meta<-read.csv("Data2analyze/Metadata.csv", header=T);
meta$RespGM<-meta$Resp
meta$RespGM[meta$RespGM=="MOD" & !is.na(meta$RespGM)]<-"GOOD"
meta$Resp[meta$Resp=="MOD" & !is.na(meta$RespGM)]<-NA

donor.data<-read.table("Data2analyze/Donor2RNA_Codes.txt",head=T)

# wk0
df.wk0<-donor.data[,c(1,2)];  colnames(df.wk0)[2]<-"Sample"
df.wk0$Week<-rep("0",nrow(df.wk0))
df.wk0<-merge(df.data[,c("Sample","Reads.All","ChainUsage.TRA","ChainUsage.TRB","ChainUsage.TRD","ChainUsage.TRG","ChainUsage.IGH","ChainUsage.IGL","ChainUsage.IGK")], df.wk0, by="Sample")
df.wk0<-df.wk0[,c("Sample","Donor","Week","Reads.All","ChainUsage.TRA","ChainUsage.TRB","ChainUsage.TRD","ChainUsage.TRG","ChainUsage.IGH","ChainUsage.IGL","ChainUsage.IGK")]
donor.data.filt<-donor.data[(donor.data$RNA1 %in% df.wk0$Sample),]

# wk12
df.wk12<-donor.data[,c(1,3)];  colnames(df.wk12)[2]<-"Sample"
df.wk12$Week<-rep("1",nrow(df.wk12))
df.reads<-readRDS("Data2analyze/ChainReadsCountWeek12.rds")
df.reads$Reads.All<-df.reads$TRA+df.reads$TRB+df.reads$TRD+df.reads$TRG+df.reads$IGH+df.reads$IGL+df.reads$IGK
df.wk12<-merge(df.wk12, df.reads[,c("Sample","Reads.All")], by="Sample")
df.chain.usage<-readRDS("Data2analyze/ChainReadsPropWeek12.rds")
df.wk12<-merge(df.wk12, df.chain.usage[,c("Sample","TRA","TRB","TRD","TRG","IGH","IGL","IGK")], by="Sample")
colnames(df.wk12)[5:ncol(df.wk12)]<-c("ChainUsage.TRA","ChainUsage.TRB","ChainUsage.TRD","ChainUsage.TRG","ChainUsage.IGH","ChainUsage.IGL","ChainUsage.IGK")
df.wk12<-df.wk12[(df.wk12$Sample %in% donor.data.filt$RNA2),]

# ctrls
df.ctrl<-meta[meta$IMID=="CTRL",c("Sample","Sample")];   colnames(df.ctrl)<-c("Donor","Sample")
df.ctrl$Week<-"3"
df.reads<-readRDS("Data2analyze/ChainReadsCountWeek12.rds")
df.reads$Reads.All<-df.reads$TRA+df.reads$TRB+df.reads$TRD+df.reads$TRG+df.reads$IGH+df.reads$IGL+df.reads$IGK
df.ctrl<-merge(df.ctrl, df.reads[,c("Sample","Reads.All")], by="Sample")
df.chain.usage<-readRDS("Data2analyze/ChainReadsPropWeek12.rds")
df.ctrl<-merge(df.ctrl, df.chain.usage[,c("Sample","TRA","TRB","TRD","TRG","IGH","IGL","IGK")], by="Sample")
colnames(df.ctrl)[5:ncol(df.wk12)]<-c("ChainUsage.TRA","ChainUsage.TRB","ChainUsage.TRD","ChainUsage.TRG","ChainUsage.IGH","ChainUsage.IGL","ChainUsage.IGK")

# merging wk0-wk12-ctrl
data.lgtd<-rbind(df.wk0, df.wk12, df.ctrl)
data.lgtd<-merge(meta[,-5], data.lgtd, by="Sample")
data.lgtd<-data.lgtd[order(data.lgtd$Donor),]
data.lgtd$Week<-as.factor(data.lgtd$Week)

df.lgtd<-gather(data.lgtd, condition, measurement, ChainUsage.TRA:ChainUsage.IGK, factor_key=T)
df.lgtd$condition<-gsub("ChainUsage.","",df.lgtd$condition)

df.lgtd$Week<-as.character(df.lgtd$Week)
df.lgtd$Week[df.lgtd$Week=="0"]<-"Baseline"
df.lgtd$Week[df.lgtd$Week=="1"]<-"Treated"
df.lgtd$Week[df.lgtd$Week=="3"]<-"Control"
df.lgtd$Week<-as.factor(df.lgtd$Week)

vars<-c("ChainUsage.TRA","ChainUsage.TRB","ChainUsage.TRD","ChainUsage.TRG","ChainUsage.IGH","ChainUsage.IGL","ChainUsage.IGK")
pvals.df<-data.frame("Variable"=vars, "Beta"=rep(0,length(vars)), "Pval"=rep(0,length(vars)), stringsAsFactors=F)
count=0
for (v in vars) {
  count=count+1
  var.col<-which(colnames(data.lgtd)==v)
  week.col<-which(colnames(data.lgtd)=="Week")
  covar.col<-which(colnames(data.lgtd)=="Reads.All")
  age.col<-which(colnames(data.lgtd)=="Age")
  donor.col<-which(colnames(data.lgtd)=="Donor")
  gender.col<-which(colnames(data.lgtd)=="Gender")
  data2test<-data.lgtd[data.lgtd$IMID=="RA",c(var.col,week.col,covar.col,donor.col,age.col,gender.col)]
  colnames(data2test)[1]<-"ChainUsage"
  model1<-glmer(Week~ChainUsage+Reads.All+Age+Gender+(1|Donor),  family=binomial, data=data2test)
  pval<-summary(model1)$coefficients["ChainUsage","Pr(>|z|)"]
  beta<-summary(model1)$coefficients["ChainUsage","Estimate"]
  pvals.df[count,"Beta"]<-beta
  pvals.df[count,"Pval"]<-pval
}
colnames(pvals.df)[1]<-"Chain"
pvals.df$Chain<-gsub("ChainUsage.","",pvals.df$Chain)
write.csv(pvals.df, "Output_Data/Table_ChainUsage_wk0wk12.csv", row.names=F)
