Nanostring_analysis.R

# NanoString Analisys
library(dplyr)
library(psych)
library(ggplot2)
library(ggpubr)

wdir <- "squadrito_livertumor2022_nanostring"

rawdata <- paste(wdir, "Raw_data", sep = "/")
plots <- paste(wdir, "Plots", sep = "/")
dir.create(plots, showWarnings = F)

# Leggi tutti i file CSV nella cartella e concatenali
setwd(rawdata)
all_data <- NULL
for (file_name in list.files(pattern = "*.RCC")) {
    data <- read.csv(file_name, skip = 26)
    data <- data[1:(nrow(data)-3),]
    data$file <- file_name
    data$patient <- substr(file_name,1,5)
    data$Type_sample <-substr(file_name,7,7)
    all_data <- rbind(all_data, data)
}

# TAKING ALL TOGETHER
A_data <- reshape(
    data = all_data,
    direction = "wide",
    idvar = "Name",
    timevar = "file",
    v.names = "Count",
    drop = c("Type_sample","patient"))

# Normalize quantiles
df_q <- A_data[rowSums(is.na(A_data))==0 & rowSums(A_data>4)>20,]
rbind(dim(A_data),dim(df_q))
df_q <- rbind(df_q,A_data[481,])#add IL10 data, because of importance
df_q[,4:ncol(df_q)] <- limma::normalizeQuantiles(df_q[,4:ncol(df_q)])
df_q[,4:ncol(df_q)] <- log10(df_q[,4:ncol(df_q)])

# Split in 2, according to Tumor (T) o peritumor (P) areas
typeS <- paste("type" , substr(colnames(df_q),13,13),sep=".")
df_qT <- df_q[,typeS=="type.T"]
rownames(df_qT) <- df_q[,2]
df_qP <- df_q[,typeS=="type.P"]
rownames(df_qP) <- df_q[,2]
df_qA <- df_q[,4:ncol(df_q)]
rownames(df_qA) <- df_q[,2]

# Differential expression
Sigs <- readxl::read_excel(paste(wdir, "Signatures_Table5S.xlsx", sep = "/"))
IFNa_genes<-intersect(df_q[,2],Sigs[[1]][-1])
TR1_genes<- intersect(df_q[,2],c(Sigs[[2]][2:17],"IL10RA"))
HLA_genes <- grep(pattern="HLA",df_q[,2],value = T)
MS_genes <- intersect(df_q[,2], c("EPCAM","CDH1","SOCS1","STAT1",
                                  "NLRC5","TGFB1","CD86","TAP1","TAP2",
                                  "CD3G","CD8A","EOMES","PDCD1","GZMK",
                                  "CD40","CD74","TREM2","LAG3","IL10RA",
                                  "HAVCR2","CTLA4","CASP3","CD69","ITGAE","ITGA1",
                                  "CD7","IL2","TNF","IL10"))

#Figures: P vs T
#Diff gene expression P vs T, selected genes
#Reshape wide to long matrix with expression values
df_qA2 <- reshape(df_qA, 
                  direction = "long",
                  v.names="value",
                  varying = names(df_qA),
                  times= names(df_qA),
                  timevar = "file",
                  idvar="genes",
                  ids = rownames(df_qA))
df_qA2$patient <- substr(df_qA2$file,7,11)
df_qA2$type <- substr(df_qA2$file,13,13)

## Plots
setwd(plots)

# Preparing_wNormByGene
GiveSigs <-  function(x1,x2){
    nor.genes <- x1[x2,]
    mean1 <- apply(nor.genes,1,geometric.mean)
    df1<- as.data.frame(t(nor.genes))
    mean2<- as.data.frame(t(mapply(`/`, df1,mean1)))
    result1<- apply(mean2,2,geometric.mean)
    names(result1) <- names(x1)
    return(result1)}
IFNa_Score.P <-GiveSigs(df_qP,IFNa_genes)
TR1_Score.P <- GiveSigs(df_qP,TR1_genes)
IFNa_Score.T <- GiveSigs(df_qT,IFNa_genes)
TR1_Score.T <- GiveSigs(df_qT,TR1_genes)

# Add IFNa and TR1 scores to df_qA
df_qT2 <- df_qT
colnames(df_qT2) <- substr(colnames(df_qT2), 7, 11)
df_qT3 <- data.frame(t(df_qT2))
df_qT3$type <- "T"
df_qT3$patient <- rownames(df_qT3)
df_qT3$TR1_Score <- TR1_Score.T
df_qT3$IFNa_Score <- IFNa_Score.T

df_qP2 <- df_qP
colnames(df_qP2) <- substr(colnames(df_qP2), 7, 11)
df_qP3 <- data.frame(t(df_qP2))
df_qP3$type <- "P"
df_qP3$patient <- rownames(df_qP3)
df_qP3$TR1_Score <- TR1_Score.P
df_qP3$IFNa_Score <- IFNa_Score.P
df_qAC <- rbind(df_qP3, df_qT3)

# Figure: Correlation TR1score and IFNascore
pdf(file = "Plot_TR1score_vs_IFNascore.pdf",width = 4, height = 4)
ggplot(df_qAC, aes(x=IFNa_Score, y=TR1_Score,color=type,fill=type))+
    geom_point()+
    geom_smooth(method=lm, se=T, fullrange=TRUE)+
    stat_cor(method = "spearman")+
    labs(title = "TR1score vs IFNascore")+
    xlab(label="IFNa signature score")+
    ylab(label="Tr1 signature score")
dev.off()

df.nano <- data.frame(rownames(df_qAC), df_qAC$IFNa_Score, df_qAC$TR1_Score, df_qAC$type)
openxlsx::write.xlsx(df.nano,"Plot_TR1score_vs_IFNascore.xlsx", rowNames = F)	


# Figure: Correlation CTLA4 and IFNascore
pdf(file = "Plot_CTLA4_vs_IFNascore.pdf",width = 4, height = 4)
ggplot(df_qAC, aes(x=IFNa_Score, y=CTLA4,color=type,fill=type))+
    geom_point()+
    geom_smooth(method=lm, se=T, fullrange=TRUE)+
    stat_cor(method = "spearman")+
    labs(title = "CTLA4 vs IFNascore")+
    xlab(label = "IFNa signature score")+
    ylab(label = "Log10 CTLA4 expression")
dev.off()

df.nano <- data.frame(rownames(df_qAC), df_qAC$IFNa_Score, df_qAC$CTLA4, df_qAC$type)
openxlsx::write.xlsx(df.nano, "Plot_CTLA4_vs_IFNascore.xlsx", rowNames = F)

# Figure: Correlation HLA-C and IFNascore
pdf(file = "Plot_HLA-C_vs_IFNascore.pdf",width = 4, height = 4)
ggplot(df_qAC, aes(x=IFNa_Score, y=HLA.C,color=type,fill=type))+
    geom_point()+
    geom_smooth(method=lm, se=T, fullrange=TRUE)+
    stat_cor(method = "spearman")+
    labs(title = "HLA-C vs IFNascore")+
    xlab(label = "IFNa signature score")+
    ylab(label = "Log10 HLA-C expression")
dev.off()

df.nano <- data.frame(rownames(df_qAC), df_qAC$IFNa_Score, df_qAC$HLA.C, df_qAC$type)
openxlsx::write.xlsx(df.nano, "Plot_HLA-C vs IFNascore.xlsx", rowNames = F)	

#Correlation all Genes vs IFNa
x1 <- df_qA
genesTotest <- subset(rownames(x1), !(rownames(x1)%in%c(IFNa_genes)))
tt1 <- x1[IFNa_genes,]
dfx <- data.frame()
for (i in genesTotest){
    AllCor <- apply(tt1,1,function(x) cor(x,as.numeric(x1[i,]),method="spearman"))
    TR1corr<- data.frame(value=AllCor,tested.Gene=rep(i,nrow(tt1)))
    TR1corr$Type <- "IFNa_genes"
    TR1corr$Median1 <- median(AllCor)
    dfx <- rbind(TR1corr,dfx)
}

dfx <- dfx[rev(order(dfx$Median1)),]
dfx2 <- dfx[1:(length(IFNa_genes)*100),]
dfx2 <- dfx2[order(dfx2$Median1, decreasing = T), ]
dfx2 <- head(dfx2, n = 20*16)
dfx2$tested.Gene <- factor(dfx2$tested.Gene, levels = unique(dfx2$tested.Gene))

pdf(file = "Plot_Top.corr.IFNaGenes.All.pdf",width = 12, height = 4)
ggplot(dfx2, aes(x = tested.Gene, y = value)) + 
    geom_boxplot() + geom_dotplot(binaxis='y', stackdir='center',
                                  position=position_dodge(0.75),dotsize = 0.5) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
    labs(title =  "Correlation against IFNa Genes")+
    xlab(label="Top genes")+
    ylab(label="rho coefficient")
dev.off()

openxlsx::write.xlsx(dfx2, "Plot_Top.corr.IFNaGenes.All.xlsx", rowNames = T)