张泽民老师课题组的新成果，值得细细学习。以下是一刷的学习笔记：

---

summary

---

• In melanoma patients, CD8+ tumor-infiltrating T cells exhibit a linear and continuous progression from predysfunctional cell state to dysfunction , but in lung cancer patients, there are two pre-exhaustion states that could develop to exhaustion

• cancer types such as multiple myeloma do not show notable exhausted T cell populations

• CD8+ cells in the tumor were featured by the emergence of exhausted T cells, whereas among CD4+ T cell populations, the most abundant population was the TNFRSF9+ Treg cell, which showed significantly lower frequencies in both blood and normal tissues.

  
  

• To distinguish the T cells reacting to tumors from bystander T cells, we jointly analyzed the features of tissue distribution, transcriptional phenotypes, proliferation, and clonal expansion. The characteristics of proliferation and clonal expansion of tumor-enriched T cells have been viewed as evidence of their tumor reactivity
  

  

• Taken together, potentially tumor-reactive T cells emerged—TFH/TH1 cells, TNFRSF9+ Treg cells, CD8+ISG+ T cells, and four Tex cell populations, representing a local antitumor immune response—whereas the expanded CD8+ Temra cells might also harbor tumor-specific TCRs, which is consistent with the notion of a systemic immune response.

• CD8+ T细胞的分化轨迹

1. on a global scale

   
   
   

   整体来说Tn有两条分化路径，一条是Tn--Temra，一条是Tn--Tex

2. for the Tex branch

   
   
   

   在Tn--Temra的路径中，又有两条路径，一条是经过Tem，一条是经过Trm

先总体看，再看看这样的分化轨迹是否是存在在各个样本或者各个癌种中

• 耗竭亚群的标志基因

  
  

• 关于CD4+ T细胞

1. In the CD4+ compartment, the major potentially tumor-reactive metaclusters were IFNG+ TFH/TH1 and TNFRSF9+ Treg cells.

2. The global diffusion map and RNA velocity analysis revealed that CD4+ T cells could develop from naive T cells to Temra cells, TFH/TH1 or
   TNFRSF9+ Treg seperately.

   
   
3. 都是先做总体的轨迹推断，然后再在每一个方向做进一步的轨迹推断
4. The RUNX3 exhibited elevated expression at a point at which a high density of TFH/TH1 cells emerged, which is consistene with a previous report that RUNX3 regulated the cytotoxic phenotype in the CD4+ cytotoxic T cells. 另外还有一个很重要的调节因子是TP73
5. For Treg cells, a trajectory from TNFRSF9- Treg cell to TNFRSF9+ Treg cell emergef, indicating a gradual transistion from the resting state to a activated state.
6. High inteferon-stimulated genes(ISG) as a feature of intermediate state during CD4+ T cell activation.
7. Various conventional CD4+ T cell populations had conversion relationships with Treg cells, but such conversion patterns were diverse and vaired among cancer types.

• The T cell composition in the tumor could be affected not only by the number of potential neoantigens, reflected by TMB but also by specific somatic mutations of cancer cells.

• 他们根据免疫细胞亚群的不同，把癌症分为了不同的五种：

  
  
  

  C1: high frequency of terminal Tex and TNFRSF9+ Treg，如食管癌，鼻咽癌
  C2: high frequency of terminal Tex
  C3-C8: 终末Tex少，high frequency of CD8+ZNF683+CXCR6+Trm
  C7：dominated by naive T cell
  C8：enriched naive T cell
  C6: enriched Temra cell
  C4: enriched Tc17 or TH17 cell 基底细胞癌，肝细胞癌
  C5: with a low frequency of TNFRSF9+ Treg

Although each immune type included mixed cancer types, certain cancer types exhibited clear preferences.

---

方法学笔记

• 这篇文章很厉害的一点就是能够把不同来源、不同测序方法、不同时间的单细胞样本很好得整合在一起。对于有这样需求的实验，可以仔细阅读这篇文章的methods。

• Then 40 CD3+CD8+CD4- and CD3+CD8-CD4+ T cells were isolated based on average expression of CD3 genes (CD3D, CD3G), CD8 genes (CD8A, CD8B), and CD4.

• To reduce noise in the data, we performed two-run clustering on each dataset. Examining the result from the first run clustering, we identified contamination clusters and clusters that arose from unwanted factors. The second run clustering without such noise prepared the basis for data integration across multiple datasets. 这种方法都可以学习一下，感觉很适合平时在分析时使用。

• For all data sets, the resolution parameter of clustering was set to 2.0, which could produce sufficient fine clustering according to our experiments. 我的习惯是先低精度聚类，注释大亚群，再根据需要进一步聚类。他们是直高精度聚类，再合并。

• marker gene的选择：All other parameters were kept as the default values. After clustering, limma (version 3.42.2) was used for detecting differentially expressed genes. For each cluster, cells were compared with cells from all other clusters. Genes with log2 fold change larger than a specified threshold (1.00 and 0.25 for SmartSeq2-based datasets and droplet-based datasets 20 respectively)** and false discovery rate (FDR) < 0.01 were defined as the marker genes of the cluster.
  10x: LogFC>1
  smartseq2:LogFC>0.25

• disassociation induced genes (DIG): 消化实体肿瘤时可能出现的干扰基因，比如HSP90

• 如何评价数据整合的好坏？——We used the local inverse Simpson’s Index (LISI) to quantitatively evaluate the integration. The LISI defines the effective number of datasets in a neighborhood of a cell. While in bad integration, neighborhoods represented by only a single dataset will give LISI value 1; higher LISI value indicates neighborhoods presented by more datasets and thus lower batch effect.

• TCR diversity的评价指标：Shannon 40 equitability index (normalized Shannon diversity index) EH.

• TCR assembly and clonotype identification pipeline

1. Only those assembled chains that were highly confident, of full-length, productive, and with a valid cell barcode and an unambiguous chain type (alpha/beta/gamma/delta) assignment were qualified and kept.
2. If a cell had at least one pair of qualified alpha and beta chains, this cell was annotated as an alpha/beta T cell. If a cell had at least one pair of qualified gamma and delta chains, this cell was annotated as a gamma/delta T cell. If a cell had both qualified alpha/beta pair(s) and gamma/delta pair(s), this cell was annotated using the TCR pair with the highest UMI counts.
3. If a cell had more than one qualified chain of the same chain type, only the two chains with the highest UMI counts were kept, and the one with the higher UMI count was determined as the dominant chain.
4. For each patient, cells with identical dominant alpha/beta chains (or gamma/delta chains) were considered to originate from the same clonal expansion, therefore they were assigned the same clonotype ID.

• In a few cases, CD4+ and CD8+ T cells shared the same clonotype. To reduce noise and confusion, for each clonotype with both CD4+ T cells and CD8+ T cells, we calculated the cell number ratio between CD4+ T cells and CD8+ T cells, then we removed cells using the following criteria: 1) if the ratio > 5, the CD8+ T cells were removed; 2) if the ratio < 0.2, the CD4+ T cells were removed; 3) otherwise, both CD4+ T cells and CD8+ T cells were removed.（这一点我还没考虑过，感觉可以注意一下）

• proliferating score was calculated as the average z-score scaled expression values of the common proliferation markers (ZWINT, E2F1, FEN1, FOXM1, H2AFZ, HMGB2, MCM2, MCM3, MCM4, MCM5, MCM6, MKI67, MYBL2, PCNA, PLK1, CCND1, AURKA, BUB1, TOP2A, TYMS, DEK, CCNB1, and CCNE1)

• Identification of individual potentially tumor-reactive T cells
  Step one: With "REACTOME_TCR_SIGNALING" from MSigDB and common proliferation markers as genesets, GSEA implemented in R package clusterProfiler was applied to the z-score expression profiles of mini-clusters to obtain the normalized enrichment scores (NES) and p-values. Mini-clusters with NES > 0 and p-value < 0.05 were considered as having high TCR signaling or 35 proliferation, and tumor infiltrated T cells in those mini-clusters were used as “baits”.
  Step two: The expanded clonotypes to which the baits belong were considered as potentially tumor-reactive clonotypes.
  Step three:All cells of those potentially tumor-reactive clonotypes were considered as pTRTs, including both cells that exhibited signature of TCR signaling or proliferation and cells that did not.

Although whether those cells were truly tumor-reactive required further experimental validation, it was helpful for the characterization of meta-clusters to analyze the distribution patterns of those pTRTs.

• Gene Set Enrichment Analysis (GSEA) was performed to investigate the pathway activities of meta-clusters. For each meta-cluster, the vector of combined effect size was extracted and then used as input for the software GSEA (version 4.0.3) in the pre-ranked mode. The annotation information of gene sets was downloaded from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp ). The scoring scheme was set to "classic", and the permutation number was set to 1000. Over-representation analysis was performed to annotate given gene lists of interest using the R package clusterProfiler (version 3.15.2).

-关于RNA velocity
In the global analyses of CD8+ T cells and CD4+ T cells, and the state transition analyses of CD4+Tfh related populations and CD4+ Treg populations, to reduce computational consumption, we applied mini-cluster level analyses, implemented by averaging the normalized counts per mini-cluster.
It should be noted that gene selection impacted the result of the RNA velocity and the explanation. In the global analyses of CD8+ T cells and CD4+ T cells, we set the highly variable genes to the informative genes which were involved in the meta-cluster analysis (see section "Data integration and meta-cluster identification"). Those genes should be more powerful in distinguishing different cell states than the genes identified in an unsupervised manner such as the default implementation of scvelo. In the analysis of exhaustion, we set the highly variable genes 30 to the informative genes among the meta-clusters of the two major exhaustion paths. We defined the “exhaustion program” as the top 50 signature genes of terminal Tex. To visualize the transition potential of the “exhaustion program”, those genes were embedded into the UMAP. In other cases, the RNA velocities of the highly variable genes were embedded into the UMAP or diffusion map, and should be explained as the transition potential of “the overall transcriptomic state”.

To investigate the exhaustion process, we first performed diffusion map analysis using the informative genes. The result clearly showed that naive T cells were in one branch, and Temra and Tex were in two different branches , indicating that T cells could differentiate into two major different fates
利用informative genes来做diffusion map和RNA velocity!!!

• 文章中用过的轨迹推断的方法：

1. diffusion map
2. monocle3
3. RNA velocity
4. PAGA
5. TCR tracking
6. Slingshot

-他们用了Nichnet来进行分析