R-script
1.# Differential expression analysis of immune-related genes between 45,Xm patients and normal 46,XX females
read.csv('/Volumes/T7\ Touch/R/turner’s gene/mRNA in N and TSm groups.csv', header=T)->a
rownames(a) <- a[,1]
a <- a[,-1] 
Probe<-rownames(a) 
a<-data.frame(cbind(Probe,a))
read.csv('/Volumes/T7\ Touch/R/turner’s gene / gene\ annotation.csv',header=T)->b
c<-merge(a,b,by='Probe')#
head(c)
c<-c[,-1]
aggregate(c[,-27],by=list(c$Gene.Symbol),mean)->cc
names(cc)[1]='Gene.Symbol'
write.csv(cc,file = '/Volumes/T7\ Touch/R/turner’s gene/genesymbol origine.csv')
read.csv('/Volumes/T7\ Touch/R/turner’s gene /immuGeneList remained1.csv',header=T)->d
f<-merge(cc,d,by='Gene.Symbol')
write.csv(f,file='/Volumes/T7\ Touch/R/turner’s gene/sample(m-n)immnology gene.csv')
rownames(f)<-f[,1]
f<-f[,-1]
log_f=log2(f)
head(log_f)
write.csv(log_f,file="/Volumes/T7\ Touch/R/turner’s gene/（m-n）Standardized genetic data.csv")

group<-factor(rep(c('N','Tm'),c(10,16)))
design <- model.matrix(~0+group)
colnames(design) <- c('N','Tm')
rownames(design)<-colnames(log_f)
library(limma)
fit <- lmFit(log_f, design)
cont.matrix<-makeContrasts(Tm-N,levels=design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
final<-topTable(fit2, coef=1, number=dim(log_f)[1], adjust.method="BH", sort.by="B", resort.by="M")
write.csv(final,file='/Volumes/T7\ Touch/R/turner’s gene/TS(m-n)DEGmRNA-immugene_data-log.csv')


2.# Differential expression analysis of immune-related genes between 45,Xp patients and normal 46,XX females
read.csv('/Volumes/T7\ Touch/R/turner’s gene/mRNA\ in\ N\ and\ TSp\ groups.csv',header=T)->a1
rownames(a1) <- a1[,1]
a1 <- a1[,-1] 
Probe<-rownames(a1)  
a1<-data.frame(cbind(Probe,a1))
read.csv('/Volumes/T7\ Touch/R/turner’s gene/ gene\ annotation.csv',header=T)->b
c1<-merge(a1,b,by='Probe')
c1<-c1[,-1]
aggregate(c1[,-21],by=list(c1$Gene.Symbol),mean)->cc1
names(cc1)[1]='Gene.Symbol'
read.csv('/Volumes/T7\Touch/R/turner’s gene/immuGeneList remained1.csv',header=T)->d
f1<-merge(cc1,d,by='Gene.Symbol')
write.csv(f1,file='/Volumes/T7\ Touch/R/turner’s gene/sample(p-n)immnology gene.csv')
rownames(f1)<-f1[,1]
f1<-f1[,-1]
log_f1=log2(f1)
write.csv(log_f1,file="/Volumes/T7\ Touch/R/turner’s gene/（p-n）Standardized genetic data.csv")
group<-factor(rep(c('N','Tp'),c(10,10)))
design1 <- model.matrix(~0+group)
colnames(design1) <- c('N','Tp')
rownames(design1)<-colnames(log_f1)
library(limma)
fit1 <- lmFit(log_f1, design1)
cont.matrix<-makeContrasts(Tp-N,levels=design1)
fit3 <- contrasts.fit(fit1, cont.matrix)
fit3 <- eBayes(fit3)
final1<-topTable(fit3, coef=1, number=dim(log_f1)[1], adjust.method="BH", sort.by="B", resort.by="M")
write.csv(final1,file='/Volumes/T7\ Touch/R/turner’s gene/TS(p-n)DEGmRNA-immugene_data-log.csv')

3.# Volcano plot
3.1# Volcano map of differentially expressed Immune-Related genes (DEIRGs) between Xm TS patients and normal individuals.
library(ggplot2)
diff <- read.csv("/Volumes/T7\ Touch/R/turner’s gene/TS(m-n)DEGmRNA-immugene_data-log.csv",header = T)
logFC <-diff$logFC
adj <- diff$P.Value
data <- data.frame(logFC=logFC,padj=adj)
dim(data)
data$sig[(data$padj > 0.05|data$padj=="NA")|(data$logFC < 1)& data$logFC > -1] <- "no"
data$sig[data$padj <= 0.05 & data$logFC >= 1] <- "up"
data$sig[data$padj <= 0.05 & data$logFC <= -1] <- "down"
x_lim <- max(logFC,-logFC)

library(ggplot2)
library(RColorBrewer)
pdf(file = "immugene_volcano.pdf",width=8,height=8)
theme_set(theme_bw())
p <- ggplot(data,aes(logFC,-1*log10(padj),
                     color = sig))+geom_point() +  labs(x="log2(FoldChange)",y="-log10(P.Val)")
p <- p + scale_color_manual(values =c("#0072B5","grey","#BC3C28"))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)
p <- p +theme(panel.grid =element_blank())+
  theme(axis.line = element_line(size=0))+ylim(0,5)+xlim(-2,2)
p <- p  +guides(colour = "none")
p <- p +theme(axis.text=element_text(size=20),axis.title=element_text(size=20))
p
dev.off()
print(p)

3.2# Volcano map of differentially expressed Immune-Related genes (DEIRGs) between Xp TS patients and normal individuals.
library(ggplot2)
diff1 <- read.csv("/Volumes/T7\ Touch/R/turner’s gene/TS(p-n)DEGmRNA-immugene_data-log.csv",header = T)

logFC1<-diff1$logFC
adj1 <- diff1$P.Value
data1 <- data.frame(logFC=logFC1,padj=adj1)
dim(data1)
data1$sig[(data1$padj > 0.05|data1$padj=="NA")|(data1$logFC < 1)& data1$logFC > -1] <- "no"
data1$sig[data1$padj <= 0.05 & data1$logFC >= 1] <- "up"
data1$sig[data1$padj <= 0.05 & data1$logFC <= -1] <- "down"
write.csv(data1,file="/Volumes/T7\ Touch/R/turner’s gene/data1.csv")
x_lim <- max(logFC,-logFC)

library(ggplot2)
library(RColorBrewer)
pdf(file = "immugene_volcano.pdf",width=8,height=8)
theme_set(theme_bw())
p1<- ggplot(data1,aes(logFC,-1*log10(padj),
                      color = sig))+geom_point() +  labs(x="log2(FoldChange)",y="-log10(P.Val)")
p1 <- p1 + scale_color_manual(values =c("#0072B5","grey","#BC3C28"))+
  geom_hline(yintercept=-log10(0.05),linetype=4)+
  geom_vline(xintercept=c(-1,1),linetype=4)
p1 <- p1 +theme(panel.grid =element_blank())+
  theme(axis.line = element_line(size=0))+ylim(0,5)+xlim(-2,2)
p1 <- p1  +guides(colour = "none")
p1 <- p1 +theme(axis.text=element_text(size=20),axis.title=element_text(size=20))
p1
dev.off()
print(p1)

4. Heatmap
4.1 #Heatmap of upregulated and downregulated genes in Xm TS patients and normal individuals.
library(pheatmap)
aa<-read.csv("/Volumes/T7\ Touch/R/turner’s gene/TS(m-n)DEGmRNA-immugene_data-log-final.csv")
bb<-read.csv('/Volumes/T7\ Touch/R/turner’s gene/（m-n）Standardized genetic data.csv',header = T)
names(bb)[1]<-'Gene.Symbol'
cc<-merge(aa,bb,by='Gene.Symbol')
rownames(cc)<-cc[,1]
cc<-cc[,-1]
ann_colors = list(CFH = c(Mutated = "#CD2626", Wildtype = "#00008B"))
pheatmap(cc, cluster_row = T,cluster_col = FALSE,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50),
         gaps_col = 10,
         annotation_colors = ann_colors,border="white",frontsize_row=10,treeheight_row = 50)

4.2#Heatmap of upregulated and downregulated genes in Xp TS patients and normal individuals.
library(pheatmap)
aa1<-read.csv("/Volumes/T7\ Touch/R/turner’s gene/TS(p-n)DEGmRNA-immugene_data-log-final.csv",header = T)
bb1<-read.csv('/Volumes/T7\ Touch/R/turner’s gene/（p-n）Standardized genetic data.csv',header = T)
names(bb1)[1]<-'Gene.Symbol'
cc1<-merge(aa1,bb1,by='Gene.Symbol')
rownames(cc1)<-cc1[,1]
cc1<-cc1[,-1]
ann_colors = list(CFH = c(Mutated = "#CD2626", Wildtype = "#00008B"))
pheatmap(cc1, cluster_row = T,cluster_col = FALSE,
         color = colorRampPalette(c("navy", "white", "firebrick3"))(50),
         gaps_col = 10,
         annotation_colors = ann_colors,border="white",frontsize_row=10,treeheight_row = 50)

5. #Functional Enrichment Analysis of DEIRGs
5.1 # The top 10 terms in DEIRGs of Xm TS patients compared with normal individuals.
x<-read.csv('/Volumes/T7\ Touch/R/turner’s gene/GO.csv',stringsAsFactors=FALSE)
x<-x[x$PValue<0.05,]
x_go<-x[,1:5]
xbp<-x_go[grep("BP",x_go$Category),]
xcc<-x_go[grep("CC",x_go$Category),]
xmf<-x_go[grep("MF",x_go$Category),]
xkegg<-x_go[grep("KEGG",x_go$Category),]
xbp$Term=gsub(".*\\~","",xbp$Term)#Biological Process
xcc$Term=gsub(".*\\~","",xcc$Term)#cell Component
xmf$Term=gsub(".*\\~","",xmf$Term)#Molecular Function
xkegg$Term=gsub(".*\\:","",xkegg$Term)#KEGG pathway
library(ggplot2)
make_GO_bubble<-function(go_data,tern_name){
  GO_DATA<-go_data[order(go_data$Count,decreasing = T),]
  GO_DATA<-head(GO_DATA,10)
  p<-ggplot(GO_DATA,aes(X.,Term))
  bubble<-p+geom_point(aes(size=Count,color=-log10(PValue)))
  bubble=bubble+scale_color_gradient(low="green",high="red")
  pr=bubble+theme_test(base_size = 16,base_rect_size = 1)
  pr=pr+labs(x="Gene Ratio",y=tern_name,title = "")
  return(pr)
}
#BP
make_GO_bubble(xbp,tern_name="Biological Process")#HEIGHT 550
make_GO_bubble(xcc,tern_name="Cell Compont")
make_GO_bubble(xmf,tern_name="Molecular Function")
make_GO_bubble(xkegg,tern_name="KEGG pathway")

5.2  # The top 10 terms in DEIRGs of Xp TS patients compared with normal individuals.
x1<-read.csv('/Volumes/T7\Touch/R/turner’s gene/GO1.csv',stringsAsFactors=FALSE)
x1<-x1[x1$PValue<0.05,]
x1_go<-x1[,1:5]
xbp1<-x1_go[grep("BP",x1_go$Category),]
xcc1<-x1_go[grep("CC",x1_go$Category),]
xmf1<-x1_go[grep("MF",x1_go$Category),]
xkegg1<-x1_go[grep("KEGG",x1_go$Category),]
xbp1$Term=gsub(".*\\~","",xbp1$Term)#Biological Process
xcc1$Term=gsub(".*\\~","",xcc1$Term)#cell Component
xmf1$Term=gsub(".*\\~","",xmf1$Term)#Molecular Function
xkegg1$Term=gsub(".*\\:","",xkegg1$Term)#KEGG pathway

library(ggplot2)
make_GO_bubble1<-function(go1_data,term_name){
  GO_DATA1<-go1_data[order(go1_data$Count,decreasing = T),]
  GO_DATA1<-head(GO_DATA1,10)
  p<-ggplot(GO_DATA1,aes(X.,Term))
  bubble<-p+geom_point(aes(size=Count,color=-log10(PValue)))
  bubble=bubble+scale_color_gradient(low="green",high="red")
  pr=bubble+theme_test(base_size = 16,base_rect_size = 1)
  pr=pr+labs(x="Gene Ratio",y=term_name,title = "")
  return(pr)
}
#BP
make_GO_bubble1(xbp1,term_name="Biological Process")#HEIGHT 550
make_GO_bubble1(xcc1,term_name="Cell Compont")
make_GO_bubble1(xmf1,term_name="Molecular Function")
make_GO_bubble1(xkegg1,term_name="KEGG pathway")  

6. # Immune Cells Assessment
6.1# cibersort analysis
a<-read.table('/Volumes/T7\Touch/R/turner’s gene /cibersort/DATA1.txt',header=T,sep="\t")
rownames(a)<-a[,1]
setwd("/Users/spine/Desktop/cibersort")
source('Cibersort.R')
result1 <- CIBERSORT('LM22.txt','DATA1.txt',  perm = 1000, QN = T)
6.2 # The relative percentage of immune cells between monosomy X TS patients
and normal individuals.
pkgs <- c("matrixStats", "pheatmap", "RColorBrewer", "tidyverse", "cowplot","ggpubr","bslib","ggthemes")
lapply(pkgs, library, character.only = T)
install.packages("tidyverse")
library('tidyverse')
# Read in results ---------------------------------------------------------

#m-n Cibersort result visualization
cibersort_raw <- read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/CIBERSORT_result.csv",row.names = 1,header = T)
library(dplyr)
library(tidyr)
res<-cibersort_raw [,colSums(cibersort_raw )>0]
typeof(res1)
is.matrix(res1)
a=t(res)
mycol<-ggplot2::alpha(rainbow(ncol(res)),0.7)
par(bty="o",mgp=c(2.5,0.3,0),mar=c(2.1,4.1,2.1,10.1),tcl=-.25,las=1,xpd=T)
res1<-as.matrix(t(res))
barplot(res1,
        border=NA,
        #names.arg=rep("",nrow(res)),
        names.arg=colnames(res1),
        yaxt="n",
        ylab = "Relative percentage",
        col=mycol)
axis(side = 2,at=c(0,0.2,0.4,0.6,0.8,1),
     labels = c("0%","20%","40%","60%","80%","100%"))
legend(par("usr")[2]-1,
       par("usr")[4],legend = colnames(res),
       xpd=T,
       fill=mycol,
       cex=0.7,
       border = NA,
       y.intersp = 1,
       x.intersp = 0.2,
       bty="n")
#dev.off()

#p-n Cibersort result visualization
cibersort_raw1 <- read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/CIBERSORT_result_2.csv",row.names = 1,header = T)
library(dplyr)
library(tidyr)
res1<-cibersort_raw1[,colSums(cibersort_raw1)>0]
mycol<-ggplot2::alpha(rainbow(ncol(res1)),0.7)
par(bty="o",mgp=c(2.5,0.3,0),mar=c(2.1,4.1,2.1,10.1),tcl=-.25,las=1,xpd=T)
barplot(as.matrix(t(res1)),
        border=NA,
        names.arg=rep("",nrow(res1)),
        yaxt="n",
        ylab = "Relative percentage",
        
        col=mycol)
axis(side = 2,at=c(0,0.2,0.4,0.6,0.8,1),
     labels = c("0%","20%","40%","60%","80%","100%"))
legend(par("usr")[2]-1,
       par("usr")[4],legend = colnames(res1),
       xpd=T,
       fill=mycol,
       cex=0.7,
       border = NA,
       y.intersp = 1,
       x.intersp = 0.2,
       bty="n")
#dev.off()

6.3 # Boxplot of Cibersort results
library(tidyverse)
#The differences in immune cell proportions between Xm TS patients (red color) and the normal individuals (blue color).

a<- read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/cibersort\ boxplot/CIBERSORT_result.csv",row.names = 1,header = T)
b<-read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/cibersort\ boxplot/group.csv",row.names = 1,header = T)
a<-a [,colSums(a )>0] identical(rownames(a),rownames(b))
b$group<-as.factor(b$group)
class(b$group)
a$group<-b$group
a<-a %>%
  rownames_to_column("sample")
#install.packages("ggsci")
#install.packages("ggpubr")
library(ggsci)
library(tidyr)
library(ggpubr)
b<-gather(a,key=CIBERSORT,value = Proportion,-c(group,sample))
ggboxplot(b,x="CIBERSORT",y="Proportion",
          fill = "group",palette = "lancet",outlier.colour = NA)+
  scale_y_continuous(breaks = seq(0, 0.6, 0.1))+    
  stat_compare_means(aes(group=group),
                     method = "wilcox.test",
                     label = "p.signif",
                     label.y = 0.6,
                     symnum.args=list(cutpoints=c(0,0.001,0.01,0.05,1),
                                      symbols=c("***","**","*","")))+
  theme(text=element_text(size=10),
        axis.text.x = element_text(angle=45,hjust=1))

library(tidyverse)
#The differences in immune cell proportions between Xp TS patients (yellow color) and the normal individuals (blue color). 
a1<- read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/cibersort\ boxplot/CIBERSORT_result_2.csv",row.names = 1,header = T)
b1<-read.csv("/Volumes/T7\ Touch/R/turner’s gene/cibersort/cibersort\ boxplot/group_2.csv",row.names = 1,header = T)
a1<-a1 [,colSums(a1 )>0
        identical(rownames(a1),rownames(b1))
        b1$group<-as.factor(b1$group)
        class(b1$group)
        a1$group<-b1$group
        a1<-a1 %>%
          rownames_to_column("sample")
        #install.packages("ggsci")
        #install.packages("ggpubr")
        library(ggsci)
        library(tidyr)
        library(ggpubr)
        b1<-gather(a1,key=CIBERSORT,value = Proportion,-c(group,sample))
        ggboxplot(b1,x="CIBERSORT",y="Proportion",
                  fill = "group",palette = "jco",outlier.colour = NA)+
          scale_y_continuous(breaks = seq(0, 0.6, 0.1))+
          stat_compare_means(aes(group=group),
                             method = "wilcox.test",
                             label = "p.signif",
                             symnum.args=list(cutpoints=c(0,0.001,0.01,0.05,1),
                                              symbols=c("***","**","*","")))+
          theme(text=element_text(size=10),
                axis.text.x = element_text(angle=45,hjust=1))
        
        
        
        
        
        
        
        
        
        
        
        
        