hello，大家好，今天我們來分享一個10X單細(xì)胞和VDJ聯(lián)合分析的方法，PlatypusDB，關(guān)于轉(zhuǎn)錄組和VDJ聯(lián)合分析的方法，之前也分享了不少了，這篇文章參考于Platypus: an open-access software for integrating lymphocyte single-cell immune repertoires with transcriptomes，2021年6月發(fā)表于NAR Genomics and Bioinformatics，方法值得一看，今天我們就來分享這個方法。首先強(qiáng)調(diào)一句，一定要關(guān)注VDJ氨基酸的序列。

當(dāng)然，這個方法最新才升級，由于最近 Cellranger（第 5 版和第 6 版）中默認(rèn)克隆分型策略的變化，目前重建 Platypus 的 v3 以圍繞 VDJ_GEX_matrix 函數(shù)。 該函數(shù)集成了VDJ repertoire和轉(zhuǎn)錄組信息，并將作為包未來迭代中所有次要函數(shù)的輸入。 這樣做的好處是在每個細(xì)胞水平上擁有所有的repertoire和轉(zhuǎn)錄組信息。

Platypus 是一個包，旨在促進(jìn)單細(xì)胞免疫組庫測序?qū)嶒灥姆治觥?該軟件包可用于單獨分析基因表達(dá) (GEX) 或免疫受體庫 (VDJ) 測序數(shù)據(jù)，此外還可以集成兩個數(shù)據(jù)集以將表型特征與repertoire分析相結(jié)合。 該軟件包旨在主要分析 10x 基因組細(xì)胞管理器的輸出（來自 GEX 和 VDJ 計數(shù)的輸出，用于豐富的免疫受體庫）。 假設(shè)正確添加了輸入列，這些函數(shù)可以與其他基于條形碼的 scSeq 技術(shù)一起使用。 基因表達(dá)分析在很大程度上依賴于 Seurat，這是一種常用的單細(xì)胞測序 (scSeq) R 包。

圖片.png

### Removing any previous versions of the package
#First can ensure that there is no previous version installed locally
#detach("package:Platypus", unload=TRUE)
#remove.packages("Platypus")

### Dependencies 
#install.packages("stringr")

### Downloading and installing Platypus

# First we need to download the most recent version from the master branch at https://github.com/alexyermanos/Platypus we can install the package using the following command. 
# WARNING: This needs to be replaced with your own directory where the downloaded package is found

# For MacOS users it may look like this
#install.packages("~/Downloads/Platypus_3.0.tar.gz", repos = NULL, type="source")

# For windows it will likely look something like this. 
# WARNING: You will need to replace 'YourPCName' with your user name for the windows account in the directory. 
# install.packages("C:\Users\YourPCName\Downloads\Platypus_3.0.tar.gz", repos = NULL, type="source")

# Now we can load the installed package into the R environment. In case of problems with installing other R packages that are used in Platypus, please see the README file at the https://github.com/alexyermanos/Platypus, where we outline how to install the other R packages for both Windows and MacOS.
#library(Platypus)

# The individual R functions can additionally be found on the github in the Functions branch. Within this branch, there is a folder "R" which contains the individual functions. This can similarly be downloaded and loaded into the R environment in case not all functions are desired. These functions are actively updated and may include more features than the in original tar.gz file. 

Extracting and integrating repertoire data with VDJ_GEX_matrix (Platypus v3),其實這一步就是在提取免疫細(xì)胞的轉(zhuǎn)錄組和VDJ數(shù)據(jù)。

### Downloading the test data for VDJ_GEX_matrix 
# While the Platypus manuscript uses the COVID-19 data, the vignette for Platypus v3 will use the data from B cells in the aged CNS, which can be found at the following link https://polybox.ethz.ch/index.php/s/fxQJ3NrRSwiPSSo This small dataset contains VDJ (separate libraries for B) and GEX libraries from the central nervous system of two murine samples. More information can be found https://doi.org/10.1098/rspb.2020.2793

# After downloading the zip file named "Platypus_CNS_data.zip", please unzip the file and find the path to the newly formed folder. Typically this will be in the Downloads folder, so the code below should work on MacOS. For windows please uncomment the code and change the user name to match your PC.

VDJ.out.directory.list <- list() ### Set directory to the outs folder of cellranger vdj
VDJ.out.directory.list[[1]] <- c("~/Downloads/Platypus_CNS_data/VDJ_S1/")
VDJ.out.directory.list[[2]] <- c("~/Downloads/Platypus_CNS_data/VDJ_S2/")

GEX.out.directory.list <- list() ### Set directory to the outs folder of cellranger count
GEX.out.directory.list[[1]] <- c("~/Downloads/Platypus_CNS_data/GEX_S1/")
GEX.out.directory.list[[2]] <- c("~/Downloads/Platypus_CNS_data/GEX_S2/")

# For windows: 
#directory_to_covid_patients_gex[[1]] <- c("C:\Users\YourPCName\Downloads\PlatypusTestData\Patient1_GEX")
#directory_to_covid_patients_gex[[2]] <- c("C:\Users\YourPCName\Downloads\PlatypusTestData\Patient2_GEX")

# We will call the output vgm (short for Vdj_Gex_Matrix) - this object can be supplied as input to downstream functions in v3 of the package.
vgm <- Platypus::VDJ_GEX_matrix(VDJ.out.directory.list = VDJ.out.directory.list,
                               GEX.out.directory.list = GEX.out.directory.list,
                               GEX.integrate = T,
                               VDJ.combine = T,
                               integrate.GEX.to.VDJ = T,
                               integrate.VDJ.to.GEX = T,
                               exclude.GEX.not.in.VDJ = F,
                               filter.overlapping.barcodes.GEX = F,
                               filter.overlapping.barcodes.VDJ = F,
                               get.VDJ.stats = T,
                               parallel.processing = "none",
                               subsample.barcodes = F,
                               trim.and.align = F,
                               group.id = c(1,2))
## The output will be a list - 
# vgm[[1]] corresponds to the VDJ master dataframe
# vgm[[2]] corresponds to the GEX in the form of a seurat object
# vgm[[3]] corresponds to the output of VDJ_stats subfunction - which provides information about the number of chains, sequencing reads, etc
# vgm[[4]] holds the input parameters 
# vgm[[5]] holds the sessionInfo output at the time of function call

接下來首先是單獨分析，Gene expression analysis (Platypus v2)

cellranger 的 count 函數(shù)的輸出以表達(dá)矩陣、條形碼和基因標(biāo)識符的形式返回基因表達(dá)信息。函數(shù)automated_GEX 使我們能夠自動對來自 cellranger 的基因表達(dá)文庫進(jìn)行轉(zhuǎn)錄分析。輸入目錄應(yīng)設(shè)置為包含這三個文件的目錄：barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz。如果需要分析多個不同的轉(zhuǎn)錄組（例如，在單獨的 UMAP 空間上），這些目錄應(yīng)進(jìn)入單獨的列表對象。例如，如果有 20 個repertoires并且想要分別分析它們（20 個單獨的 UMAP，20 個單獨的repertoires），那么輸入目錄的每個外部列表元素將包含各個repertoires/GEX 文件的目錄。此函數(shù)的輸出是一個 Seurat 對象，類似于 Seurat 注釋中演示的標(biāo)準(zhǔn)管道。這涉及縮放、歸一化、聚類和執(zhí)行降維（默認(rèn)為 tSNE 和 UMAP）。在這個階段，沒有納入repertoires特征，需要使用 Platypus 中的后續(xù)功能進(jìn)行集成。

當(dāng)然， Platypus集成了Seurat的功能。(我們簡單看看)

### Removing any previous versions of the package
#First, ensure there is no previous version installed locally
#detach("package:Platypus", unload=TRUE)
#remove.packages("Platypus")

### Downloading and installing Platypus

# Download most recent version from master branch at https://github.com/alexyermanos/Platypus We can install the package using the following command. 
# WARNING: This needs to be replaced with your own directory where the downloaded package is found

# For MacOS users it may look like this:
#install.packages("~/Downloads/Platypus_2.0.5.tar.gz", repos = NULL, type="source")

# For windows it will likely look something like this: 
# WARNING: You will need to replace 'YourPCName' with your user name for the windows account in the directory. 
# install.packages("C:\Users\YourPCName\Downloads\Platypus_2.0.6.tar.gz", repos = NULL, type="source")

# Now we can load the installed package into the R environment. In case of problems with installing other R packages that are used in Platypus, please see the README file at the https://github.com/alexyermanos/Platypus, where we outline how to install the other R packages for both Windows and MacOS.
library(Platypus)

# Individual R functions can additionally be found on the github Functions branch. Within this branch, there is a folder "R" which contains the individual functions. This can similarly be downloaded and loaded into the R environment in case not all functions are desired. These functions are actively updated and may include more features than the in original tar.gz file. 

### Downloading the test data
# The COVID-19 data (~136 MB size of the zip file) can be found at the following link https://polybox.ethz.ch/index.php/s/fxQJ3NrRSwiPSSo This dataset contains VDJ (separate libraries for B and T cells) and GEX libraries from two convalescent COVID-19 patients.

# After downloading the zip file named "PlatypusTestData.zip", please unzip the file and find the path to the newly formed folder. Typically this will be in the Downloads folder, so the code below should work on MacOS. For Windows please uncomment the code and change the user name to match your PC.

directory_to_covid_patients_gex <- list()
directory_to_covid_patients_gex[[1]] <- c("~/Downloads/PlatypusTestData/Patient1_GEX/")
directory_to_covid_patients_gex[[2]] <- c("~/Downloads/PlatypusTestData/Patient2_GEX/")

# For Windows: 
#directory_to_covid_patients_gex[[1]] <- c("C:\Users\YourPCName\Downloads\PlatypusTestData\Patient1_GEX")

GEX_automate 允許我們在一行代碼中執(zhí)行標(biāo)準(zhǔn)的 Seurat 管道，并能夠更改 Seurat 使用的參數(shù)，包括最小讀取數(shù)、線粒體基因百分比、cluster分辨率等。此外，我們可以從 通過簡單地將 VDJ.gene.filter 參數(shù)設(shè)置為 TRUE，數(shù)據(jù)集不會讓克隆影響轉(zhuǎn)錄聚類。

covid_gex_patients_not_integrated <- Platypus::GEX_automate(GEX.outs.directory.list = directory_to_covid_patients_gex,integration.method = "scale.data",mito.filter = 20,cluster.resolution = 0.5,VDJ.gene.filter = T)

現(xiàn)在可以單獨可視化每個患者的二維圖。 默認(rèn)情況下，對象中包含 UMAP、PCA 和 TSNE 縮減。 在默認(rèn)參數(shù)下，這將顯示來自 Patient1 的 GEX 庫中的所有cell。 通過更改為第二個列表元素，我們可以查看來自 Patient2 的cell。

Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "umap")

Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "pca")

Seurat::DimPlot(covid_gex_patients_not_integrated[[1]],reduction = "tsne")

# UMAP for the second patient
Seurat::DimPlot(covid_gex_patients_not_integrated[[2]],reduction = "umap")

圖片.png

如果有興趣合并轉(zhuǎn)錄組（對應(yīng)于兩個不同的基因表達(dá)庫），則將這兩個目錄提供給單個列表元素。 GEX_automate 的輸出（這次對應(yīng)于單個 Seurat 對象列表元素）包含一個列，其中包含基于輸入目錄順序的樣本來源信息。 此外，如果存在生物重復(fù)，則可以按組標(biāo)記不同的樣本。 在這種情況下，只需將患者設(shè)置為 group1 和 group2 .

directory_to_covid_patients_integrated <- list()
directory_to_covid_patients_integrated[[1]] <- c("~/Downloads/PlatypusTestData/Patient1_GEX/",
                                                 "~/Downloads/PlatypusTestData/Patient2_GEX/")

## Here we use the previous version of automate_GEX to produce the identical UMAP seen in the manuscript. 
covid_gex <- Platypus::automate_GEX(GEX.outs.directory.list = directory_to_covid_patients_integrated[1:1],integration.method = "scale.data",mito.filter = 20,cluster.resolution = 0.5,VDJ.gene.filter = T)

Seurat::DimPlot(covid_gex[[1]],reduction = "umap")

圖片.png

當(dāng)然，后面還有標(biāo)準(zhǔn)的差異富集結(jié)果，大家應(yīng)該都很了解，我們接下來看看VDJ的分析。

VDJ Repertoire anaylsis

讀取數(shù)據(jù)

現(xiàn)在可以在集成 GEX 庫之前分析 VDJ 曲目數(shù)據(jù)。 如果只對 VDJ 文庫進(jìn)行了測序而沒有附帶的基因表達(dá)數(shù)據(jù)，這可能很有用。 首先根據(jù)默認(rèn)的 10x cellranger 克隆分型策略讀取克隆型，并增加包含克隆型信息的數(shù)據(jù)幀中的克隆信息量。 第一步是在 cellranger vdj 的輸出目錄上使用 VDJ_analyze 函數(shù)。 可以刪除較大的 BAM 文件以節(jié)省空間，因為當(dāng)前管道不需要這些文件。

#Read in the directories
VDJ.out.directory.list <- list()
VDJ.out.directory.list[[1]] <- "~/Downloads/PlatypusTestData/Patient1_BCR/"
VDJ.out.directory.list[[2]] <- "~/Downloads/PlatypusTestData/Patient2_BCR/"

#Run VDJ_analyze
covid_vdj_repertoire_bcells <- Platypus::VDJ_analyze(VDJ.out.directory =VDJ.out.directory.list, filter.1HC.1LC = T) 

VDJ_analyze 的輸出包含各種信息，包括哪些條形碼構(gòu)成克隆家族，nt_clone_ids（以防使用 VDJ_clonotype 函數(shù)更改克隆分型方法）。 此外，每個克隆家族的多數(shù)種系基因是從重疊群文件中提取的。

Changing the clonotype strategy

考慮到體細(xì)胞超突變可能發(fā)生在 CDR3 區(qū)域，通常配對的核苷酸 CDRH3 和 CDRL3 克隆分型可能不是最佳策略。因此，可能存在高度相似的克隆，它們可能結(jié)合正式屬于不同克隆家族的相同抗原。為了解決這個問題，我們添加了一個允許各種啟發(fā)式克隆分型策略的函數(shù)。這涉及通過相同的氨基酸 CDRH3 + CDRL3 序列、相同的種系使用或序列同源性要求進(jìn)行克隆分型。然后，這將通過使用用于克隆分型的新氨基酸（或其他策略）更新 clonotype_id 列來更新來自 VDJ_analyze 的原始克隆分型對象。同樣，新的克隆家族決定因素可以在 new_unique_clone 中找到。在下面的例子中，new_unique_clone 將包含粘貼在一起的重鏈和輕鏈的氨基酸序列。如果通過更改克隆分型策略合并多個核苷酸克隆，則可以在“條形碼”列中找到它們的條形碼。

covid_vdj_aminoacid_clonotype <- Platypus::VDJ_clonotype(clonotype.list=covid_vdj_repertoire_bcells,
                                                         clone.strategy = "cdr3.aa")

covid_vdj_germline_clonotype <- Platypus::VDJ_clonotype(clonotype.list=covid_vdj_repertoire_bcells,
                                                         clone.strategy = "hvj.lvj")

We can see that the clonotype_id column now contains the new ids based on this clonotyping strategy. We can always recover the nucleotide ids by looking at the nt_clones_ids.

Getting clonal information at the level of the single cell

到目前為止，這些函數(shù)都處于克隆水平，因此忽略了可能在克隆家族中有所不同的細(xì)胞特定特征，例如同種型、序列變異等。 VDJ_per_clone 函數(shù)將通過返回嵌套列表來提供此信息，其中外部list 元素對應(yīng)于曲目（例如，list[[1]] 是第一個目錄，作為 VDJ_analyze 函數(shù)的輸入）。內(nèi)部列表元素對應(yīng)于每個克隆的數(shù)據(jù)幀。在此數(shù)據(jù)框中，可以提取諸如重鏈和輕鏈的種系基因使用等信息（例如，HC_Vgene 是抗體 IGH 的 V 基因或 T 細(xì)胞 TRB 的 V 基因）。此外，可以根據(jù) cellranger 的輸出提取每個細(xì)胞條形碼的重鏈和輕鏈的完整序列。一個問題是從 cellranger 的 fasta 文件返回的序列超出了 FR1 到 FR4 區(qū)域（例如恒定區(qū)、信號肽），這可能對下游實驗驗證/表達(dá)造成問題。germline序列也會發(fā)生同樣的情況 - 再次由 cellranger 確定。germline信息可以在內(nèi)部列表元素的 full_HC_germline 和 full_LC_germline 列中找到。要獲得可克隆的 B/TCR 序列，可以繼續(xù)使用 call_MIXCR 函數(shù)或?qū)С鲂蛄胁⑹褂貌煌谋葘ぞ摺?/h5>

covid_single_cell <- Platypus::VDJ_per_clone(clonotype.list = covid_vdj_repertoire_bcells,VDJ.out.directory =VDJ.out.directory.list)

covid_single_clone_tcells <- Platypus::VDJ_per_clone(clonotype.list = covid_vdj_repertoire_tcells,VDJ.out.directory =VDJ.out.directory.list_TCR)

Extracting full-length sequences from the VDJRegion

為了量化體細(xì)胞變異的數(shù)量或提取用于表達(dá)的全長序列，從框架區(qū) 1 (FR1) 到框架區(qū) 4 (FR4) 的核苷酸序列通常很有用。使用 call_MIXCR 函數(shù)，可以將全長 VDJRegion 序列添加到克隆信息中，然后輕松提取。此功能適用于 UNIX/Mac，此外還要求 mixcr 已在本地下載（有許可協(xié)議）。只需在 call_MIXCR 函數(shù)中為可執(zhí)行文件提供目錄，如下所示。分別用于小鼠和人類的“mmu”或“hsa”。同樣，格式與輸入類似，其中外部列表對應(yīng)于單個曲目，內(nèi)部列表是具有各種信息的數(shù)據(jù)幀，包括全長 VDJ 序列（例如，VDJ.AA.LC 和 VDJ.AA .HC 表示輕鏈和重鏈氨基酸序列）。人們會注意到生殖系序列仍然很長（例如，在下面的示例中，“full_HC_germline”長度超過 600 個核苷酸）。這將使用單獨的函數(shù) VDJ_extract_germline 填充。

### WARNING: You will need to download MiXCR and change the mixcr.directory to the location of MiXCR
covid_vdj_region <- Platypus::call_MIXCR(VDJ.per.clone = covid_single_cell,mixcr.directory = "~/Downloads/mixcr-3.0.12/mixcr",species = "hsa")

Extracting full-length germline sequence corresponding to the VDJRegion

提取germline column 可以使用 VDJ_extract_germline 來完成，它采用由 cellranger 確定的germline column 。 輸出包含給定的數(shù)據(jù)幀格式的所有種系序列。 原始克隆型標(biāo)識符可以在 descrR1 列中找到repertoire。

extracted_covid_germline <- Platypus::VDJ_extract_germline(VDJ.per.clone=covid_single_cell,mixcr.directory="~/Downloads/mixcr-3.0.12/mixcr",extract.VDJRegion=T,species = "hsa")

Extracting trimmmed VDJ region sequence and germline

在 all_contig_annotations.json 文件可用的情況下，函數(shù) VDJ_per_clone 額外提供修剪過的 VDJ 序列和相應(yīng)的修剪過的種系序列。 這對于計算克隆型中單細(xì)胞的 VDJ 區(qū)域的體細(xì)胞超突變率非常有用。 名為“full_seq_”、“sequence_”和“trimmed_ref_”的列代表從cellranger的fasta文件返回的序列，在VDJ區(qū)域修剪的序列，在VDJ區(qū)域修剪的生殖系序列（后面是HC for Heavy Chain 和輕鏈的 LC）。 此數(shù)據(jù)不需要 MiXCR 對齊，但是，它不一定從 Framework1 開始，因此如果未來計劃涉及克隆 BCR/TCR，請謹(jǐn)慎使用。 已經(jīng)對此進(jìn)行了注釋，因為當(dāng)時測試數(shù)據(jù)不包括 JSON 文件。

#covid_single_clone_with_JSON <- VDJ_per_cell(clonotype.list = covid_vdj_repertoire_bcells,VDJ.out.directory =VDJ.out.directory.list,JSON = T)

# covid_single_cell was the object from VDJ_per_clone

#print(covid_single_clone_with_JSON[[1]][[1]]$trimmed_HC_sequence[1])  # trimmmed sequence 
#print(covid_single_clone_with_JSON[[1]][[1]]$trimmed_HC_germline[1])  # trimmmed germline sequence 

#stringdist::stringdist(covid_single_clone_with_JSON[[1]][[1]]$sequence_HC[1], VDJ.per.cell[[1]][[1]]$trimmed_ref_HC[1]) 

Organizing full-length sequences into clonal lineages easily exportable for phylogenetics.

在提取germline序列和全長 VDJRegion 序列后，我們可以結(jié)合這些信息并將序列分組為克隆譜系，這些譜系可以單獨導(dǎo)出用于下游系統(tǒng)發(fā)育學(xué)，或?qū)?Platypus 中的 VDJ_tree 函數(shù)。 輸出列表類似于其他 VDJ 函數(shù)的輸出，其中外部列表元素對應(yīng)于曲目，內(nèi)部列表元素對應(yīng)于各個克隆/克隆家族。 此處重鏈和輕鏈已粘貼在一起，但如果您想將它們分開，可以在 VDJ_extract_germline 函數(shù)的輸出中找到該信息。 例如，上面 VDJ_extract_germline 函數(shù)的輸出是extracted_covid_germline[[1]][[1]]$VDJ.AA.HC.LC[2]。 在這個序列中，重鏈和輕鏈序列的末端有一個“”。所以可以將“”處的字符串分開，并使用這個序列來創(chuàng)建系統(tǒng)發(fā)育樹。

covid_clonal_lineages <- Platypus::VDJ_clonal_lineages(call_MIXCR.output=covid_vdj_region, VDJ_extract_germline.output=extracted_covid_germline,as.nucleotide=F,with.germline=T)

Neighbor-joining phylogenetic trees

VDJ_clonal_lineage 函數(shù)的輸出可用作 VDJ_tree 的輸入，以生成快速、相鄰的系統(tǒng)發(fā)育樹。 考慮 VDJRegion（或來自 10x 的序列已被修剪）很重要，因為此函數(shù)不涉及多字符串對齊 - 而是從所有序列創(chuàng)建距離矩陣。 可以指定獨特序列的最小數(shù)量（例如，如果克隆譜系中有 50 個細(xì)胞，但只有一個獨特的抗體序列，則不會生成此樹）。 類似地，可以指定最大序列數(shù)，這將導(dǎo)致從所有序列中隨機(jī)采樣。 繪制第一棵樹，可以看到有一個變體（名為 Seq_1）有 90 個細(xì)胞（由 Freq_90 看到）并且它非常接近生殖系（Seq_6_Freq_1）。 我們知道它是種系，因為它是最高序列的變體，它被設(shè)置為根。

covid_trees <- Platypus::VDJ_tree(clonal.lineages = covid_clonal_lineages,with.germline=T,min.sequences = 5,max.sequences = 30,unique.sequences = T)

圖片.png

Plotting per-cell isotype information

其他 VDJ 功能包括繪制每個克隆的同種型分布。 這可以通過 VDJ_isotypes_per_clone 函數(shù)并設(shè)置要顯示的克隆數(shù)來執(zhí)行。 以下代碼提取并繪制了前 30 個克隆的同種型分布。 當(dāng)使用默認(rèn)的克隆分型策略時，可以看到第一個患者的前四個克隆中大部分有明確的 IgA 同種型。 還可以使用新的克隆分型策略來比較克隆定義的變化如何影響克隆擴(kuò)增譜。 只是提供 VDJ_clonotype 的輸出。

covid_isotypes <- Platypus::VDJ_isotypes_per_clone(VDJ_clonotype_output = covid_vdj_repertoire_bcells, VDJ_per_clone_output = covid_single_cell, clones = 30)

covid_isotypes_aa_clonotype <- Platypus::VDJ_isotypes_per_clone(VDJ_clonotype_output = covid_vdj_germline_clonotype, VDJ_per_clone_output = covid_single_cell, clones = 30)

covid_isotypes_germline_clonotype <- Platypus::VDJ_isotypes_per_clone(VDJ_clonotype_output = covid_vdj_aminoacid_clonotype, VDJ_per_clone_output = covid_single_cell, clones = 30)

print(covid_isotypes[[1]])

圖片.png

In the second patient (second repertoire entered in original VDJ_analyze) we also see a similarly expanded IgA clone with over 60 cell barcodes.

圖片.png

Sequence similarity networks

其他功能是專門為庫分析量身定制的 - 例如 VDJ_network，它通過將具有序列相似性的克隆連接起來，在庫之間或庫內(nèi)創(chuàng)建序列相似性網(wǎng)絡(luò)。 該函數(shù)依賴 igraph 來可視化地顯示和構(gòu)建圖形——這意味著具有大量序列的網(wǎng)絡(luò)將不容易顯示。 在以下示例中，從 VDJ_analyze 的輸出中取出前 60 個擴(kuò)展最多的克隆，并將其用作網(wǎng)絡(luò)構(gòu)建函數(shù)的輸入。 將 per.mouse 參數(shù)設(shè)置為 false 表示應(yīng)該為多個曲目生成一個網(wǎng)絡(luò)。

## Take the output from VDJ_analyze (or subsample as in this case the top 60 clones)
network_clones_covid <- list()
network_clones_covid[[1]] <- covid_vdj_repertoire_bcells[[1]][1:60,]
network_clones_covid[[2]] <- covid_vdj_repertoire_bcells[[2]][1:60,]
network_clones_covid <- list()
network_clones_covid[[1]] <- covid_vdj_repertoire_tcells[[1]][1:60,]
network_clones_covid[[2]] <- covid_vdj_repertoire_tcells[[2]][1:60,]

covid_bcell_igraph <- Platypus::VDJ_network(network_clones_covid[1:2],per.sample = F,distance.cutoff = 10,connected = F)
covid_bcell_igraph_d14 <- Platypus::VDJ_network(network_clones_covid[1:2],per.sample = F,distance.cutoff = 14,connected = F)

igraph::plot.igraph(covid_bcell_igraph[[4]],vertex.label=NA,vertex.size=7+(.06*covid_bcell_igraph[[2]]$frequency),vertex.color=factor(covid_bcell_igraph[[2]]$mouse))
igraph::plot.igraph(covid_bcell_igraph_d14[[4]],vertex.label=NA,vertex.size=7+(.06*covid_bcell_igraph_d14[[2]]$frequency),vertex.color=factor(covid_bcell_igraph_d14[[2]]$mouse))

圖片.png

Germline gene usage heatmaps

還可以在重鏈 V 基因和輕鏈 V 基因的背景下生成種系基因使用的熱圖。 VDJ_Vgene_usage 函數(shù)的輸出是每個曲目的矩陣，對應(yīng)于 VDJ_analyze 指定的順序。 外表對應(yīng)樣本，內(nèi)表對應(yīng)矩陣，行對應(yīng)重鏈V基因，列對應(yīng)V基因輕鏈。 因此，輸出[[1]][i,j] 對應(yīng)于使用 IGH-Vgene[i] 和 IGK/L-Vgene[j] 組合的克隆數(shù)。

#First calculate adjacency matrix for V gene usage
covid_Vgene_usage <- Platypus::VDJ_Vgene_usage(VDJ.clonotype.output = covid_vdj_repertoire_bcells)
library(pheatmap)

pheatmap::pheatmap(covid_Vgene_usage[[1]],show_rownames = F,show_colnames = F)

圖片.png

然后可以輕松地將其繪制為熱圖以觀察曲目之間的模式，或者可用于使用“pheatmap”包計算 V 基因相關(guān)性。

Platypus 還允許單獨分析 HC 和 LC 的 V 基因使用情況。 VDJ_Vgene_usage_barplot 允許用戶繪制最常用的 IgH 或 IgK/LV 基因。 默認(rèn)情況下，此功能僅提供 HC V 基因的可視化，但如果是 LC，也可以提供 LC。 Vgene 設(shè)置為 TRUE。 還可以選擇要描繪的最常用基因的數(shù)量。

covid_Vgene_usage_barplot <- Platypus::VDJ_Vgene_usage_barplot(covid_vdj_repertoire_bcells, HC.gene.number = 10, LC.Vgene = T, LC.gene.number = 10)
print(covid_Vgene_usage_barplot[[1]])

圖片.png

example.vdj.vgene_usage <- Platypus::VDJ_Vgene_usage_stacked_barplot(clonotype.list = covid_vdj_repertoire_bcells, LC.Vgene = F,HC.gene.number = 10, Fraction.HC = 1)
example.vdj.vgene_usage[[1]]

圖片.png

此外，我們還可以生成 V 和 J 基因如何在整個曲目中組合的圓形可視化。 在接下來的示例中，我們使用 VDJ_VJ_usage_circos 來查看每個擴(kuò)展克隆型的 V 基因和相應(yīng)的 J 基因。

vj_circos_bcells <- Platypus::VDJ_VJ_usage_circos(covid_vdj_repertoire_bcells[2:2], c.threshold = 1,label.threshold=50,cell.level = T)

圖片.png

vj_circos_tcells <- Platypus::VDJ_VJ_usage_circos(covid_vdj_repertoire_tcells[1:1], c.threshold = 1,label.threshold=50,cell.level = T)

圖片.png

Assessing CDR3 sequence similarity

最后，還可以查看不同克隆中出現(xiàn)的任何特定 HC 和 LC CDR3 氨基酸模式。 使用 VDJ_logoplot 函數(shù)，用戶可以繪制特定長度的 CDR3 區(qū)域的徽標(biāo)圖，由 length_cdr3 參數(shù)指定。 例如，下面的標(biāo)志圖對應(yīng)于所有長度為 25 的 CDR3 氨基酸序列

covid_CDR3_logoplot <- Platypus::VDJ_logoplot(VDJ.object = covid_vdj_repertoire_bcells, length_cdr3 = 25)

圖片.png

Integrating repertoire and gene expression

當(dāng)前 5' 測序的優(yōu)勢在于基因表達(dá) (GEX) 和庫 (VDJ) 庫是從同一樣本中提取的，然后可以將其鏈接回以證明給定的 T 細(xì)胞具有特定的基因表達(dá)模式和 還有某種T細(xì)胞受體序列。 以下函數(shù)旨在整合這兩條信息。

Integrating transcriptional clusters to the VDJ objects

可能會問的一件事是給定克隆家族中的 B 或 T 細(xì)胞在轉(zhuǎn)錄水平上有多相似。 我們可以做到這一點的一種方法是利用在automatic_GEX 函數(shù)中執(zhí)行的轉(zhuǎn)錄聚類，然后將此信息提供給曲目信息。 這可以在克隆級別（例如，以 VDJ_analyze 的輸出格式）或單元級別（例如，以 VDJ_per_clone 的輸出格式）完成。 下面是克隆水平的例子。 輸出與 VDJ_analyze 輸出相同，除了現(xiàn)在有對應(yīng)于多數(shù)轉(zhuǎn)錄組簇的列（例如，如果克隆中的大多數(shù)細(xì)胞在簇 1 中，則該列的值為 1）。 在此旁邊有一個列，告訴您該克隆中所有單元格的群集百分比。 最后，GEX 對象中有單元格索引，因此您可以查看 GEX 對象中找到了哪些單元格并可以顯式調(diào)用它們。

covid_integrating_clonal_level <- Platypus::VDJ_GEX_integrate(GEX.object = covid_gex[[1]], 
                                                 clonotype.list =  covid_vdj_repertoire_bcells,
                                                 VDJ.per.clone = covid_single_cell,
                                                 clonotype.level = TRUE)
covid_integrating_clonal_level_tcell <- Platypus::VDJ_GEX_integrate(GEX.object = covid_gex[[1]], 
                                                 clonotype.list =  covid_vdj_repertoire_tcells,
                                                 VDJ.per.clone = covid_single_cell,
                                                 clonotype.level = TRUE)

covid_integrating_cell_level <- Platypus::VDJ_GEX_integrate(GEX.object = covid_gex[[1]], 
                                                 clonotype.list =  covid_vdj_repertoire_bcells[1:1],
                                                 VDJ.per.clone = covid_single_cell[1:1],
                                                 clonotype.level = FALSE)

Relating clonal expansion to transcriptional cluster membership

VDJ_GEX_expansion 函數(shù)采用先前集成的 VDJ 和 GEX 數(shù)據(jù)集，然后繪制指定克隆的轉(zhuǎn)錄簇成員分布。 在下面的示例中，繪制了前 20 個克隆及其集群成員的分布。

covid_clonotype_clusters_plot <-  Platypus::VDJ_GEX_expansion(GEX.list=covid_gex[[1]],
                                                              VDJ.GEX.integrate.list = covid_integrating_clonal_level,
                                                              highlight.isotype = "None",
                                                              highlight.number=1:10)

covid_clonotype_clusters_plot_tcells <-  Platypus::VDJ_GEX_expansion(GEX.list=covid_gex[[1]],
                                                              VDJ.GEX.integrate.list = covid_integrating_clonal_level_tcell,
                                                              highlight.isotype = "None",
                                                              highlight.number=1:10)

圖片.png

We see that some clones did not have any barcodes overlapping between the GEX and VDJ data sets for this particular patient and therefore there is no bar at that given clonal index on the X axis.

Visualizing clones on the 2 dimensional landscape

我們想要在我們之前使用 automatic_GEX 計算的 tSNE/UMAP 空間上覆蓋克隆，這可以使用visualize_clones_GEX 函數(shù)來完成。 以下代碼將在使用automation_GEX 函數(shù)計算的 UMAP 空間上覆蓋前 10 個克隆。 該組對應(yīng)于克隆順序 - 例如，組 1 是集成 VDJ 列表中的第一個克隆家族。 在這里我們可以看到 T 細(xì)胞簇中有幾個細(xì)胞，這表明在細(xì)胞捕獲過程中存在雙胞胎，因為在這個實驗中，B 和 T 細(xì)胞匯集在一個反應(yīng)中。

covid_top10_umap <- Platypus::GEX_visualize_clones(GEX.list=covid_gex,
                     VDJ.GEX.integrate.list=covid_integrating_clonal_level_tcell[2:2],
                     highlight.type="clonotype",
                     highlight.number=1:10,
                     reduction="umap")
print(covid_top10_umap[[1]]) 

圖片.png

Specific gene expression information on the clonal level

之前將 GEX 信息集成到 VDJ 輸出的格式中。 然而，可能想知道基因表達(dá)如何尋找某些克隆型（例如，頂級克隆中有多少細(xì)胞表達(dá) Fas/GL7 或生發(fā)中心細(xì)胞標(biāo)記）。 因此，需要將 VDJ 信息和條形碼集成到 GEX 對象中。 這可以通過使用 clonotype_GEX 函數(shù)和 GEX_heatmap 函數(shù)來完成。 為了節(jié)省空間，可以將此函數(shù)的輸出分配給原始的automatic_GEX 輸出，因為它將是另一個 Seurat 對象。 通過將 b.or.t 設(shè)置為“b”，將顯示一組預(yù)先選擇的 B 細(xì)胞基因。 還可以提供定制基因來測試特定的克隆型。 組/列顏色對應(yīng)于給定的克隆家族。

covid_gex_integrate <- Platypus::GEX_clonotype(GEX.object=covid_gex[[1]],VDJ.per.clone = covid_single_clone_tcells)

covid_gex_integrate_bcells <- Platypus::GEX_clonotype(GEX.object=covid_gex[[1]],VDJ.per.clone = covid_single_cell)

GEX_phenotype_per_clone_plot <- Platypus::GEX_phenotype_per_clone(seurat.object = covid_gex_integrate, clonotype.ids= c(1,2,3,4,5))

print(GEX_phenotype_per_clone_plot)

圖片.png

covid_tcell_gene_heatmap <- Platypus::GEX_heatmap(covid_gex_integrate,b.or.t = "t",clone.rank.threshold = 20,sample.index = 1)

covid_bcell_gene_heatmap <- Platypus::GEX_heatmap(covid_gex_integrate_bcells,b.or.t = "b",clone.rank.threshold = 20,sample.index = 1)

Integrating clonal lineage information and transcriptional cluster information

我們可以將轉(zhuǎn)錄簇信息提供到更早組織的克隆譜系中。 這可以使用 VDJ_GEX_clonal_lineage_clusters 函數(shù)來完成，該函數(shù)會將轉(zhuǎn)錄簇信息放入 clonal_lineage 對象的名稱中。 由于只在 covid_integrating_cell_level 對象中整合了第一個患者，因此將只獲取第一批患者的克隆譜系。 現(xiàn)在，如果打印第一個克隆譜系中第一個序列的名稱，會看到一個 cluster5 已添加到名稱的末尾。 這表明這個特定的 B 細(xì)胞是在原始自動化 GEX 的 cluster5 中發(fā)現(xiàn)的。

clonal_lineage_integrate <- Platypus::VDJ_GEX_clonal_lineage_clusters(covid_integrating_cell_level,covid_clonal_lineages[1:1])

最后再強(qiáng)調(diào)一句，一定要關(guān)注VDJ的氨基酸序列

生活很好，有你更好