一个R包玩转单细胞免疫组库分析，还能与Seurat无缝对接
單細胞免疫組庫數據分析
NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。
在介紹clonotypr(https://github.com/mattfemia/clonotypr)包之前，我們需要先去解決幾個問題：
（1）為什么要進行單細胞免疫組庫的分析？
應用方向一：探索腫瘤免疫微環境，輔助免疫治療。
每個人都擁有一個自己的適應性免疫組庫，TCR和BCR通過基因重組和體細胞突變取得多樣性，使得我們身體可以識別和抵御各種內部和外部的入侵者。而腫瘤的發生往往躲避了人體T淋巴細胞而產生、增殖和轉移。
使用10X Genomics ChromiumTM Single Cell Immune Profiling Solution可以捕捉腫瘤發生時的免疫微環境變化，尋找免疫治療的靶點，從而輔助免疫治療更好地抗擊腫瘤。
應用方向二：探索自身免疫性疾病和炎癥性疾病發生機制，輔助疫苗的研究
自身免疫性疾病發生起始和發展的中心環節被認為是抗原特異性T細胞激活導致的，使用10X Genomics ChromiumTM Single Cell Immune Profiling Solution，可以解析自身免疫性疾病的發病機制，從而為疾病的診療提供依據。
應用方向三：移植和免疫重建
器官或者骨髓移植時，經常會誘發宿主的排斥反應，從而發生慢性移植抗宿主病。同種異體反應隨機分布在整個T細胞組庫的交叉反應，因此延遲T細胞恢復和限制的T細胞受體多樣性與異體移植后感染和疾病復發風險增加相關。
而我比較注意的是在疫苗接種前后BCR/TCR CDR3免疫組庫的分析，最近medRxiv上發表的有關新冠的文獻Immune Cell Profiling of COVID-19 Patients in the recovery stage by Single-cell sequencing中對不同BCR/TCR的VDJ重排進行分析，揭示針對新冠特異的克隆擴增。
（2）免疫組庫主要包括哪幾個方面？
T淋巴細胞（T cell）和B淋巴細胞（B cell）主要負責適應性免疫應答，其抗原識別主要依賴于T細胞受體（T cell recptor, TCR）和B細胞受體（B cell recptor, BCR），這兩類細胞表面分子的共同特點是其多樣性，可以識別多種多樣的抗原分子。BCR的輕鏈和TCRβ鏈由V、D、J、C四個基因片段組成，BCR的重鏈和TCRα鏈由V、J、C三個基因片段組成，這些基因片段在遺傳過程中發生重組、重排，組合成不同的形式，保證了受體多樣性。其中變化最大的就是CDR3區。
（3）10X Genomics在進行VDJ測序時為什么從5’端開始？
我們都知道3’端首先捕獲的是polyA尾，然后其實就是C段恒定區，由于長度限制難以對VDJ段進行測序。（但其實我們在進行10X Genomics3’測序時，也可以發現一些IGH類的基因表達。。。）
（4）10× Genomics進行cellranger后的輸出形式是什么樣的？
Outputs: - Run summary HTML: /home/jdoe/runs/sample345/outs/web_summary.html - Run summary CSV: /home/jdoe/runs/sample345/outs/metrics_summary.csv - All-contig FASTA: /home/jdoe/runs/sample345/outs/all_contig.fasta - All-contig FASTA index: /home/jdoe/runs/sample345/outs/all_contig.fasta.fai - All-contig FASTQ: /home/jdoe/runs/sample345/outs/all_contig.fastq - Read-contig alignments: /home/jdoe/runs/sample345/outs/all_contig.bam - Read-contig alignment index: /home/jdoe/runs/sample345/outs/all_contig.bam.bai - All contig annotations (JSON): /home/jdoe/runs/sample345/outs/all_contig_annotations.json - All contig annotations (BED): /home/jdoe/runs/sample345/outs/all_contig_annotations.bed - All contig annotations (CSV): /home/jdoe/runs/sample345/outs/all_contig_annotations.csv - Filtered contig sequences FASTA: /home/jdoe/runs/sample345/outs/filtered_contig.fasta - Filtered contig sequences FASTQ: /home/jdoe/runs/sample345/outs/filtered_contig.fastq - Filtered contigs (CSV): /home/jdoe/runs/sample345/outs/filtered_contig_annotations.csv - Clonotype consensus FASTA: /home/jdoe/runs/sample345/outs/consensus.fasta - Clonotype consensus FASTA index: /home/jdoe/runs/sample345/outs/consensus.fasta.fai - Clonotype consensus FASTQ: /home/jdoe/runs/sample345/outs/consensus.fastq - Concatenated reference sequences: /home/jdoe/runs/sample345/outs/concat_ref.fasta - Concatenated reference index: /home/jdoe/runs/sample345/outs/concat_ref.fasta.fai - Contig-consensus alignments: /home/jdoe/runs/sample345/outs/consensus.bam - Contig-consensus alignment index: /home/jdoe/runs/sample345/outs/consensus.bam.bai - Contig-reference alignments: /home/jdoe/runs/sample345/outs/concat_ref.bam - Contig-reference alignment index: /home/jdoe/runs/sample345/outs/concat_ref.bam.bai - Clonotype consensus annotations (JSON): /home/jdoe/runs/sample345/outs/consensus_annotations.json - Clonotype consensus annotations (CSV): /home/jdoe/runs/sample345/outs/consensus_annotations.csv - Clonotype info: /home/jdoe/runs/sample345/outs/clonotypes.csv - Barcodes that are declared to be targeted cells: /home/jdoe/runs/sample345/out/cell_barcodes.json - Loupe V(D)J Browser file: /home/jdoe/runs/sample345/outs/vloupe.vloupeclonotypr
clonotypr是一個R包，旨在使用Seurat（單細胞分析Seurat使用相關的10個問題答疑精選！）作為“骨干”，將V（D）J-seq/單細胞免疫譜數據與scRNA-seq數據的分析連結在一起。除了將克隆型數據附加到對象之外，clonotypr還提供了分析CDR3aa長度、化學組成和頻率，以及將單鏈克隆型重新分配給其'productive'-relatives的策略。
安裝
if (!require("devtools")) {install.packages("devtools") }install.packages("devtools") devtools::install_github("mattfemia/clonotypr")library(clonotypr)library(kableExtra) # 數據表格格式化clonotypr建立在Seurat對象基礎上，其中seurat示例對象如pbmc_vdj?和pbmc_null都已內置于包中。
pbmc_null?->不包含克隆型數據。
pbmc_vdj?->在meta.data中包含克隆型數據。
用于構建“pbmc_null”和“pbmc_vdj”的數據均來自人PBMC，由10X Genomics提供。可在https://support.10xgenomics.com/single-cell-vdj/datasets/3.1.0/vdj_v1_hs_pbmc3找到有關樣本、scRNA-seq和VDJ庫的更多信息。
處理數據
AddClonotypes（）和LoadClonotypes（）是clonotypr的典型函數：
LoadClonotypes（）：
創建一個新的數據框，不需要Seurat對象。
AddClonotypes（）：
將克隆型數據追加到Seurat對象中。
后面的函數將以Seurat為對象進行輸入，并且都需要CellRanger輸出中的兩個文件：
filtered_contig_annotations.csv
clonotypes.csv
可以使用兩個參數將文件進行調用：
上傳帶有directory參數的目錄
使用clonotypeFile（clonotypes.csv）和filteredContigFile（filtered_contigs_annotations.csv）上傳文件
我們從Seurat對象開始：
head(pbmc_null@meta.data)object@meta.data?——?增加clonotype data之前:
head(pbmc_null@meta.data) %>% # <-- Can also be accessed with pbmc_vdj[[]]kable() %>%kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 13) %>%row_spec(0, bold = T, color = "black", background = "white") %>%scroll_box(width = "100%", height = "250px") #生成數據表AddClonotypes（）?——?加上clonotype data：
clonoFile <- system.file("extdata", "vdj_v1_hs_pbmc3_t_clonotypes.csv", package = "clonotypr") fcaFile <- system.file("extdata", "vdj_v1_hs_pbmc3_t_filtered_contig_annotations.csv", package = "clonotypr") clono <- AddClonotypes(pbmc_null,clonotypeFile = clonoFile,filteredContigFile = fcaFile) # Individual file uploadhead(clono@meta.data) #可以看到已將兩個csv表加入clono中；orig.ident nCount_RNA nFeature_RNA barcode clonotype_id v_gene d_gene AAACCTGTCAACGGGA SeuratProject 5433 1470 AAACCTGTCAACGGGA clonotype104 TRAV9-2 None AAACGGGGTTAAAGAC SeuratProject 3226 1147 AAACGGGGTTAAAGAC clonotype32 TRBV9 None AAAGATGAGTCACGCC SeuratProject 7432 1958 AAAGATGAGTCACGCC clonotype113 TRBV7-6 None AAAGATGCAGCTGCAC SeuratProject 7105 2014 AAAGATGCAGCTGCAC clonotype114 TRAV8-3 None AAAGATGCAGTAGAGC SeuratProject 6313 1504 AAAGATGCAGTAGAGC clonotype115 TRAV19 None AAAGATGCATCCGTGG SeuratProject 3541 1364 AAAGATGCATCCGTGG clonotype116 TRAV1-2 Nonej_gene c_gene cdr3s_aa AAACCTGTCAACGGGA TRAJ38 TRAC TRA:CALSGFHNAGNNRKLIW;TRB:CASSLILRGEQFF AAACGGGGTTAAAGAC TRBJ1-1 TRBC1 TRB:CASKGETNTEAFF AAAGATGAGTCACGCC TRBJ1-1 TRBC1 TRA:CALSEARAAGNKLTF;TRA:CAVRDFIGFGNVLHC;TRB:CASSVGQITEAFF AAAGATGCAGCTGCAC TRAJ42 TRAC TRA:CAAMDSNYQLIW;TRA:CAVDPDYGGSQGNLIF;TRB:CASSLDYEQYF;TRB:CASSLSSGANVLTF AAAGATGCAGTAGAGC TRAJ30 TRAC TRA:CALSEASRDDKIIF;TRB:CASSPDWRGDTDTQYF AAAGATGCATCCGTGG TRAJ33 TRAC TRA:CYSMDSNYQLIW;TRB:CAISSGRAADIQYFcdr3s_nt AAACCTGTCAACGGGA TRA:TGTGCTCTGAGTGGTTTCCATAATGCTGGCAACAACCGTAAGCTGATTTGG;TRB:TGTGCCAGCAGTTTAATACTCAGGGGTGAGCAGTTCTTC AAACGGGGTTAAAGAC TRB:TGTGCCAGCAAGGGGGAAACGAACACTGAAGCTTTCTTT AAAGATGAGTCACGCC TRA:TGTGCTCTGAGTGAGGCAAGGGCTGCAGGCAACAAGCTAACTTTT;TRA:TGTGCTGTGAGAGACTTTATAGGCTTTGGGAATGTGCTGCATTGC;TRB:TGTGCCAGCAGCGTCGGACAGATCACTGAAGCTTTCTTT AAAGATGCAGCTGCAC TRA:TGTGCTGCCATGGATAGCAACTATCAGTTAATCTGG;TRA:TGTGCTGTGGACCCCGATTATGGAGGAAGCCAAGGAAATCTCATCTTT;TRB:TGCGCCAGCAGCTTGGATTACGAGCAGTACTTC;TRB:TGTGCCAGCAGTCTCTCTTCTGGGGCCAACGTCCTGACTTTC AAAGATGCAGTAGAGC TRA:TGTGCTCTGAGTGAGGCTAGCAGAGATGACAAGATCATCTTT;TRB:TGTGCCAGCAGTCCCGACTGGCGGGGGGACACAGATACGCAGTATTTT AAAGATGCATCCGTGG TRA:TGCTATTCCATGGATAGCAACTATCAGTTAATCTGG;TRB:TGTGCCATCTCTAGCGGGCGAGCCGCAGACATTCAGTACTTChead(clono@meta.data) %>%kable() %>%kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 13) %>%row_spec(0, bold = T, color = "black", background = "white") %>%scroll_box(width = "100%", height = "250px")克隆型鏈
DescribeClonotypes可以快速檢查樣本中克隆型鏈類型的分布。
DescribeClonotypes(clono)頻率分析
與filtered_contig_annotations.csv文件中的數據相似，clonotypr計算相對于總樣本群體的每種克隆型的原始頻率和百分比。
freqs <- MakeFrequencyDF(clono) head(freqs)freqs[freqs$clonotypeVector[35:45],] nrow(freqs)freqs[freqs$clonotypeVector[35:45],] %>%kable() %>%kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 13) %>%row_spec(0, bold = T, color = "black", background = "white") %>%scroll_box(width = "100%", height = "250px") nrow(freqs)單鏈克隆
clonotypr還可以處理遇到僅包含α或β鏈的克隆型的可能情況。那么我們就該想，不是說好的α和β鏈嘛，為什么只有一條，于是這是10X官方的解釋：
在克隆型頻率要求樣本內具有高精度的情況下，JoinSingleChains可以解決此問題。在此，首先分離各個CDR3alpha和CDR3beta鏈。然后將那些單鏈克隆型的頻率添加到productive克隆型的頻率中，最后僅返回productive克隆型和頻率。
prod_freqs <- JoinSingleChains(clono)prod_freqs[prod_freqs$clonotypeVector[35:45],] nrow(prod_freqs)prod_freqs <- JoinSingleChains(clono) prod_freqs[prod_freqs$clonotypeVector[35:45],] %>%kable() %>%kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"), font_size = 13) %>%row_spec(0, bold = T, color = "black", background = "white") %>%scroll_box(width = "100%", height = "250px") nrow(prod_freqs)CDR3 Lengths
CDR3鏈長度可以使用clonotypr進行分析：
AnalyzeCDR3Length（）->表格輸出
PlotLengths（）->可視化輸出
長度表
生成將細胞barcodes與頻率數據合并的表格：
barcodefreqs <- GetBarcodeFreqs(clono)lengths <- AnalyzeCDR3Lengths(barcodefreqs, productive=TRUE)head(lengths)Plotting Lengths
length_dist <- PlotLengths(lengths) length_dist[1] # CDR3alpha lengths length_dist[2] # CDR3beta lengths length_dist[3] # Combined lengthsCDR3氨基酸組成
要分析樣品的氨基酸組成，我們可以按CDR3長度將其分解，然后關注每個氨基酸位置以研究其多樣性。
clonotypr提供了兩種用于CDR3組成分析的方法：
AnalyzeCDR3aa（）->表格輸出
PlotCDR3aa（）->可視化輸出
CDR3位置矩陣
“AnalyzeCDR3aa”創建一個位置矩陣，并返回一個數據框列表（而不是矩陣）。這些提供了不同長度的CDR3序列之間氨基酸組成的定量測量。在此處詳細了解詳情：https://en.wikipedia.org/wiki/Position_weight_matrix
三個“輸出”選項：
位置頻率矩陣（PFM）
位置概率矩陣（PPM）
位置權重矩陣（PWM）
繪制CDR3組成
我們可以使用PlotCDR3aa將位置權重矩陣（PWM）作為seqlogos進行分析（R語言 - 繪制seq logo圖）。此輸出可用于可視化整個樣本中所有長度序列的序列概率：
comp <- PlotCDR3aa(lengths) comp[1] # <- CDR3alpha chains comp[2] # <- CDR3beta chains繪制CDR3位置矩陣
盡管位置計算已經集成到了函數“AnalyzeCDR3aa”中，但clonotypr允許以下方式來簡化每種位置矩陣的計算：
BuildPFM（）->位置頻率矩陣列表
BuildPPM（）->位置概率矩陣列表
BuildPWM（）->位置權重矩陣列表
這些函數和“AnalyzeCDR3aa”輸出的表格都可以繪制位置矩陣的熱圖，并且每種類型的矩陣都提供有單獨的函數：
PlotPFM（）-> PFM熱圖列表
PlotPPM（）-> PPM熱圖列表
PlotPWM（）-> PWM熱圖列表
校對：生信寶典
參考
10X Genomics單細胞測序全新升級，一次實驗同時玩轉單細胞轉錄組和單細胞免疫組庫：
http://www.bioon.com.cn/news/showarticle.asp?newsid=74423
clonotypr：
https://github.com/mattfemia/clonotypr
你可能還想看
Seurat亮點之細胞周期評分和回歸
單細胞分析Seurat使用相關的10個問題答疑精選！
綜述：變溫動物的適應性免疫
TISIDB：這個科研工具能助你順利發高分腫瘤免疫文章
往期精品(點擊圖片直達文字對應教程)
后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
總結
以上是生活随笔為你收集整理的一个R包玩转单细胞免疫组库分析，还能与Seurat无缝对接的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 分享清华大学鲁志教授实验室生物信息学教程
- 下一篇： 【NGS接龙】薛宇：漫谈生物信息圈儿的那