前言

  上一次《TCR-seq数据分析一》分享了如何初步处理TCR-seq的数据，这一次接着来说如何进一步处理数据。上一章介绍了免疫受体的V(D)J结构，当发生免疫反应时，VDJ区域的基因重排会形成很多不同的免疫受体，即免疫组库。TCR-seq数据分析的主要目的就是统计各区域基因的出现频率，即geneUsage。这个时候就轮到今天的主角上场了——immunarch是一个R包，可以用来对很多软件的TCR-seq数据如mixcr、10X等做后续的数据分析。该R包含有丰富的处理函数以及多样性的数据展示类型，用起来也是相当方便，那话不多说来看一下如何使用该R包。

代码

install.packages("immunarch")
library(immunarch)
data(immdata)
gene_stats()
##     alias                 species ighd ighj ighv igij igkj igkv iglj iglv traj
## 1      bt               BosTaurus   21    4   25    0    1    6    5   26   46
## 2      cd      CamelusDromedarius    0    0    0    0    0    0    0    0    0
## 3     clf    CanisLupusFamiliaris    0    0    0    0    0    0    0    0    0
## 4      dr              DanioRerio    7    7    0    3    0    0    0    0    0
## 5      hs             HomoSapiens   30   13  248    0    5   64    7   69   57
## 6  macmul           MacacaMulatta   24    7   19    0    4   83    5    0    0
## 7     mmc    MusMusculusCastaneus    0    0    0    0    0    4    0    0    0
## 8     mmd   MusMusculusDomesticus    0    0    0    0    0    2    0    0    0
## 9  musmus             MusMusculus   32    8  225    0    8  109    3    5   42
## 10     oa OrnithorhynchusAnatinus    3   10    0    0    0    0    0    0    0
## 11     oc    OryctolagusCuniculus   10   11   39    0    8   26    2   20    0
## 12     om      OncorhynchusMykiss    9    7    6    0    0    0    0    0    0
## 13     rn        RattusNorvegicus   30    4  113    0    6  132    2    8    0
## 14   smth   MusMusculusMolossinus    0    0    0    0    0    1    0    0    0
## 15   smth     MusMusculusMusculus    0    0    0    0    0    1    0    0    0
## 16   smth              MusSpretus    0    0    0    0    0    2    0    2    0
## 17     ss               SusScrofa    5    5   15    0    8   19    4   14    0
##    trav trbd trbj trbv trdd trdj trdv trgj trgv
## 1     0    0    0    0    5    3    0    6   15
## 2     0    0    0    0    0    0    7    2    2
## 3     0    2    8   19    0    0    0    7    8
## 4     0    0    0    0    0    0    0    0    0
## 5    60    3   14   64    3    4    6    4   10
## 6     0    2   15   58    0    0    0    0    0
## 7     0    0    0    0    0    0    0    0    0
## 8     0    0    0    0    0    0    0    0    0
## 9   145    2   14   23    2    3    7    0   11
## 10    0    0    0    0    0    0    0    0    0
## 11    0    0    0    0    0    0    0    0    0
## 12    0    1    9    0    0    0    0    0    0
## 13    0    0    0    0    0    0    0    0    0
## 14    0    0    0    0    0    0    0    0    0
## 15    0    0    0    0    0    0    0    0    0
## 16    0    0    0    0    0    0    0    0    0
## 17    0    0    0    0    0    0    0    0    0

imm_gu <- geneUsage(immdata$data[c(1, 2)], "hs.trbv", .norm = T)
vis(imm_gu)

vis(imm_gu, .grid = T)

vis(imm_gu, .by = "Status", .meta = immdata$meta, .plot = "box")

imm_gu_cor <- geneUsageAnalysis(imm_gu, .method = "cor", .verbose = F)
vis(imm_gu_cor, .title = "Gene usage correlation", .leg.title = "Cor", .text.size = 1.5)

imm_cl_pca <- geneUsageAnalysis(imm_gu, "js+pca+kmeans", .verbose = F)
vis(imm_cl_pca, .plot = "clust")

imm_cl_pca2 <- geneUsageAnalysis(imm_gu, "js+pca+kmeans", .k = 3, .verbose = F)
vis(imm_cl_pca2)

p1 <- vis(spectratype(immdata$data[[1]], .quant = "id", .col = "nt"))
p2 <- vis(spectratype(immdata$data[[1]], .quant = "count", .col = "aa+v"))
p1 + p2

  现在是不是觉得immunarch包很nice，同时值得提醒的是immunarch提供的可视化是基于ggplot2的vis函数来展示结果，这样对可视化做起来修改那也是相当的方便。该包还有更多的分析功能，我这里就不一一展示了，留个大家自己去探索吧。

最后

  TCR-seq关键要理解这些数据结构以及其背后能够反应的生物学意义，这一点我也在继续学习中。。。