【打算一步步按照文章来做里面的分析，做到哪写到哪，记录下~】
说先要清楚几个定义：
diversity：文中是用的observed diversity，以及inverse Simpson's index(哈哈哈，去年弄过新冠宏转录组的分析，这一套我熟悉~)
clonality：1-Pielou‘s evenness (就是1-Pielou均匀度，J = Shannon.Wiener / log(S)，S=specnumber(data))
开始一步步分析
首先对31个样本的asites、blood、tumor做γδT检测，看看各个部位里面γδT多寡（用的流式细胞仪？我不懂这个仪器能干嘛，以为只是一个数细胞个数的玩意，发现功能不止于此），然后计算下median（excel直接可以算，简单方便，上一个R的计算）：
> d<- read.table("/project/baowenjuan/Tumor_FFPE/project/draw/EOC/data.1A",header=T,sep="\t")
> d1<-d[,2:4]
> rownames(d1)<-d[,1]
> median(d1[,1])
[1] 2.12
> median(d1[,2])
[1] 3.41
> median(d1[,3])
[1] 1.83
Fig.1A
然后对8个病人的3个部位的样本的TCR的CDR3区域进行测序，我理解的就是类似细菌里面16S的V4区域测序，反正就是引物扩增子，然后上NGS，计算各个类型，算几个多样性的指标。它这里已经是转化为氨基酸序列的统计的。那么就有一个原始矩阵：每一行是各个motif的集合，每一列是一个样本（总共3*8列）。可是文章没有提供，但是呢，根据文章提供的S3表格，我把这个矩阵整理出来了（真是个聪明的小天才~）
文章的S3表格数据，两两对比的
于是把相同样本的都放一起，根据两两关系，把原始矩阵整理出来：
IM11为例，3组样本的原始矩阵。最后的others是总数据量-已列出的和，具体motif我也不知道
先测试下IM11这个样本，看对不对：启动R，开始计算：
d<-read.table("TCR.IM11",header=T)
head(d)
sequence IM11A IM11B IM11T
CAAWDYPDPTKLF 16 2 0
CALWAQGFWELGKKIKVF 2339 2410 612
CALWDGAEELGKKIKVF 24 53 0
CALWEAQELGKKIKVF 12 51 2
CALWEEIELGKKIKVF 54 119 8
CALWEEIVLGKKIKVF 71 149 0
dataR = rrarefy(d1,min(rowSums(d[,2:4]))) /对数据做一个均一化矫正/
diversity(dataR[,1],index = "inv")
发现我们算出的inverse辛普森值跟文章给的恰好相反，其实问题就出在others上。看数据，asites和blood有比较明显的单体motif expansion，而tumor里面要分散得多，这样的数据，想想都知道多样性要高，但因为我整理出的矩阵的motif有限，那些低频率的都统一归为了others,这对tumor样本影响很大，所以这个矩阵算的不做数，还是要原始的矩阵。
ins Simpson
文章给了Observed diversity数据，那我们就用这个画个图吧(excel画的，嘻嘻嘻)：
Fig.1C
Observed Diversity
Inversed Simpson
Clonality
接着画一个Krona图，文章里用的 IM11号样本的数据，但是我左算右算，得不出图片里的频率，哎，头大anyway，还是画个图玩玩
/biostack/bioconda/pkgs/krona-2.7.1-pl526_3/bin/ktImportText krona.txt
krona.txt
IM11B Krona图
文章对31个样本的 γδ T细胞的频率做了一个boxplot图，来，画一画~
d<-read.table("pct.txt",header=True)
library(ggplot2)
library(dplyr)
library(showtext)
showtext_auto(enable = TRUE)
font.add('SimSun', '/usr/share/fonts/zh_CN/simsun.ttf')
ggplot(d,aes(x=group,y=PCT))+geom_boxplot(outlier.shape = NA,col=c("lightpink","darkorange","cyan"))+geom_jitter(width = 0.1)+theme(axis.text.x=element_text(angle=45,hjust = 1,colour="black",family="Times",size=10),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"))+xlab("")+ylab("%γδT cells of total T cells")
d1<-as.matrix(read.table("pct.txt1",header=T))
dimnames(d1)<-list(1:31,c("Blood","Ascites","Tumor"))
friedman.test(d1)
Firedman结果
我画出的样子
论文里的样子
再画一个各个样本种类中各TCR种类的累积频率图（画图渣渣一个，临时百度~）
累积分布图
library(Rmisc)
d<-read.table("TCR",header=T)
d$TCR<-as.factor(d$TCR)
head(d)
PCT SampleType TCR
31.85 Blood 1
51.03 Blood 2
61.86 Blood 3
64.00 Blood 4
66.09 Blood 5
67.71 Blood 6
d2<- summarySE(d,measurevar = "PCT",groupvars =c("TCR","SampleType"))
head(d2)
TCR SampleType N PCT sd se ci
1 Ascites 8 24.46000 12.444657 4.399851 10.403994
1 Blood 8 15.62875 8.209323 2.902434 6.863166
1 Tumour 8 14.36125 12.152796 4.296662 10.159991
2 Ascites 8 35.66375 16.114574 5.697362 13.472121
2 Blood 8 26.32500 13.006231 4.598397 10.873482
2 Tumour 8 21.22125 13.134351 4.643694 10.980592
ggplot(d2,aes(x=TCR,y=PCT,group=SampleType,color=SampleType))+geom_line()+geom_point()+ geom_errorbar(aes(ymin=PCT - 0.5*sd, ymax=PCT + 0.5*sd),width=.2, position=position_dodge(0.05))+labs( x="Distinct TRG CDR3(top10)", y = "% of repertoire (accumulated)")+geom_jitter(width = 0.1)+theme(axis.text.x=element_text(angle=5,hjust = 1,colour="black",family="Times",size=10),panel.grid.major = element_blank(),panel.grid.minor = element_blank(),panel.background = element_blank(),axis.line = element_line(colour = "black"))
画图好累啊，要调整，这一步暂时跳过~
文章对各个频率分层，然后统计“频率密度”，发现tumor中主要集中在低频，ascites在>5%的区段密度有一个峰，意思就是说特定一些clones在ascites这有expansion。文章接着对TCR氨基酸的长度分布，以及种类做了一个统计，发现blood里主要是14个aa，ascites里14，tumor里12。然后对这3个最高丰度的长度的motif类型做个统计，在WebLogo（WebLogo 3 - About (threeplusone.com)）上画了一个图（Fig.1H）。
接着对gene segment usage(不懂这个概念)进行统计，就是看看不同segment在不同样本类型中的频率。（流式细胞仪会产生这个数据？就是这个仪器可以分清楚细胞类型，是Vδ1+还是Vδ2+，以及segment也分得清？TRGV系列~对这个实在是没一个理性、感性认识）结论是tumor里以Vδ1+为主，blood里以Vδ2+为主。
文章分别对Vδ1+和Vδ2+中CD69+CD103+的T细胞进行统计（数据源：S4：Fig.2D），发现来源于tumor的要多。也发现CD69+CD103+的Vδ1+频率与clonality负相关，但Vδ2+没有。因为clonality是TCR测序数据出来的，反应的是TCR种类的多寡以及是否富集，这里的意思说表型的富集并不是由于TCR种类的富集导致的。也对 Vδ1 + CD27 lo/neg T 表型（该表型跟转向记忆细胞有关，且是TCR依赖的）进行统计，发现tumor里较少，且与clonality不相关，说明tumor浸润γδ T细胞的反应是非TCR依赖的。但ascites里的 Vδ1 + CD27 lo/neg T 与clonality正相关，说明ascites的里是TCR导致的适应性扩增。【这里就是对表型跟TCR型做一个相关性研究，看表型是不是TCR型导致的结果，大致这个意思吧~tumor里更多是非TCR依赖，即innate；ascites里是TCR-driven的，即adaptive】
然后文章用adaptive-（anti-γδ TCR 、anti-CD3、IL-15） 和innate-like stimuli（anti-MICA/B），发现tumor和ascites的差异：tumor对anti-MICA/B刺激又反应，较ascites，tumor的MICA + 要多，然后anti-MICA/B在αβT里没反应，所以说名对anti-MICA/B反应的是tumor里的γδ T。
Clonality Vδ1
0.418959677 8.11
0.403540671 2.77
0.367390573 10
0.280934036 43.2
0.292384207 42.1
0.329681098 17.2
0.328588516 49.1
0.200026765 64.5
> cor.test(d$Clonality,d$Vδ,method="spearman")
Spearman's rank correlation rho
data: d$Clonality and d$Vδ
S = 160, p-value = 0.004563
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-0.9047619
看文章的数据应该是选的spearman method。
excel就可以画
Ascites 和blood-derived γδT cells的TCR序列一致性高，但tumor的要差异大点
然后文章对各个样本（这里换到微生物里来说，就是菌落），做一个聚类，距离用的F pairwise similarity（relative overlap frequencies的几何平均值，我感觉这里是我们平常算的β多样性指标，但这里具体怎么算的没明白~）。然后相同样本类型之间的相似性指标，用的是Morisita-Horn。那么这里就有一个矩阵（行数：各个TCR种类，列数：样本数，即3x8，针对每一个样本类型，假如tumor，两两相比就有C2/8中，即8x7/2=28，文章直接把数据给出来了，然后画了个图）：
!Morisita-Horn 数据](https://img.haomeiwen.com/i21340631/0ded42aa59982f47.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
Morisita-Horn咋计算的呢，举个例子（数据随便找的一对）：具体见：R: The Morisita index and Horn-Morisita index (r-project.org)
> library("abdiv")
> d<-read.table("Morisita.txt",header=T)
> head(d)
sampleA_frequency sampleB_frequency
1 0.182778141 0.114534161
2 0.030481167 0.007453416
3 0.006422817 0.000000000
4 0.011430438 0.002360248
5 0.005443066 0.000124224
6 0.005225343 0.000745342
> horn_morisita(d$sampleA_frequency, d$sampleB_frequency)
[1] 0.04119227
总得来说，这里就是先把24个样本按照样本类型来分（3类），计算类内人与人的相似性（8个人之间的, 有8x7/2=28次比较）；又对类内人与人之间共享的clonotype占比分析（8个人，3种比较，共计24个点，图Fig3.C）；又分别两两种样本类型比对：举个例子：tumor里median（8个人的中位数）>1%的TCR和blood里的比，也有blood里median（8个人的中位数）>1%的TCR和tumor比，两个的区别就是候选TCR是按照谁>1%的来选。所以理论上总共有3*2种比较，但因有的缺少，所以文章只有4对比较的。发现非Vγ9的类型: blood<ascites<tumor。
tumor和ascites derived γδ T不产生IL-17A
文章给的数据：
Fig5B
但文章说“but the production was significantly higher among ascites-derived αβ T cells compared to their γδ T cell counterparts (0.2% versus 0.1%,P = 0.042; Fig. 5B)”里的不一样，这里是IL-17A+的T细胞比例啊。不明白用的什么数据~
后面是数据是对一个well里细胞数来做均一化（γδ：6500；αβ：100000）来计算。为什么原始excel里除的6000？有点疑惑。
PMA/I刺激，IFN-γ，TNF-α，MIP-1β+的γδ T细胞多且在各组间差异大
画生存分析图：
Fig.5F
> head(d)
Time status group
1 15.7 1 Low TNF α
2 31.6 1 Low TNF α
3 31.3 0 Low TNF α
4 14.8 1 Low TNF α
5 8.9 1 Low TNF α
6 1.7 1 Low TNF α
library(survminer)
library(survival)
ggsurvplot(fit,palette = c("#E7B800", "#2E9FDF"),title="TNF-α response by γδ T cells(PMA/I)",ylab="Overall survival",xlab ="Time(months)",pval = TRUE,legend.title = "",legend = c(0.8,0.15),legend.labs = c("Low TNF α","High TNF α"))
微小残留与cytokine 的相关性分析，残留越大，细胞因子越少，即γδ T对抑制肿瘤生长、进展。γδ T（TNF-α+）与CD39+的γδ T细胞呈负相关。说明γδ T对CD39更敏感。
TCGA数据库 表达量与OS
TRDV1 TRDV2：γδ T
TRDV2: Vδ2+T
TRBC2: αβT
用基因表达量分别代表对应细胞的含量。对CD39（ENTPD1产物）表达高低做生存分析。
感觉
CD39是不是像PDL1一样的东西啊，当然不是作用机制上，就是治疗方向上类似。
文章包括基本统计、NGS数据分析（扩增子测序数据分析、多样性分析[α、β多样性]）、相关性分析、TCGA数据库整理、基因表达分析、生存分析。对生信的要求是基本功要扎实，有一定的综合经验要求，绘图能力要加强。
网友评论