#Open seurat and load the Seurat samples (can be downloaded from https://drive.google.com/drive/folders/1aMo8RgZt9bZy5tc94Kr1x-wuzUoRbOCz)
library(Seurat)
load("~/Desktop/CD4T6_seurat.rdata")
mycells = UpdateSeuratObject(object = mycells)
virus<-mycells@meta.data

#Also load the annotation files including sample source and TCR information (can be downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152522)
library(readr)
cd4t6_annotation <- read_table2("Desktop/cd4t6_annotation.txt")

TRA.aa.chains.tag<-cd4t6_annotation$TRA.aa.chains.tag
TRB.aa.chains.tag<-cd4t6_annotation$TRB.aa.chains.tag

mycells<-AddMetaData(mycells,metadata =TRA.aa.chains.tag,col.name = "TRA.aa.chains.tag")
mycells<-AddMetaData(mycells,metadata =TRB.aa.chains.tag,col.name = "TRB.aa.chains.tag")

#Extract cells that belongs to influenza vaccine cohort
mycells<-subset(x=mycells, subset=orig.virus2=="FLU")

#Extract singlets
mycells<-subset(x=mycells, subset=orig.MULTI_classification.global=="Singlet")

#Recluster the cells
mycells <- RunPCA(mycells, verbose = FALSE)
ElbowPlot(mycells, ndims =50)
mycells <- FindNeighbors(mycells, dims = 1:9, reduction="pca")
mycells <- RunUMAP(mycells, reduction = "pca", dims = 1:9)
mycells <- FindClusters(mycells, resolution = 0.7, reduction.type="pca")
FeaturePlot(mycells,"CXCR5",order = T,pt.size = 2)
FeaturePlot(mycells,"CD69",order = T,pt.size = 2,min.cutoff = 3) #There seems to have a more activated population, perhaps due to the antigen stimulation
DimPlot(mycells,label = T)

#Extract cTfh population
mycells <- subset(x = mycells, idents=c("2","3","5","4"))

#Calculate Tfh1/Tfh17 signature gene scores on each cell, the genes were showed in Figure S5A
Tfh1geneset <- read.csv("~/Desktop/Tfh1_up.csv")
mycells <- AddModuleScore(object = mycells,features = Tfh1geneset,name = 'Tfh1score', assay = "RNA")
Tfh17geneset <- read.csv("~/Desktop/Tfh17_up.csv")
mycells <- AddModuleScore(object = mycells,features = Tfh17geneset,name = 'Tfh17score')

#Exporting values used for plotting
aCDR3<-mycells@meta.data[["TRA.aa.chains.tag"]]
bCDR3<-mycells@meta.data[["TRB.aa.chains.tag"]]
donor<-mycells@meta.data[["orig.donor"]]
Tfh1score<-mycells@meta.data[["Tfh1score1"]]
Tfh17score<-mycells@meta.data[["Tfh17score1"]]
pre_post<-mycells@meta.data[["orig.virus"]]

#Generation of table used for plotting
TCRanalysis<-cbind(Tfh1score,Tfh17score,aCDR3,bCDR3,donor,pre_post)
write.csv(TCRanalysis,file="TCRanalysis.csv")

#Calculate the clonal abundance and to identify most abundant clones
TCRanalysis<-as.data.frame(TCRanalysis)
table1<-table(TCRanalysis$aCDR3,TCRanalysis$donor)
View(table1)
write.csv(table1,file="TCRanalysisabundance.csv")

#Then use excel to separate different blood donors; separate before/after vaccination; select the Tfh1/Tfh17 scores of the highly abundant clones and use Prism for plotting
