单细胞免疫组库数据分析||Seurat整合单细胞转录组与VDJ数据

生信技能树jimmy

发布于 2020-06-05 10:03:22

3.6K0

发布于 2020-06-05 10:03:22

文章被收录于专栏：单细胞天地单细胞天地

分享是一种态度

作者 | 周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

在做10X单细胞免疫组库分析的是往往是做一部分BCR、TCR做一部分5‘转录组，那么怎样才能把两者结合到一起呢？

今天我们尝试用我们的趁手工具做一下整合分析。

首先是下载数据，我们从10X官方的dataset中下载数据：https://support.10xgenomics.com/single-cell-vdj/datasets/3.1.0/vdj_v1_hs_pbmc3

在下载页面有关于这个样本的基本介绍，如这个数据集根据单细胞V(D)J试剂试剂盒使用指南和细胞表面蛋白特征条形码技术(CG000186)，从标记的细胞中扩增出cDNA，生成5’基因表达、细胞表面蛋白、TCR富集和Ig富集文库。这里是不是应该把本文的题目改成：Seurat整合单细胞转录组与VDJ数据与膜蛋白数据啊。

读入数据：

1rm(list =ls())
2library(Seurat)
3data <-   Read10X(data.dir = "F:\\VDJ\\filtered_feature_bc_matrix")
410X data contains more than one type and is being returned as a list containing matrices of each type.

上来就报错：

1seurat_object = CreateSeuratObject(counts = data)
2Error in CreateAssayObject(counts = counts, min.cells = min.cells, min.features = min.features) : 
3  No cell names (colnames) names present in the input matrix
4>

原来我们的feature中不仅有基因表达还有膜蛋白的信息，分开来读好了：

1seurat_object = CreateSeuratObject(counts = data$`Gene Expression`)
2seurat_object[['Protein']] = CreateAssayObject(counts = data$`Antibody Capture`)

我们熟悉的seurat对象又回来了:

1seurat_object
2
3An object of class Seurat 
433555 features across 7231 samples within 2 assays 
5Active assay: RNA (33538 features)
6 1 other assay present: Protein

注意这里是两个数据， RNA ，Protein：

当我们需要分析哪一部分的时候，用DefaultAssay(seurat_object)来指定，这里我们只分析RNA的数据啊，茂密的森林伸出两条路，我们也看一眼膜蛋白的数据吧：

1 seurat_object@assays$Protein
2Assay data with 17 features for 7231 cells
3First 10 features:
4 CD3-UCHT1-TotalC, CD19-HIB19-TotalC, CD45RA-HI100-TotalC, CD4-RPA-T4-TotalC, CD8a-RPA-T8-TotalC, CD14-M5E2-TotalC, CD16-3G8-TotalC, CD56-QA17A16-TotalC, CD25-BC96-TotalC,
5CD45RO-UCHL1-TotalC

接着我们读入BCR的数据：

 1bcr <- read.csv("F:\\VDJ\\filtered_contig_annotations.csv")
 2bcr[1:4,]
 3 barcode is_cell                   contig_id high_confidence length chain    v_gene d_gene j_gene c_gene full_length productive             cdr3
 41 AAACCTGTCACTGGGC-1    True AAACCTGTCACTGGGC-1_contig_1            True    556   IGK IGKV1D-33   None  IGKJ4   IGKC        True       True      CQQYDNLPLTF
 52 AAACCTGTCACTGGGC-1    True AAACCTGTCACTGGGC-1_contig_2            True    516   IGH  IGHV4-59   None  IGHJ5   IGHM        True       True     CARGGNSGLDPW
 63 AAACCTGTCACTGGGC-1    True AAACCTGTCACTGGGC-1_contig_3            True    569   IGK IGKV2D-30   None  IGKJ1   IGKC        True      False CWGLLLHARYTLAWTF
 74 AAACCTGTCAGGTAAA-1    True AAACCTGTCAGGTAAA-1_contig_1            True    548   IGK  IGKV1-12   None  IGKJ4   IGKC        True       True      CQQANSFPLTF
 8                                           cdr3_nt reads umis raw_clonotype_id       raw_consensus_id
 91                TGTCAACAGTATGATAATCTCCCGCTCACTTTC  7410   61       clonotype7 clonotype7_consensus_1
102             TGTGCGAGAGGCGGGAACAGTGGCTTAGACCCCTGG  1458   18       clonotype7 clonotype7_consensus_2
113 TGTTGGGGTTTATTACTGCATGCAAGGTACACACTGGCCTGGACGTTC  2045   17       clonotype7                   None
124                TGTCAACAGGCTAACAGTTTCCCGCTCACTTTC   140    1       clonotype2 clonotype2_consensus_2

为了保持barcode一致，处理一下barcode的命名：

 1bcr$barcode <- gsub("-1", "", bcr$barcode)
 2bcr <- bcr[!duplicated(bcr$barcode), ]
 3names(bcr)[names(bcr) == "raw_clonotype_id"] <- "clonotype_id"
 4table(bcr$v_gene)
 5
 6    IGHV1-18      IGHV1-2     IGHV1-24      IGHV1-3     IGHV1-45     IGHV1-46     IGHV1-58     IGHV1-69    IGHV1-69D      IGHV1-8 IGHV1/OR15-1 IGHV1/OR21-1     IGHV2-26      IGHV2-5     IGHV2-70 
 7          30           48            5           14            1           18            2            1           26            8            1            1            9           28            8 
 8    IGHV3-11     IGHV3-13     IGHV3-15     IGHV3-16     IGHV3-21     IGHV3-23     IGHV3-30     IGHV3-33     IGHV3-43     IGHV3-48     IGHV3-49     IGHV3-53     IGHV3-64    IGHV3-64D     IGHV3-66 
 9          20            7           29            1           43           95           40           50           10           25           10           13            5            4            2 
10     IGHV3-7     IGHV3-72     IGHV3-73     IGHV3-74   IGHV4-30-2     IGHV4-31     IGHV4-34     IGHV4-39      IGHV4-4     IGHV4-59     IGHV4-61   IGHV5-10-1     IGHV5-51      IGHV6-1    IGHV7-4-1 
11          15            6            7           25            9            4           30           55           22           76            9            5           21            4            8 
12    IGKV1-12     IGKV1-16     IGKV1-17     IGKV1-27     IGKV1-33     IGKV1-37     IGKV1-39      IGKV1-5      IGKV1-6      IGKV1-8      IGKV1-9    IGKV1D-13    IGKV1D-33    IGKV1D-37    IGKV1D-39 
13          29           10            8            9            9            2            6           42            3           16           11            1           15            6           65 
14   IGKV1D-43     IGKV1D-8     IGKV2-24     IGKV2-28     IGKV2-30    IGKV2D-26    IGKV2D-28    IGKV2D-29    IGKV2D-30    IGKV2D-40     IGKV3-11     IGKV3-15     IGKV3-20      IGKV3-7    IGKV3D-11 
15           5            4            6           26           14            5            1            9            1            2           39           58           89            2            1 
16   IGKV3D-15    IGKV3D-20      IGKV4-1      IGKV5-2     IGKV6-21    IGKV6D-21     IGLV1-40     IGLV1-44     IGLV1-47     IGLV1-51    IGLV10-54     IGLV2-11     IGLV2-14     IGLV2-18     IGLV2-23 
17           9            2           60            8            2            1           31           37           26           24            3           24           60            1           20 
18     IGLV2-8      IGLV3-1     IGLV3-10     IGLV3-12     IGLV3-16     IGLV3-19     IGLV3-21     IGLV3-25     IGLV3-27      IGLV3-9      IGLV4-3     IGLV4-60     IGLV4-69     IGLV5-37     IGLV5-39 
19          20           35           15            1            2           15           27           16            4            4            2            5           15            7            1 
20    IGLV5-45     IGLV5-48     IGLV5-52     IGLV6-57     IGLV7-43     IGLV7-46     IGLV8-61     IGLV9-49         None 
21           9            1            2           14            5           10           15            9          173 
22>

读入克隆型信息：

 1clono <- read.csv("F:\\VDJ\\vdj_v1_hs_pbmc3_b_clonotypes.csv")
 2head(clono)
 3
 4 clonotype_id frequency  proportion                            cdr3s_aa                                                                                cdr3s_nt
 51   clonotype1         7 0.009067358                    IGK:CQQYDNWPPYTF                                                IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
 62   clonotype2         3 0.003886010    IGK:CQQANSFPLTF;IGK:CQQYDNWPPYTF          IGK:TGTCAACAGGCTAACAGTTTCCCGCTCACTTTC;IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
 73   clonotype5         2 0.002590674                     IGL:CQARDSSTVVF                                                   IGL:TGTCAGGCGCGGGACAGCAGCACTGTGGTATTC
 84   clonotype6         2 0.002590674 IGH:CARPGTTGTTGLKNW;IGK:CQQYNNWPLTF IGH:TGTGCGAGACCCGGTACAACTGGAACGACGGGTTTAAAAAACTGG;IGK:TGTCAGCAGTATAATAACTGGCCTCTCACCTTC
 95   clonotype4         2 0.002590674 IGH:CAKSFFTGTGQFHYW;IGK:CQQYSTYSQTF IGH:TGTGCGAAATCATTTTTTACCGGGACAGGACAGTTTCACTATTGG;IGK:TGCCAACAGTATAGTACTTATTCTCAAACGTTC
106   clonotype3         2 0.002590674                    IGL:CQSYDSSLSGVF                                                IGL:TGCCAGTCCTATGACAGCAGCCTGAGTGGGGTGTTC

整合到seurat之中：

 1# Slap the AA sequences onto our original table by clonotype_id.
 2bcr <- merge(bcr, clono)
 3
 4# Reorder so barcodes are first column and set them as rownames.
 5rownames(bcr) <- bcr[,2]
 6#tcr[,1] <- NULL
 7
 8# Add to the Seurat object's metadata.
 9clono_seurat <- AddMetaData(object=seurat_object, metadata=bcr)
10
11head(clono_seurat@meta.data)
12                    orig.ident nCount_RNA nFeature_RNA nCount_Protein nFeature_Protein clonotype_id barcode is_cell contig_id high_confidence length chain v_gene d_gene j_gene c_gene full_length
13AAACCTGAGATCTGAA SeuratProject       4386         1206           5092               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
14AAACCTGAGGAACTGC SeuratProject       6258         1723           6480               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
15AAACCTGAGGAGTCTG SeuratProject       4243         1404           3830               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
16AAACCTGAGGCTCTTA SeuratProject       5205         1622           3589               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
17AAACCTGAGTACGTTC SeuratProject       7098         2123           5250               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
18AAACCTGGTATCACCA SeuratProject       1095          604           3573               17         <NA>    <NA>    <NA>      <NA>            <NA>     NA  <NA>   <NA>   <NA>   <NA>   <NA>        <NA>
19                 productive cdr3 cdr3_nt reads umis raw_consensus_id frequency proportion cdr3s_aa cdr3s_nt
20AAACCTGAGATCTGAA       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>
21AAACCTGAGGAACTGC       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>
22AAACCTGAGGAGTCTG       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>
23AAACCTGAGGCTCTTA       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>
24AAACCTGAGTACGTTC       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>
25AAACCTGGTATCACCA       <NA> <NA>    <NA>    NA   NA             <NA>        NA         NA     <NA>     <NA>

为什么这么多NA？因为转录组中的细胞很多，而免疫细胞（BCR）只有一部分。

 1 table(rownames(clono_seurat@meta.data)  %in% rownames(bcr))
 2
 3FALSE  TRUE 
 4 6469   762 
 5> clono_seurat<- subset(clono_seurat,cells =  rownames(bcr))
 6> head(clono_seurat@meta.data)
 7                    orig.ident nCount_RNA nFeature_RNA nCount_Protein nFeature_Protein clonotype_id          barcode is_cell                   contig_id high_confidence length chain   v_gene d_gene
 8CTTTGCGCAAGGACTG SeuratProject       1266          642            696               17   clonotype1 CTTTGCGCAAGGACTG    True CTTTGCGCAAGGACTG-1_contig_1            True    479   IGK IGKV3-15   None
 9CAGCTAACACCGATAT SeuratProject        935          484            545               17   clonotype1 CAGCTAACACCGATAT    True CAGCTAACACCGATAT-1_contig_1            True    549   IGK IGKV3-15   None
10TATTACCGTCGCTTCT SeuratProject       1352          594            422               17   clonotype1 TATTACCGTCGCTTCT    True TATTACCGTCGCTTCT-1_contig_1            True    565   IGK IGKV3-15   None
11CTGAAGTGTGAAGGCT SeuratProject       1259          626            450               17   clonotype1 CTGAAGTGTGAAGGCT    True CTGAAGTGTGAAGGCT-1_contig_1            True    554   IGK IGKV3-15   None
12CCATGTCAGACCGGAT SeuratProject       1088          555            741               17   clonotype1 CCATGTCAGACCGGAT    True CCATGTCAGACCGGAT-1_contig_1            True    556   IGK IGKV3-15   None
13CGTTCTGAGCTCAACT SeuratProject        808          461           4896               17   clonotype1 CGTTCTGAGCTCAACT    True CGTTCTGAGCTCAACT-1_contig_1            True    565   IGK IGKV3-15   None
14                 j_gene c_gene full_length productive         cdr3                              cdr3_nt reads umis       raw_consensus_id frequency  proportion         cdr3s_aa
15CTTTGCGCAAGGACTG  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   197    2 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
16CAGCTAACACCGATAT  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   274    2 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
17TATTACCGTCGCTTCT  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   299    2 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
18CTGAAGTGTGAAGGCT  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   250    3 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
19CCATGTCAGACCGGAT  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   324    2 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
20CGTTCTGAGCTCAACT  IGKJ2   IGKC        True       True CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   187    3 clonotype1_consensus_1         7 0.009067358 IGK:CQQYDNWPPYTF
21                                                 cdr3s_nt
22CTTTGCGCAAGGACTG IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
23CAGCTAACACCGATAT IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
24TATTACCGTCGCTTCT IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
25CTGAAGTGTGAAGGCT IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
26CCATGTCAGACCGGAT IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT
27CGTTCTGAGCTCAACT IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT

seurat标准过程：

 1library(tidyverse)
 2clono_seurat %>% NormalizeData() %>% FindVariableFeatures() %>% 
 3 ScaleData(assay="RNA")  %>% 
 4 RunPCA(npcs = 50,assay="RNA")%>% 
 5 FindNeighbors(assay="RNA")%>%
 6 FindClusters(assay="RNA")%>%
 7 RunUMAP(1:30) -> clono_seurat
 8
 9head(clono_seurat@meta.data)
10
11                  orig.ident nCount_RNA nFeature_RNA nCount_Protein nFeature_Protein clonotype_id          barcode is_cell
12CTTTGCGCAAGGACTG SeuratProject       1266          642            696               17   clonotype1 CTTTGCGCAAGGACTG    True
13CAGCTAACACCGATAT SeuratProject        935          484            545               17   clonotype1 CAGCTAACACCGATAT    True
14TATTACCGTCGCTTCT SeuratProject       1352          594            422               17   clonotype1 TATTACCGTCGCTTCT    True
15CTGAAGTGTGAAGGCT SeuratProject       1259          626            450               17   clonotype1 CTGAAGTGTGAAGGCT    True
16CCATGTCAGACCGGAT SeuratProject       1088          555            741               17   clonotype1 CCATGTCAGACCGGAT    True
17CGTTCTGAGCTCAACT SeuratProject        808          461           4896               17   clonotype1 CGTTCTGAGCTCAACT    True
18                                  contig_id high_confidence length chain   v_gene d_gene j_gene c_gene full_length productive
19CTTTGCGCAAGGACTG CTTTGCGCAAGGACTG-1_contig_1            True    479   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
20CAGCTAACACCGATAT CAGCTAACACCGATAT-1_contig_1            True    549   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
21TATTACCGTCGCTTCT TATTACCGTCGCTTCT-1_contig_1            True    565   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
22CTGAAGTGTGAAGGCT CTGAAGTGTGAAGGCT-1_contig_1            True    554   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
23CCATGTCAGACCGGAT CCATGTCAGACCGGAT-1_contig_1            True    556   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
24CGTTCTGAGCTCAACT CGTTCTGAGCTCAACT-1_contig_1            True    565   IGK IGKV3-15   None  IGKJ2   IGKC        True       True
25                        cdr3                              cdr3_nt reads umis       raw_consensus_id frequency  proportion
26CTTTGCGCAAGGACTG CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   197    2 clonotype1_consensus_1         7 0.009067358
27CAGCTAACACCGATAT CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   274    2 clonotype1_consensus_1         7 0.009067358
28TATTACCGTCGCTTCT CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   299    2 clonotype1_consensus_1         7 0.009067358
29CTGAAGTGTGAAGGCT CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   250    3 clonotype1_consensus_1         7 0.009067358
30CCATGTCAGACCGGAT CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   324    2 clonotype1_consensus_1         7 0.009067358
31CGTTCTGAGCTCAACT CQQYDNWPPYTF TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT   187    3 clonotype1_consensus_1         7 0.009067358
32                        cdr3s_aa                                 cdr3s_nt RNA_snn_res.0.8 seurat_clusters
33CTTTGCGCAAGGACTG IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5
34CAGCTAACACCGATAT IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5
35TATTACCGTCGCTTCT IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5
36CTGAAGTGTGAAGGCT IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5
37CCATGTCAGACCGGAT IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5
38CGTTCTGAGCTCAACT IGK:CQQYDNWPPYTF IGK:TGTCAGCAGTATGATAACTGGCCTCCGTACACTTTT               5               5

1p1 <- DimPlot(clono_seurat)
2p2 <-DimPlot(clono_seurat,group.by = "chain")
3p3 <-FeaturePlot(clono_seurat,features="frequency")
4
5library(patchwork)
6p1+p2+p3

1library(sankeywheel)
2sankeywheel(
3  from = clono_seurat@meta.data$chain, to = clono_seurat@meta.data$v_gene,
4  weight = clono_seurat@meta.data$nCount_Protein,#type = "sankey",
5  title = "sk",
6  subtitle = "chain &  v_gene",
7  seriesName = "", width = "100%", height = "600px"
8)

计算不同链的差异基因：

1DT::datatable(FindMarkers(clono_seurat,ident.1 = 'IGK',group.by = "chain"))
2mk<-FindMarkers(clono_seurat,ident.1 = 'IGK',group.by = "chain")

1DotPlot(clono_seurat,features=rownames(mk),group.by = "chain") + RotatedAxis()

1RidgePlot(clono_seurat,features=rownames(mk),group.by = "chain")

References

[1] Tutorial: Integrating VDJ sequencing data with Seurat: https://links.jianshu.com/go?to=https%3A%2F%2Fwww.biostars.org%2Fp%2F384640%2F [2] Tutorial: Associating VDJ clonotyping data with scRNA-seq in Seurat: https://links.jianshu.com/go?to=https%3A%2F%2Fwww.biostars.org%2Fp%2F383217%2F