細(xì)胞免疫浸潤網(wǎng)站使用
1.cibersort
數(shù)據(jù)輸入格式：輸入數(shù)據(jù)的要求： 1.不可以有負(fù)值和缺失值 2.不要取log 3.如果是芯片數(shù)據(jù)摊阀，昂飛芯片使用RMA標(biāo)準(zhǔn)化人断，Illumina 的Beadchip 和Agilent的單色芯片吭从，用limma處理。 4.如果是RNA-Seq表達(dá)量恶迈，使用FPKM和TPM都很合適涩金。芯片的要求可能把你唬住了，GEO常規(guī)的表達(dá)矩陣都是這樣得到的暇仲，直接下載使用即可步做。注意有的表達(dá)矩陣下載下來就已經(jīng)取過log，需要逆轉(zhuǎn)回去奈附。有的經(jīng)過了標(biāo)準(zhǔn)化或者有負(fù)值全度，需要處理原始數(shù)據(jù)，前面寫過介紹文：
http://www.reibang.com/p/d7035ba8347b
http://www.reibang.com/p/e3d734b2c404
做成cibersort要求的輸入文件
這個(gè)算法并沒有被寫成R包斥滤，而是只有一個(gè)放著函數(shù)的腳本–CIBERSORT.R将鸵，把它下載下來放在工作目錄即可。
需要兩個(gè)輸入文件:
一個(gè)是表達(dá)矩陣文件
一個(gè)是官網(wǎng)提供的LM22.txt佑颇，記錄了22種免疫細(xì)胞的基因表達(dá)特征數(shù)據(jù)咨堤。三個(gè)輸入文件的要求及形式，參考以下網(wǎng)址：
http://www.reibang.com/p/0baac4c52ac8
由于CIBERSORT.R讀取文件的代碼比較粗暴漩符，為了適應(yīng)它，導(dǎo)出文件之前需要把行名變成一列驱还。不然后面就會(huì)有報(bào)錯(cuò)嗜暴。
exp2=as.data.frame(tpms)
exp2=rownames_to_column(exp2)
write.table(exp2,file="exp.txt",row.names=F,quote=F,sep="\t")
source("CIBERSORT.R")
if(F){
? TME.results = CIBERSORT("LM22.txt",
? ? ? ? ? ? ? ? ? ? ? ? ? "exp.txt" ,
? ? ? ? ? ? ? ? ? ? ? ? ? perm = 1000,
? ? ? ? ? ? ? ? ? ? ? ? ? QN = T)
? save(TME.results,file = "ciber_CHOL.Rdata")
}
load("ciber_CHOL.Rdata")
TME.results[1:4,1:4]
## B cells naive B cells memory Plasma cells
## TCGA-W5-AA36-01A-11R-A41I-07 0.00000000 0.002351185 0.02550133
## TCGA-W5-AA2H-01A-31R-A41I-07 0.04512086 0.354414124 0.01961627
## TCGA-ZU-A8S4-11A-11R-A41I-07 0.00203370 0.000000000 0.04582565
## TCGA-WD-A7RX-01A-12R-A41I-07 0.15785229 0.000000000 0.01847074
## T cells CD8
## TCGA-W5-AA36-01A-11R-A41I-07 0.07766099
## TCGA-W5-AA2H-01A-31R-A41I-07 0.14262301
## TCGA-ZU-A8S4-11A-11R-A41I-07 0.09962641
## TCGA-WD-A7RX-01A-12R-A41I-07 0.13769951
re <- TME.results[,-(23:25)]
運(yùn)行有些慢。計(jì)算出來的結(jié)果包含了22種免疫細(xì)胞的豐度议蟆，還有三列其他統(tǒng)計(jì)量闷沥，不管它們。
3.5. 經(jīng)典的免疫細(xì)胞豐度熱圖
那些在一半以上樣本里豐度為0的免疫細(xì)胞咐容，就不展示在熱圖里了舆逃。我看了一下這個(gè)熱圖，從聚類的情況來看戳粒，normal和tumor沒有很好的分開路狮。
library(pheatmap)k<-apply(re,2,function(x){sum(x==0)<nrow(TME.results)/2})table(k)
## k## FALSE? TRUE ##? ? 8? ? 14
re2<-as.data.frame(t(re[,k]))an=data.frame(group=Group,row.names=colnames(exp))pheatmap(re2,scale="row",show_colnames=F,annotation_col=an,color=colorRampPalette(c("navy","white","firebrick3"))(50))
3.6. 直方圖
可以展示出每個(gè)樣本的免疫細(xì)胞比例
library(RColorBrewer)
mypalette<-colorRampPalette(brewer.pal(8,"Set1"))
dat<-re%>%as.data.frame()%>%
rownames_to_column("Sample")%>%
gather(key=Cell_type,value=Proportion,-Sample)
ggplot(dat,aes(Sample,Proportion,fill=Cell_type))+
geom_bar(stat="identity")+labs(fill="Cell Type",x="",y="Estiamted Proportion")+
theme_bw()+
theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),legend.position="bottom")+
scale_y_continuous(expand=c(0.01,0))+
scale_fill_manual(values=mypalette(22))
3.7 箱線圖
展示免疫細(xì)胞之間的比較。
ggplot(dat,aes(Cell_type,Proportion,fill=Cell_type))+
geom_boxplot(outlier.shape=21,color="black")+
theme_bw()+labs(x="Cell Type",y="Estimated Proportion")+
theme(axis.text.x=element_blank(),axis.ticks.x=element_blank(),legend.position="bottom")+
scale_fill_manual(values=mypalette(22))
亂了點(diǎn)蔚约？那就讓箱線圖擁有順序吧奄妨。
a = dat %>%? group_by(Cell_type) %>%? summarise(m = median(Proportion)) %>%? arrange(desc(m)) %>%? pull(Cell_type)dat$Cell_type = factor(dat$Cell_type,levels = a)
ggplot(dat,aes(Cell_type,Proportion,fill = Cell_type)) +? geom_boxplot(outlier.shape = 21,color = "black") +??
?theme_bw() +??
?labs(x = "Cell Type", y = "Estimated Proportion") +? ?
?theme(axis.text.x = element_blank(),? ? ? ??
axis.ticks.x = element_blank(),? ? ? ??
legend.position = "bottom") +?
?scale_fill_manual(values = mypalette(22))
既然我們已經(jīng)把正常樣本也算了，那就做個(gè)比較：
dat$Group = ifelse(as.numeric(str_sub(dat$Sample,14,15))<10,"tumor","normal")
library(ggpubr)
ggplot(dat,aes(Cell_type,Proportion,fill = Group)) +? geom_boxplot(outlier.shape = 21,color = "black") +??
?theme_bw() +? labs(x = "Cell Type", y = "Estimated Proportion") +? theme(legend.position = "top") +??
?theme(axis.text.x = element_text(angle=80,vjust = 0.5))+? scale_fill_manual(values = mypalette(22)[c(6,1)])+ stat_compare_means(aes(group = Group,label = ..p.signif..),method = "kruskal.test")
2.TImer2.0
腫瘤微環(huán)境中免疫細(xì)胞的組成和豐度是影響腫瘤進(jìn)展和免疫治療效果的重要因素苹祟。由于直接測(cè)量方法的局限性砸抛，計(jì)算算法通常用于從腫瘤轉(zhuǎn)錄組圖譜推斷免疫細(xì)胞組成评雌。這些估計(jì)的腫瘤免疫浸潤人群與腫瘤的基因組和轉(zhuǎn)錄組學(xué)變化相關(guān)，提供了對(duì)腫瘤免疫相互作用的見解直焙。今天我們就一起來看這個(gè)腫瘤免疫分析的又一神器——TIMER 2.0數(shù)據(jù)庫（開始的頁面就已經(jīng)有數(shù)據(jù)格式要求景东，注意看）。具體操作如下：
https://zhuanlan.zhihu.com/p/367363990
http://www.sci666.net/59270.html
https://www.jingege.wang/2021/09/27/%E7%A0%94%E7%A9%B6%E5%85%8D%E7%96%AB%E6%B5%B8%E6%B6%A6%E7%9A%84%E7%A5%9E%E5%99%A8timer%E6%9B%B4%E6%96%B0-timer2-0/
3.xcell
http://www.reibang.com/p/a20a534233aa
4.ESTIMATE
腫瘤微環(huán)境中奔誓，免疫細(xì)胞和基質(zhì)細(xì)胞是兩種主要肺腫瘤細(xì)胞類型斤吐；而ESTIMATE，利用表達(dá)譜數(shù)據(jù)來預(yù)測(cè)基質(zhì)細(xì)胞和免疫細(xì)胞評(píng)分丝里，進(jìn)而預(yù)測(cè)這兩種細(xì)胞的含量曲初；最后可以計(jì)算每個(gè)腫瘤樣本里面的腫瘤純度了，也就是說杯聚，基質(zhì)細(xì)胞和免疫細(xì)胞含量多臼婆，那腫瘤純度低，反之腫瘤純度就高了幌绍。
可能對(duì)RNA-Seq數(shù)據(jù)的預(yù)測(cè)不是很準(zhǔn)颁褂，芯片數(shù)據(jù)還可以。
代碼參考：http://www.reibang.com/p/b64e9f338f8b
另一種分析方法網(wǎng)頁+代碼參考：https://zhuanlan.zhihu.com/p/404707114
https://zhuanlan.zhihu.com/p/136747705
整理的代碼及注意事項(xiàng)
http://www.reibang.com/p/ec5307256ca5
library(utils)rforge<-"http://r-forge.r-project.org"
install.packages("estimate",repos=rforge,dependencies=TRUE)
library(estimate)
初步運(yùn)行自帶的測(cè)試數(shù)據(jù)
library(estimate)
OvarianCancerExpr<system.file("extdata","sample_input.txt",package="estimate")
filterCommonGenes(input.f=OvarianCancerExpr,output.f="OV_10412genes.gct",id="GeneSymbol")
estimateScore(input.ds="OV_10412genes.gct",output.ds="OV_estimate_score.gct",platform="affymetrix")
plotPurity(scores="OV_estimate_score.gct",samples="s519",platform="affymetrix")
從測(cè)試數(shù)據(jù)可以看出傀广，設(shè)置好工作路徑后颁独，將.txt文件直接放在當(dāng)前工作目錄下，就可以不跑OvarianCancerExpr伪冰，直接跑filterCommonGenes誓酒，因?yàn)榭梢栽O(shè)置input.f=xxx.txt。
數(shù)據(jù)格式如圖1：
完整代碼
setwd("E:/1.RStudio/geo_ovarian_GPL96_570/geo_array_ovarian/OVCP3")
library(utils)
rforge <- "http://r-forge.r-project.org"
install.packages("estimate", repos=rforge, dependencies=TRUE)
library(estimate)
###提前將Exp0放在當(dāng)前工作目錄下
filterCommonGenes(input.f="Exp0.txt",
? ? ? ? ? ? ? ? ? output.f="OV_6934genes.gct",
? ? ? ? ? ? ? ? ? id="GeneSymbol")
estimateScore(input.ds = "OV_6934genes.gct",
? ? ? ? ? ? ? output.ds="OV_estimate_score.gct",
? ? ? ? ? ? ? platform="affymetrix")
plotPurity(scores="OV_estimate_score.gct",
? ? ? ? ? platform="affymetrix")
###########結(jié)果讀取和導(dǎo)出方式
scores=read.table("OV_estimate_score.gct",skip = 2,header = T)
rownames(scores)=scores[,1]
scores=t(scores[,3:ncol(scores)])
View(scores)
save(scores,file = 'BRCA_estimate_score.Rdata')
