胃癌單細胞數據集GSE163558複現(十)：CopyKAT分析

前言

Hello小夥伴們大家好，我是生信技能樹的小學徒”我才不吃蛋黃“。今天是胃癌單細胞數據集複現系列第十期。第九期我們用TCGA-STAD數據進行生存分析。本期，我們将回到高級分析（分析），通過計算CNV評估上皮細胞良惡性。

1.背景介紹

在第六期推文中，我們聯合TCGA-STAD數據給上皮添加了惡性評分，從而協助區分腫瘤細胞與非腫瘤細胞。在高級分析中，我們可以利用和對單細胞進行CNV(Copy , 拷貝數變異)分析，區分腫瘤細胞與非腫瘤細胞。

CNV是基因組結構變異( ，SV)的重要組成部分，主要原因是基因組發生重排,一般指長度為1kb以上的基因組大片段的拷貝數增加或者減少,主要表現為亞顯微水平的缺失和重複。CNV在腫瘤的發生和發展研究中扮演重要的角色。

本期我們使用 ( of )進行CNV分析，它是一種集成貝葉斯分割方法。利用scRNA-seq數據推斷人類腫瘤的基因組拷貝數，然後通過對拷貝數數據進行分層聚類，以識别非整倍體腫瘤細胞和二倍體基質細胞之間的最大距離。其原理同樣是通過單細胞轉錄組數據來推斷細胞的染色體倍數，進而推斷是正常細胞（）還是腫瘤細胞（）。

2.數據分析2.1 導入數據

首先清除系統環境變量，加載R包，設置新文件夾和工作目錄，導入惡性上皮細胞數據：

rm(list=ls())
options(stringsAsFactors=F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(infercnv)
library(copykat)
library(tidyverse)
#devtools::install_github("broadinstitute/infercnv")
dir.create("8-CNV")
getwd()
setwd("../8-CNV")
sce=readRDS("../6-TCGA_STAD/malignant.rds")
table(sce$celltype)
Idents(sce)=sce$celltype

避免後面流程運行的太長了對細胞進行抽樣：

seurat_object=sce
seurat_object<-subset(seurat_object,downsample=200)
table(Idents(seurat_object))

預測腫瘤/正常細胞狀态的基本原理是非整倍體在人類癌症中很常見 (90%)。具有廣泛全基因組拷貝數畸變（非整倍體）的細胞被認為是腫瘤細胞，而基質正常細胞和免疫細胞通常具有2N二倍體或接近二倍體的拷貝數分布。通過單細胞轉錄組數據來推斷細胞的染色體倍數，進而推斷是正常細胞（）還是腫瘤細胞（）。它還可以進一步對腫瘤細胞進行聚類，找出不同的亞群。

運行

ngene.chr參數是過濾細胞的一個标準，它要求被用來做CNV預測的細胞，一個染色體上至少有5個基因。

sam.name定義樣本名稱 ( name)，會給出來的文件加前綴

scRNA<-seurat_object
counts<-as.matrix(scRNA@assays$RNA$counts)
table(scRNA$celltype)
if(T){cnv<-copykat(rawmat=counts,ngene.chr=5,sam.name="test",n.cores=8)
saveRDS(cnv,"cnv.rds")}

添加結果到對象meta.data信息中：

cnv=readRDS("cnv.rds")
table(rownames(cnv$CNAmat))
a=cnv$prediction$copykat.pred
table(a)
scRNA$CopyKAT=a

結果可視化：正常細胞（，藍色）還是腫瘤細胞（，紅色）

p1<-DimPlot(scRNA,group.by="celltype",label=T,reduction='tsne')
p2<-DimPlot(scRNA,group.by="CopyKAT",reduction='tsne')+scale_color_manual(values=c("#F8766D",'#02BFC4',"gray"))
pc<-p1+p2
pc

可視化maker基因'TPPP3','KRT17'的表達情況，檢查結果準确性：

cols=c("gray","coral2")
plot<-FeaturePlot(scRNA,features=c('TPPP3','KRT17'),cols=cols,pt.size=1,reduction='tsne')+
theme(panel.border=element_rect(fill=NA,color="black",size=1,linetype="solid"))#加邊框
plot_grid(p2,plot)

結果主要有兩個：

1）預測的結果正常細胞（）還是腫瘤細胞（）；

2） 每個CNV 在每個細胞的表達量。這裡看出和的不同來了，是基于 而不是gene level的表達量。

繪制熱圖

copykat.test=cnv
pred.test<-data.frame(copykat.test$prediction)
CNA.test<-data.frame(copykat.test$CNAmat)

head(pred.test)
head(CNA.test[,1:5])
my_palette<-colorRampPalette(rev(RColorBrewer::brewer.pal(n=3,name="RdBu")))(n=999)
chr<-as.numeric(CNA.test$chrom)%%2+1
rbPal1<-colorRampPalette(c('black','grey'))
CHR<-rbPal1(2)[as.numeric(chr)]
chr1<-cbind(CHR,CHR)
rbPal5<-colorRampPalette(RColorBrewer::brewer.pal(n=8,name="Dark2")[2:1])
com.preN<-pred.test$copykat.pred
pred<-rbPal5(2)[as.numeric(factor(com.preN))]
cells<-rbind(pred,pred)
col_breaks=c(seq(-1,-0.4,length=50),seq(-0.4,-0.2,length=150),seq(-0.2,0.2,length=600),seq(0.2,0.4,length=150),seq(0.4,1,length=50))
heatmap.3(t(CNA.test[,4:ncol(CNA.test)]),dendrogram="r",distfun=function(x)parallelDist::parDist(x,threads=4,method="euclidean"),
hclustfun=function(x)hclust(x,method="ward.D2"),
ColSideColors=chr1,Colv=NA,Rowv=TRUE,
notecol="black",col=my_palette,breaks=col_breaks,key=TRUE,
keysize=1,density.info="none",trace="none",
cexRow=0.1,cexCol=0.1,cex.main=1,cex.lab=0.1,
symm=F,symkey=F,symbreaks=T,cex=1,cex.main=4,margins=c(10,10))

legend("topright",paste("pred.",names(table(com.preN)),sep=""),pch=15,col=RColorBrewer::brewer.pal(n=8,name="Dark2")[2:1],cex=0.6,bty="n")

再對腫瘤細胞再聚類并畫熱圖，又能分成兩群。

table(pred.test$copykat.pred)
tumor.cells<-pred.test$cell.names[which(pred.test$copykat.pred=="aneuploid")]
colnames(CNA.test)<-gsub(".1$","-1",colnames(CNA.test))
tumor.mat<-CNA.test[,colnames(CNA.test)%in%tumor.cells]
hcc<-hclust(parallelDist::parDist(t(tumor.mat),threads=4,method="euclidean"),method="ward.D2")
hc.umap<-cutree(hcc,2)

rbPal6<-colorRampPalette(RColorBrewer::brewer.pal(n=8,name="Dark2")[3:4])
subpop<-rbPal6(2)[as.numeric(factor(hc.umap))]
cells<-rbind(subpop,subpop)

heatmap.3(t(tumor.mat),dendrogram="r",distfun=function(x)parallelDist::parDist(x,threads=4,method="euclidean"),
hclustfun=function(x)hclust(x,method="ward.D2"),
ColSideColors=chr1,RowSideColors=cells,Colv=NA,Rowv=TRUE,
notecol="black",col=my_palette,breaks=col_breaks,key=TRUE,
keysize=1,density.info="none",trace="none",
cexRow=0.1,cexCol=0.1,cex.main=1,cex.lab=0.1,
symm=F,symkey=F,symbreaks=T,cex=1,cex.main=4,margins=c(10,10))

legend("topright",c("c1","c2"),pch=15,col=RColorBrewer::brewer.pal(n=8,name="Dark2")[3:4],cex=0.9,bty='n')

最後把CNV的結果投射到單細胞聚類結果上看一看是否合理，标準流程走一遍，聚類結果和分群結果投射到TSNE上。

standard10X=function(dat,nPCs=50,res=1.0,verbose=FALSE){
srat=CreateSeuratObject(dat)
srat=NormalizeData(srat,verbose=verbose)
srat=ScaleData(srat,verbose=verbose)
srat=FindVariableFeatures(srat,verbose=verbose)
srat=RunPCA(srat,verbose=verbose)
srat=RunTSNE(srat,dims=seq(nPCs),verbose=verbose)
srat=FindNeighbors(srat,dims=seq(nPCs),verbose=verbose)
srat=FindClusters(srat,res=res,verbose=verbose)
return(srat)
}

GC1<-standard10X(counts,nPCs=30,res=0.8)
GC1$copykat.pred<-pred.test$copykat.pred
GC1$copykat.tumor.pred<-rep("normal",nrow(GC1@meta.data))
table(hc.umap)
GC1$copykat.tumor.pred[rownames(GC1@meta.data)%in%names(hc.umap[hc.umap==1])]<-"tumorcluster1"
GC1$copykat.tumor.pred[rownames(GC1@meta.data)%in%names(hc.umap[hc.umap==2])]<-"tumorcluster2"

p1<-DimPlot(GC1,label=T)
p2<-DimPlot(GC1,group.by="copykat.pred")
p3<-DimPlot(GC1,group.by="copykat.tumor.pred")
p1+p2+p3

可以看到：2,4,8群是腫瘤細胞亞克隆

從免疫細胞和腫瘤細胞的标記基因表達來看，可以正确找出正常細胞和腫瘤細胞。

FeaturePlot(GC1,features=c("PTPRC","EPCAM"),order=T)

結語

本期，我們使用分析評估上皮細胞良惡性。下一期，我們将對T細胞亞群進行細分。順便提前預告一下，胃癌系列推文完成後，将開啟肺腺癌單細胞數據集複現系列，相關視頻已經在B站上線：

文獻推文詳見（單細胞測序+空間轉錄組描繪從癌前病變到浸潤性肺腺癌的動态演變 ）。此外，關于推文内容的提升和優化，歡迎大家提寶貴意見。謝謝！