CIBERSORT免疫浸潤(rùn)的代碼實(shí)現(xiàn) (qq.com)
數(shù)據(jù)直接從xena上下載贴捡，命名為TCGA_chol
rm(list=ls())
library(data.table)
#help(fread)
TGCA_CHOL<-fread("TCGA_chol",data.table=F)
rownames(TGCA_CHOL)<-TGCA_CHOL[,1]
TGCA_CHOL<-TGCA_CHOL[,-1]
exp<-TGCA_CHOL
exp<-2^TGCA_CHOL-1
head(exp)
按照上面的教程走豪筝，但是走了一點(diǎn)彎路遗契，事實(shí)上這些如果是從xena上下載之后，其實(shí)不需要轉(zhuǎn)換為TPM的浊伙，直接使用FPKM就行，但是還是保存下來出皇，因?yàn)橛行┐a還是廢了不少腦筋去自己寫
gencode.v22.annotation.gtf.gz
# if(T){
# library(rtracklayer)
#
# gtf = fread("gencode.v25.annotation.gff3.gz")
# class(gtf)
# gtf = as.data.frame(gtf);dim(gtf)
# table(gtf$type)
# colnames(gtf)<-c("chrom","speices","type","start","end","strand1","strand2","strand3","gene_name")
# exon = gtf[gtf$type=="exon",
# c("start","end","gene_name")]
# s=as.data.frame(strsplit(exon$gene_name,";"),fill=T)
# n=split(exon,exon$gene_name)
# n2<-as.data.frame(n)
# head(n2)
# gle = lapply(split(exon,exon$gene_name),function(x){
# #exon$name=strsplit(strsplit(exon$gene_name,";")[[1]][7],"=")[[1]][2]
# tmp=apply(x,1,function(y){
# y[1]:y[2]
# })
# length(unique(unlist(tmp)))
# })
# gle=data.frame(gene_name=names(gle),
# length=as.numeric(gle))
# head(gle)
# library(tidyverse)
# gle %>% separate(gene_name,c("key","value"),"gene_name=")%>%separate(value,c("gene","b"),";")->gle2
# gle3<-gle2[,c("gene","length")]
# head(gle3)
# save(gle3,file = "v25_gle.Rdata")
# }
#
# sub=function(x){
# strsplit(strsplit(x,";")[[1]][7],"=")[[1]][2]
# }
# bb<-unlist(lapply(exon$gene_name, sub))
# bb2<-length(unique(unlist(bb)))
# bb2
#le = gle[match(rownames(exp),gle3$gene_name),"length"]
# if(T){
# #head(gtf)
# library(rtracklayer)
# gtf = rtracklayer::import("gencode.v22.annotation.gtf.gz")
# class(gtf)
# gtf = as.data.frame(gtf);dim(gtf)
# table(gtf$type)
# exon = gtf[gtf$type=="exon",
# c("start","end","gene_name")]
# gle = lapply(split(exon,exon$gene_name),function(x){
# tmp=apply(x,1,function(y){
# y[1]:y[2]
# })
# length(unique(unlist(tmp)))
# })
# gle=data.frame(gene_name=names(gle),
# length=as.numeric(gle))
# save(gle,file = "v22_gle.Rdata")
# }
# load("v22_gle.Rdata")
# head(gle)
# ```
```{r}
# le = gle[match(rownames(exp),gle$gene_name),"length"]
# head(le)
# #這個(gè)函數(shù)是現(xiàn)成的钠署。
# countToTpm <- function(counts, effLen)
# {
# rate <- log(counts) - log(effLen)
# denom <- log(sum(exp(rate)))
# exp(rate - denom + log(1e6))
# }
# tpms <- apply(exp,2,countToTpm,le)
# tpms[1:3,1:3]
直接從這里開始，上面的可以不要管
source("CIBERSORT.R")
exp2<-as.data.frame(exp)
exp2<-data.frame(gene=rownames(exp2),exp2)
head(exp2)
write.table(exp2,file = "exp.txt",row.names = F,quote = F,sep = "\t")
if(T){
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]
re <- TME.results[,-(23:25)]
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))
library(tidyverse)
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))
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))
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))
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")
