免疫组库数据分析||immunarch教程：Kmer 与 Motif 分析

immunarch — Fast and Seamless Exploration of Single-cell and Bulk T-cell/Antibody Immune Repertoires in R

今天，我们继续我们的免疫组库数据分析的Demos，这一次我们来谈谈Kmer 与 Motif 分析。像我这样刚入门免疫组库的人首先会问什么是Kmer 与 Motif ？

不巧的是，18年初当时你们家运来小哥哥还在做宏基因组了解过Meta小课题：K-mer.

k-mer 是指将reads分成包含k个碱基的字符串，一般长短为m的reads可以分成m-k+1个k-mers.

Kmer 反应序列的特性，在我们免疫组库数据中cdr3和nt序列，自然可以计算他们的Kmer 了。

在生物学领域里，也经常会提到到motif这个词。模体（motif）是蛋白质分子中具有特定空间构象和特定功能的结构成分。其中一类就是具有特殊功能的超二级结构。一个模体总有其特征性的氨基酸序列，并发挥特殊的功能。一般而言，常见的motion有以下几种形式：有规则地聚集在一起形成全由α-螺旋、全由β-片层或α-螺旋与β-片层混合、具体说，形成相对稳定的αα、βββ、βαβ、β2α和αTα等超二级结构，又称模体(motif)或模序。这百科。

Kmer 分析

在immunarch计算Kmer 很简单。

library(immunarch); data(immdata)       # Load the package and the test dataset
?getKmers

kmers <- getKmers(immdata$data[[1]], 3)
kmers

# A tibble: 4,405 x 2
   Kmer  Count
   <chr> <int>
 1 AAA       6
 2 AAD       5
 3 AAE       9
 4 AAF       2
 5 AAG      37
 6 AAI       2
 7 AAK       5
 8 AAL       4
 9 AAM       1
10 AAN      13
# ... with 4,395 more rows

如果有可能计算一批免疫序列中出现的kmers。为了做到这一点，您只需要向函数提供一个免疫程序列表。NA表示样品中没有这样的kmer，具体由列名决定。

kmers <- getKmers(immdata$data, 5)
kmers

# A tibble: 172,926 x 13
   Kmer  `A2-i129` `A2-i131` `A2-i133` `A2-i132` `A4-i191`
   <chr>     <int>     <int>     <int>     <int>     <int>
 1 AAAAG        NA        NA        NA        NA        NA
 2 AAAAL        NA        NA        NA        NA        NA
 3 AAAAW        NA        NA        NA        NA        NA
 4 AAADE        NA        NA        NA        NA        NA
 5 AAADN        NA         1        NA        NA        NA
 6 AAADT        NA        NA        NA        NA        NA
 7 AAAEA        NA        NA        NA        NA        NA
 8 AAAEN        NA        NA        NA         1        NA
 9 AAAET        NA        NA        NA        NA        NA
10 AAAFE        NA        NA        NA         1        NA
# ... with 172,916 more rows, and 7 more variables:
#   `A4-i192` <int>, MS1 <int>, MS2 <int>, MS3 <int>,
#   MS4 <int>, MS5 <int>, MS6 <int>

注意，默认情况下，getKmers()会在计算kmer统计数据之前过滤掉所有非编码序列。通过设置coding参数为FALSE，你可以同时使用编码和非编码序列:

kmers <- getKmers(immdata$data[[1]], 3, .coding = F)
kmers

# A tibble: 4,645 x 2
   Kmer  Count
   <chr> <int>
 1 **G       1
 2 *~G       4
 3 *~L       1
 4 *AA       1
 5 *D~       1
 6 *EE       1
 7 *G~       1
 8 *GG       1
 9 *GP       1
10 *HL       1
# ... with 4,635 more rows

一样的啦，可视化很简单

kmers <- getKmers(immdata$data, 5)
vis(kmers)

kmers的vis()函数有许多参数来操作图形。首先，.head参数指定了要可视化的最丰富的kmers数量。

p1 <- vis(kmers, .head = 5)
p2 <- vis(kmers, .head = 10)
p3 <- vis(kmers, .head = 30)

(p1 + p2) / p3

选项“stack”将所有的条堆叠在一起，这样你就可以看到kmers的完整分布。选项“fill”堆叠所有条在彼此之上也一样，但常态化它在这样的方式，所以你看到每kmer计数的分布，也就是说，你可以清楚地看到哪个曲目比其他特定kmer有更多的kmer。选项“dodge”将不同样本的kmer条分组，以便您可以清楚地看到，哪些样本总体上有更多的kmer出现。

如果您的kmer计数的分布对于某些曲目非常不平衡，则需要另外的参数是.log。它允许使用y轴的对数变换所以你可以看到kmer计数的数量级上的差异。

p1 <- vis(kmers, .head = 10, .position = "stack")
p2 <- vis(kmers, .head = 10, .position = "stack", .log = T)
p3 <- vis(kmers, .head = 10, .position = "fill",.log = T)
p4 <- vis(kmers, .head = 10, .position = "dodge", .log = T)

(p1 + p2) /(p3+ p4 )

Motif 分析

immunarch使用常见的方法进行序列基序分析（ Motif 分析），并使用不同类型的矩阵来表示序列基序:

位置频率矩阵(PFM)——每个氨基酸在每个位置出现的矩阵;
位置概率矩阵(PPM) -一个矩阵，每个氨基酸在每个位置出现的概率;
位置权重矩阵(PWM) -一个矩阵与PPM元素的对数概率;
PWM中具有元件自信息的矩阵。

为了计算和可视化序列基序，首先你需要计算一个输入免疫序列的kmer统计，然后应用kmer_profile()函数来计算序列基序矩阵:

kmers <- getKmers(immdata$data[[1]], 5)
kmer_profile(kmers)

kmer_profile(kmers)
   [,1]  [,2]  [,3] [,4] [,5]
A  8955  8976  3791 3246 3159
C  6469    26    27   18   10
D  2129  2131  2131 2130 2085
E  3315  5868  5872 5870 5720
F   688   691   691 2600 9029
G  9194  9603  9603 9566 9326
H   405   405   405  680  659
I   750  1020  1107  896  852
K   721  1099  1100 1100 1044
L  3073  3078  4108 4097 4033
M   241   245   246  246  226
N  2421  2423  2424 2419 2355
P  2285  2564  2564 2559 2509
Q  2459  2466  7087 7086 7051
R  3487  3502  3504 3496 2917
S 14021 14029 13082 8128 3462
T  3408  5669  5671 6024 5797
V  1703  1927  1927 1561 1508
W   485   486   487  450  439
Y  2060  2061  2442 6097 6088
attr(,"class")
[1] "immunr_kmer_profile_pfm" "matrix"

目前我们不支持对多个样本进行序列基序分析，但我们正在进行研究。为了计算和可视化所有样本的序列基序矩阵，您需要逐个处理它们，这可以通过for循环或lapply()函数轻松完成。

.method = "freq" - position frequency matrix (PFM);
.method = "prob" - position probability matrix (PPM);
.method = "wei" - position weight matrix (PWM);
.method = "self" - self-information matrix.

如：

kmer_profile(kmers, "freq")
   [,1]  [,2]  [,3] [,4] [,5]
A  8955  8976  3791 3246 3159
C  6469    26    27   18   10
D  2129  2131  2131 2130 2085
E  3315  5868  5872 5870 5720
F   688   691   691 2600 9029
G  9194  9603  9603 9566 9326
H   405   405   405  680  659
I   750  1020  1107  896  852
K   721  1099  1100 1100 1044
L  3073  3078  4108 4097 4033
M   241   245   246  246  226
N  2421  2423  2424 2419 2355
P  2285  2564  2564 2559 2509
Q  2459  2466  7087 7086 7051
R  3487  3502  3504 3496 2917
S 14021 14029 13082 8128 3462
T  3408  5669  5671 6024 5797
V  1703  1927  1927 1561 1508
W   485   486   487  450  439
Y  2060  2061  2442 6097 6088
attr(,"class")
[1] "immunr_kmer_profile_pfm" "matrix"

序列基序矩阵的可视化由vis()完成。有两种图形可供选择—— sequence logo和“text logo”。参数.plot规定了plot的类型:.plot =“seq”用于 sequence logo图形，.plot =“text”(默认情况下)用于“text logo”图形。

kp <- kmer_profile(kmers, "self")
p1 <- vis(kp)
p2 <- vis(kp, .plot = "seq")

p1 + p2

序列标志是一个图形表示的氨基酸或核酸多序列。每个标志由成堆的符号,一个堆栈中的每个位置序列。堆栈的总高度表示该位置的序列守恒，而堆栈内符号的高度表示每个氨基酸或核酸在该位置的相对频率。一般来说，序列标志比一致序列提供更丰富、更精确的描述，例如对结合位点的描述。

本节课我们体验了Kmer 与 Motif 分析，这两个都是基于序列来分析的，进一步细化了CDR3的aa和nt的多样性。之前的分析都是把他们当作一个整体看多个样本之间的关系。那么，这样做的意义什么呢？这就偏向于结构了，尽管序列结构只是初级结构。

http://weblogo.threeplusone.com/create.cgi
https://immunarch.com/articles/web_only/v9_kmers.html