劉小澤寫于19.9.6-第四單元第二講：評估任意基因集在癌癥的表現(xiàn)
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53

前言

上一篇是探索兩個細(xì)胞亞群（vCAF、mCAF）特有的基因在TCGA中的表現(xiàn)，發(fā)現(xiàn)兩個亞群的基因都是和TCGA
相關(guān)的基因在內(nèi)部相關(guān)，說明了分群的效果不錯

目的就是做下面這個圖的相關(guān)性分析：

可以看到，橫坐標(biāo)的vCAF就是我們前一篇得到的vCAF基因集在TCGA數(shù)據(jù)集中的表達(dá)量，那么縱坐標(biāo)，就需要去文章里找，作者是拿到了5篇不同參考文獻(xiàn)的6個數(shù)據(jù)集

文章正文放了四張相關(guān)性的圖，是vCAF和mCAF與第27篇參考文獻(xiàn)中的兩個乳腺癌數(shù)據(jù)集進(jìn)行的比較

然后再來看看第27篇文獻(xiàn)的圖，其中列出了乳腺癌的ECM和Endothelial的基因集

然后這篇參考文獻(xiàn)的作者定義基因集的方法就是：在大部分癌癥中都存在的基因就是基因集，因此我們看到，本文使用的基因集中就3個基因

看到其中有一個奇怪的基因名CXorf36，按說人類的基因名都應(yīng)該是大寫，所以拿到GenCard查詢一下：https://www.genecards.org/cgi-bin/carddisp.pl?gene=DIPK2B&keywords=CXorf36。發(fā)現(xiàn)這個基因目前叫DIPK2B

有了基因集，就去獲取各個基因集的表達(dá)量

各個基因如下：

library(stringr)
# vCAF基因集
vCAF='Esam, Gng11, Higd1b, Cox4i2, Cygb, Gja4, Eng'
vCAF=unlist(str_split(vCAF,', '))
# mCAF基因集
mCAF='Dcn, Col12a1, Mmp2, Lum, Mrc2, Bicc1, Lrrc15, Mfap5, Col3A1, Mmp14, Spon1, Pdgfrl, Serpinf1, Lrp1, Gfpt2, Ctsk, Cdh11, Itgbl1, Col6a2, Postn, Ccdc80, Lox, Vcan, Col1a1, Fbn1, Col1a2, Pdpn, Col6a1, Fstl1, Col5a2, Aebp1'
mCAF=unlist(str_split(mCAF,', '))
# ECM基因集
ECM=c('COL1A1', 'COL1A2','COL3A1')
# Endothelial基因集
endothelial=c('CDH5', 'DIPK2B','TIE1')

還是使用上一篇的GDC乳腺癌TCGA的表達(dá)矩陣（60,489 identifiers X 1217 samples）

需要注意的是：上一篇我們單純比較多個基因相關(guān)性，所以得到多個基因的表達(dá)量然后做個熱圖就好；但是這次要比較的是兩個基因集（每一組內(nèi)都有不同數(shù)量的基因）。作者用散點圖來展現(xiàn)，其中的每一個點實際上就是一個樣本，但是同一個樣本在兩個基因集中對應(yīng)的基因數(shù)量不同，不能簡單拿任何一個基因進(jìn)行比較。作者給的方法是：

但是這里我們只是簡單對每個樣本取個均值

# 讀入數(shù)據(jù)TCGA-BRCA.htseq_counts.tsv.gz
library(data.table)
filepath <- file.choose()
a=fread(filepath ,data.table=F)

# Ensembl ID切割
library(stringr)
esid=str_split(a$Ensembl_ID,
                 '[.]',simplify = T)[,1]
# ID轉(zhuǎn)換
e2s=select(org.Hs.eg.db,keys = esid,columns = c( "ENSEMBL" ,  "SYMBOL" ),keytype = 'ENSEMBL')
vCAF=toupper(vCAF);vCAF=vCAF[vCAF %in% e2s$SYMBOL,]
mCAF=toupper(mCAF);mCAF=mCAF[mCAF %in% e2s$SYMBOL,]

# 獲得表達(dá)量
rownames(a)=esid
a=a[,-1]

ng=e2s[match(vCAF,e2s$SYMBOL),1]
vCAF_value=colMeans(a[ng,])
ng=e2s[match(mCAF,e2s$SYMBOL),1]
mCAF_value=colMeans(a[ng,])
ng=e2s[match(ECM,e2s$SYMBOL),1]
ECM_value=colMeans(a[ng,])
ng=e2s[match(endothelial,e2s$SYMBOL),1]
endothelial_value=colMeans(a[ng,])

dat=data.frame(vCAF_value=vCAF_value,
                 mCAF_value=mCAF_value,
                 ECM_value=ECM_value,
                 endothelial_value=endothelial_value )

最后繪制散點圖

library(ggpubr)
colnames(dat)
ggscatter(dat, x = "vCAF_value", y = "endothelial_value",
          color = 'black', shape = 21, size = 0.5, # Points color, shape and size
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE, # Add confidence interval
          cor.coef = TRUE,  
          cor.coeff.args = list(method = "pearson",  label.sep = "\n")
)