hello，大家好垮卓，隔離的第六天垫桂，一個人的日子總是很難熬，孤獨的時候粟按，還是很希望大家多與我交流交流單細(xì)胞空間的內(nèi)容诬滩，以此來排解孤獨的情緒，今天周六灭将，我們不講新的內(nèi)容疼鸟，回顧一下推斷CNV的分析軟件CopyCAT，之前呢庙曙，分享過一次空镜，文章在copyKAT推斷單細(xì)胞轉(zhuǎn)錄組腫瘤細(xì)胞CNV（自動識別腫瘤normal和tumor），這一篇我們來回顧一下這個軟件，包括算法和示例代碼吴攒。

單細(xì)胞轉(zhuǎn)錄組分析廣泛用于研究腫瘤张抄。然而，將腫瘤微環(huán)境中的正常細(xì)胞類型與惡性細(xì)胞區(qū)分開來并解析腫瘤內(nèi)的克隆亞結(jié)構(gòu)仍然具有挑戰(zhàn)性（目前我們運用最多的還是inferCNV）洼怔。為了應(yīng)對這些挑戰(zhàn)署惯，開發(fā)了一種集成的貝葉斯分割方法，稱為CopyKAT镣隶，以從高通量單細(xì)胞 RNA 測序數(shù)據(jù)中的讀取深度估計平均基因組分辨率為 5 Mb 的基因組拷貝數(shù)譜（這個地方的數(shù)據(jù)我們還是要詳細(xì)說明一下）极谊。

Overview of CopyKAT workflow

scRNA-seq 方法已成為描繪腫瘤微環(huán)境 (TME) 中正常細(xì)胞類型和了解各種癌癥中腫瘤細(xì)胞表達(dá)程序的有力工具。然而矾缓，分析大規(guī)模數(shù)據(jù)集的一個主要挑戰(zhàn)是將腫瘤細(xì)胞與 TME 中的基質(zhì)細(xì)胞和免疫細(xì)胞區(qū)分開來怀酷，以便可以獨立研究它們。區(qū)分腫瘤細(xì)胞與正常細(xì)胞的一種有效方法是鑒定非整倍體拷貝數(shù)譜嗜闻，這在大多數(shù)腫瘤中很常見（88%）蜕依，而在具有二倍體基因組的基質(zhì)細(xì)胞類型中不存在。以前的方法琉雳，例如 inferCNV 和 HoneyBadger样眠，已經(jīng)表明可以根據(jù)足夠大的基因組區(qū)域的 RNA 讀取計數(shù)來估計基因組拷貝數(shù)譜。然而翠肘，這些方法是為分析來自第一代 scRNA-seq 技術(shù)的數(shù)據(jù)而設(shè)計的檐束，這些技術(shù)具有較低的細(xì)胞通量和較高的覆蓋深度（這個地方要注意一下）。這些方法不適用于分析來自新開發(fā)的高通量 scRNA-seq 平臺（微滴和納米孔平臺）的數(shù)據(jù)束倍，這些平臺執(zhí)行全轉(zhuǎn)錄組擴(kuò)增和僅在非常稀疏的覆蓋深度下對 mRNA 的 3' 或 5' 端進(jìn)行測序被丧。此外，以前的方法不能準(zhǔn)確地解決特定染色體斷點的基因組位置或從非整倍體拷貝數(shù)譜中對腫瘤和正常細(xì)胞進(jìn)行分類绪妹。為了應(yīng)對這些挑戰(zhàn)甥桂，開發(fā)了 CopyKAT 并將其應(yīng)用于各種腫瘤樣本，以識別非整倍體腫瘤細(xì)胞并描繪腫瘤塊內(nèi)共存的不同亞群的克隆亞結(jié)構(gòu)邮旷。

總結(jié)一下inferCNV的分析缺點（這個是文章中作者的觀點）

• 為分析來自第一代 scRNA-seq 技術(shù)的數(shù)據(jù)而設(shè)計的黄选，技術(shù)具有較低的細(xì)胞通量和較高的覆蓋深度
• 不適用于分析來自新開發(fā)的高通量 scRNA-seq 平臺（微滴和納米孔平臺）的數(shù)據(jù)，這些平臺執(zhí)行全轉(zhuǎn)錄組擴(kuò)增和僅在非常稀疏的覆蓋深度下對 mRNA 的 3' 或 5' 端進(jìn)行測序（10X的單細(xì)胞技術(shù)具有這個特點）
• 不能準(zhǔn)確地解決特定染色體斷點的基因組位置或從非整倍體拷貝數(shù)譜中對腫瘤和正常細(xì)胞進(jìn)行分類（這個地方的難度相當(dāng)高婶肩，我們來看看CopyCAT如何解決這個問題）办陷。

Overview of CopyKAT workflow.

圖片.png

為了檢測染色體斷點（chromosome breakpoints,）椎瘟，我們整合了泊松伽馬模型和馬爾可夫鏈蒙特卡洛 (MCMC) 迭代來生成每個基因窗口的后驗均值，然后應(yīng)用 Kolmogorov-Smirnov (KS) 檢驗來加入在它們之間沒有顯著差異的相鄰窗口方法侄旬。為了加快計算速度肺蔚，將數(shù)千個單細(xì)胞分成clusters，找到一致的染色體斷點并將它們合并在一起儡羔，形成樣本中整個細(xì)胞群的基因組斷點的聯(lián)合宣羊。然后將每個窗口的最終拷貝數(shù)值計算為跨越每個細(xì)胞中相鄰染色體斷點的所有基因的后驗平均值。我們通過將基因重新排列到 220-kb 可變基因組bin中汰蜘，進(jìn)一步將得到的拷貝數(shù)值從基因空間轉(zhuǎn)換為基因組位置仇冯，從而以大約 5 Mb 的分辨率獲得每個單細(xì)胞的全基因組拷貝數(shù)譜。基因組分辨率是根據(jù)整個基因組的中位相鄰基因距離（~20 kb）乘以基因窗口的大屑ā（25個基因）來估計的（精度高于inferCNV）。然后我們對單細(xì)胞拷貝數(shù)數(shù)據(jù)進(jìn)行層次聚類澈缺，以確定非整倍體腫瘤細(xì)胞和二倍體基質(zhì)細(xì)胞之間的最大距離坪创；但是，如果基因組距離不顯著姐赡，我們切換到 GMM 定義模型來逐個預(yù)測單個腫瘤細(xì)胞莱预。最后，我們對單細(xì)胞拷貝數(shù)數(shù)據(jù)進(jìn)行聚類以識別克隆亞群并計算代表亞克隆基因型的共有譜项滑，以進(jìn)一步分析它們的基因表達(dá)差異依沮。

算法原理

至于前面的數(shù)據(jù)預(yù)處理和推斷參考細(xì)胞類型的部分相恃，我們不過多討論韩脏，主要就是CNV的推斷。

Copy number segmentation.

為了執(zhí)行segmentation，首先通過應(yīng)用 R 包 MCMCpack 中的 MCpoissongamma 函數(shù)狼纬，使用基于具有伽馬先驗的泊松似然的蒙特卡羅模擬來轉(zhuǎn)換clusters一致性以在細(xì)胞clusters內(nèi)找到染色體斷點。 假設(shè)分段均值的泊松分布和泊松參數(shù)的共軛伽馬分布如下：

圖片.png

其中 y 是反向轉(zhuǎn)換的 UMI 計數(shù)蛮位，α 和 β 是伽馬分布的形狀參數(shù)穷吮。 為每個窗口模擬 1,000 個后驗均值，然后使用 Kolmogorov-Smirnov (KS) 檢驗來確定是否應(yīng)將相鄰窗口連接在一起薄嫡。 將 KS 檢驗統(tǒng)計量 D 計算為兩個相鄰窗口的后驗均值累積分布之間的最大垂直距離：

圖片.png

其中 x 對應(yīng)于相鄰窗口 i 和 j 的后驗均值氧急。 如果 KS 檢驗統(tǒng)計量 D 大于截止值，則定義斷點毫深。 如果在第一輪中使用初始截止值定義的斷點少于 25 個吩坝，將截止值降低 50%。 我們對每個clusters最多重復(fù)這個過程兩次哑蔫。 統(tǒng)一來自每個clusters的所有斷點钉寝，以形成整個群體的共識斷點。 然后鸳址，通過將 MCpoissongamma 函數(shù)應(yīng)用于每個單個細(xì)胞的所有共識窗口來計算相鄰斷點之間每個窗口的后驗均值作為分段均值瘩蚪。 然后，對平均值進(jìn)行對數(shù)轉(zhuǎn)換稿黍，并將數(shù)據(jù)居中作為每個基因的最終相對拷貝數(shù)比率疹瘦。 最后，通過計算落入human基因組中 220kb 基因組bin的基因的平均值巡球，將單個基因拷貝數(shù)轉(zhuǎn)換為基因組bin言沐，以根據(jù) scRNA-seq 數(shù)據(jù)估計每個細(xì)胞的全基因組拷貝數(shù)譜。

Theoretical estimation of genomic copy number resolution

為了估計從單細(xì)胞 RNA 數(shù)據(jù)推斷出的拷貝數(shù)譜的預(yù)期分辨率酣栈，下載了 GRCh38 (v28) 中所有基因條目的 BED 文件险胰。因為染色體 Y 不包括在拷貝數(shù)計算中，只考慮了位于染色體 1-22 和染色體 X 上的基因矿筝，它們共有 56,051 個基因條目起便。通過取基因起始位置和基因結(jié)束位置的平均值來估計單個基因的基因組中心位置。接下來窖维，根據(jù)基因組位置對所有基因進(jìn)行排序榆综，并通過計算兩個基因中心之間的距離來估計兩個相鄰基因之間的距離≈罚總的來說鼻疮，在整個基因組中定義了 56,028 個基因區(qū)間。從染色體 1-22 和染色體 X 中琳轿，基因區(qū)間的數(shù)量如下：5,127, 3,872, 2,925, 2,430, 2,779, 2,802, 2,292, 2,189, 2,137, 3,189, 2,857, 1,279, 2,152, 2,081, 2,440, 1,133判沟、2,917耿芹、1,350、795挪哄、1,300 和 2,281吧秕。整個基因組中基因間隔的第一四分位數(shù)、中位數(shù)中燥、平均值寇甸、第三四分位數(shù)和最大值如下：9,430 bp、24,532 bp疗涉、52,806 bp拿霉、58,485 bp 和 21,765,992 bp。因為基因區(qū)間的大小分布嚴(yán)重向右傾斜咱扣，計算了中值來估計拷貝數(shù)分辨率绽淘。因為需要在pipeline中的整個單細(xì)胞群中檢測到至少 7,000 個基因，所以這個數(shù)字相當(dāng)于基因檢測率的中位數(shù) 7,000/56,051 ≈ 12.5%闹伪。最后沪铭，將分析中的最小基因間隔計算為每個基因間隔 24,532 bp ÷ 12.5% ≈ 200 kb。我們使用 25 個基因窗口啟動了我們的拷貝數(shù)分析偏瓤；因此杀怠，估計片段的最小大小為 200 kb × 25 = 5 Mb，用于檢測每個細(xì)胞基因組中的拷貝數(shù)事件的基因組分辨率厅克。

示例代碼及結(jié)果解讀

installation

library(devtools)
install_github("navinlabcode/copykat")

注意輸入的矩陣是測序的原始矩陣

copykat.test <- copykat(rawmat=exp.rawdata, sam.name=“test”)

prepare readcount input

準(zhǔn)備運行 copykat 的唯一一個直接輸入是原始基因表達(dá)矩陣赔退，其中行中包含基因，列中包含細(xì)胞名稱证舟。 基因 id 可以是基因符號或 ensemble id硕旗。 矩陣值通常是來自當(dāng)今高通量單細(xì)胞 RNAseq 數(shù)據(jù)的唯一分子標(biāo)識符 (UMI) 的計數(shù)。 早期生成的 scRNAseq 數(shù)據(jù)可以總結(jié)為 TPM 值或總讀取計數(shù)女责，這也應(yīng)該有效漆枚。 下面提供了一個從 10X 輸出生成這個 UMI 計數(shù)矩陣的示例。

An example to generate input from 10X genomics cellranger V3 output

library(Seurat)
raw <-Read10X(data.dir =data.path.to.cellranger.outs)
raw <- CreateSeuratObject(counts =raw,project ="copycat.test",min.cells = 0,min.features = 0)
exp.rawdata <- as.matrix(raw@assays$RNA@counts)

run copykat

已經(jīng)準(zhǔn)備好唯一的一個輸入抵知，原始 UMI 計數(shù)矩陣墙基，準(zhǔn)備運行 copykat。 默認(rèn)的基因 ids in cellranger 輸出是基因符號刷喜，所以輸入“符號”或“S”残制。 為了過濾掉細(xì)胞，需要每條染色體中至少有 5 個基因來計算 DNA 拷貝數(shù)吱肌。 可以將它調(diào)低到 ngene.chr=1 以保持盡可能多的細(xì)胞痘拆，但是認(rèn)為使用至少 5 個基因來代表一個染色體并不是很嚴(yán)格仰禽。 為了過濾外基因氮墨，可以調(diào)整參數(shù)以僅保留在 LOW.DR 到 UP.DR 部分細(xì)胞中表達(dá)的基因纺蛆。 這里設(shè)置默認(rèn) LOW.DR=0.05，UP.DR=0.2规揪。 可以調(diào)低這些值以在分析中保留更多基因桥氏。 不過，需要確保 LOW.DR 小于 UP.DR猛铅。

要求 copykat 每個片段至少取 25 個基因字支。 也可以使用其他選項，范圍為每箱bin 15-150 個基因奸忽。 KS.cut 是分段參數(shù)堕伪，范圍從 0 到 1。增加 KS.cut 會降低靈敏度栗菜，即減少分段/斷點欠雌。

一項艱難而有趣的觀察是，沒有一種聚類方法可以適合所有數(shù)據(jù)集疙筹。 在這個版本中富俄，為聚類添加了一個距離參數(shù)，包括“歐幾里得”距離和相關(guān)距離而咆，即霍比。 1-“pearson”和“spearman”的相似性。 一般來說暴备，校正距離傾向于有利于噪聲數(shù)據(jù)悠瞬，而歐幾里德距離傾向于有較大 CN 段的數(shù)據(jù)。

library(copykat)
copykat.test <- copykat(rawmat=exp.rawdata,id.type="S",cell.line="no",ngene.chr= 5,win.size=25,KS.cut=0.2,sam.name="test",distance="euclidean",n.cores=4)
pred.test <- data.frame(copykat.test$prediction)
CNA.test <- data.frame (copykat.test$CNAmat)

navigate prediction results

Predicted aneuploid cells are inferred as tumor cells; diploid cells arestromal normal cells.

圖片.png

The first 3 columns in the CNA matrix are the genomic coordinates. Rows are 220KB bins in genomic orders.

圖片.png

熱圖：

圖片.png

define subpopulations of aneuploid tumor cells

tumor.cells <- pred.test$cell.names[which(pred.test$copykat.pred=="aneuploid")]
tumor.mat <- CNA.test[, which(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')

圖片.png

現(xiàn)在已經(jīng)定義了兩個具有主要亞克隆差異的亞群馍驯，可以繼續(xù)比較它們的基因表達(dá)譜阁危，并評估亞克隆拷貝數(shù)變化對基因表達(dá)的影響。

最后一點汰瘫，一些有用的注釋數(shù)據(jù)和copykat：

• full.anno：56051個基因（hg38）的完整注釋狂打，包括絕對位置、chr混弥、start趴乡、end ensemble id、基因符號和Gband蝗拿。
• DNA.hg20：不包括 Y 染色體的 220kb 變量 bin 在 hg38 中的坐標(biāo)晾捏。
• cyclegenes：從CNV分析中刪除。
• exp.rawdata：來自乳腺腫瘤的 UMI 矩陣哀托。

最后惦辛，CopyKAT 難以預(yù)測具有少量 CNA 的pediatric和liquid tumors的腫瘤和正常細(xì)胞。 CopyKAT 提供了兩種方法來繞過它以提供特定的輸出：1）輸入來自同一數(shù)據(jù)集的已知正常細(xì)胞的細(xì)胞名稱向量 2）或嘗試搜索 T 細(xì)胞仓手。（看來可以用免疫細(xì)胞作為ref）胖齐。

好了玻淑，生活很好，有你更好呀伙，希望疫情趕緊結(jié)束