干货批量处理泛癌相关微卫星不稳定MSI分析
MSI（Microsatellite Instability）微卫星不稳定性是一种遗传变异现象，指的是肿瘤细胞中微卫星DNA重复序列的插入或缺失，导致某些DNA微卫星的长度发生异常改变。MSI在多种癌症中普遍存在，尤其是结直肠癌。
MSI微卫星不稳定性分析的常用方法是通过PCR扩增肿瘤组织和正常组织中的一组微卫星标记位点，并比较它们的DNA序列长度。通常选择了一组已知的高度多态性的微卫星标记位点，这些位点在正常组织中具有稳定的DNA序列长度。如果在肿瘤组织中发现了与正常组织相比显著的微卫星长度改变（即微卫星不稳定性），则说明肿瘤具有MSI。
library(limma)
library(data.table)
library(tidyverse)
library(ggsignif)
library(RColorBrewer)
library(future.apply)
library(fmsb)
载入R包
files=dir("./TG_Exp/")
#Batch read geneexp files 批量读取基因表达文件
#示例基因表达文件，这是之前根据tcga手动批量整理的，后续有人想了解，可以联系小云。
##下面是编写的一个功能函数来批量读取处理单基因表达文件，函数接受一个目录名作为输入a. 读取位于./TG_Exp/目录中指定的目录名的基因表达文件。 b. 过滤数据，只保留肿瘤样本。 c. 使用数据的第一列作为行名。 d. 去掉表达矩阵中的第一列，因为第一列通常是样本名，不是基因名。 e. 解析目录名，提取出肿瘤类型。 f. 读取MSI（微卫星不稳定）数据，MSI数据的文件路径为"./MSI.txt"。 g. 过滤MSI和基因表达数据，只保留两者共有的样本。
corStatMSI=function(dirname){
exp<-read.table(paste0("./TG_Exp/",dirname),header=T,sep='\t',check.names=F)
exp1=subset(exp,exp$sample_type=="tumor")
rownames(exp1)=exp1[,1]
exp1=exp1[,-1]
##names
tmp<- strsplit(dirname,split=".",fixed=TRUE)
dirnam<-unlist(lapply(tmp,head,1))
###read MSI读取MSI微卫星不稳定数据，下面为数据示范。
MSI=read.table("./MSI.txt", header=T,sep="\t",row.names=1,check.names=F)
rownames(MSI)=substr(rownames(MSI),1,nchar(rownames(MSI))-1) #去最后一个字符A
#pass normal
group=sapply(strsplit(row.names(exp1),"\\-"),"[",4)
group=sapply(strsplit(group,""),"[",1)
group=gsub("2","1",group)
exp1=exp1[group==0,]
#intersect
exp1=na.omit(exp1)
sameSample=intersect(row.names(MSI),row.names(exp1))
MSI=MSI[sameSample,]
exp1=exp1[sameSample,]
#COR TEST
outTab=data.frame()
fmsbTab=data.frame()
tissue=c("ACC", "BLCA", "BRCA", "CESC", "CHOL", "COAD", "DLBC", "ESCA", "GBM", "HNSC", "KICH", "KIRC", "KIRP", "LGG", "LIHC", "LUAD", "LUSC", "MESO", "OV", "PAAD", "PCPG", "PRAD", "READ", "SARC", "SKCM", "STAD", "TGCT", "THCA", "THYM", "UCEC", "UCS", "UVM")
#对肿瘤类型循环
for ( j in 1:length(tissue) ){
cancer=tissue[j]
exp2=exp1[(exp1[,2]==cancer),]
MSI1=MSI[(MSI[,2]==cancer),]
x=as.numeric(MSI1[,1])
y=as.numeric(exp2[,1])
if(length(x)!=length(y)){next}
if(length(unique(exp2[,1]))>3){
corT=cor.test(x,y,method="spearman")
cor=corT$estimate
pValue=corT$p.value
sig=ifelse(pValue<0.001,"***",ifelse(pValue<0.01,"**",ifelse(pValue<0.05,"*"," ")))
outTab=rbind(outTab,cbind(tissue=cancer,cor=cor,pValue=pValue,sig))
fmsbTab=rbind(fmsbTab,cbind(tissue=cancer,cor=cor))
}
else{next}
write.table(outTab,paste0("./corStat_MSI/",dirnam,"_","MSI_corStat.txt"),sep="\t",row.names=F,quote=F)
write.table(t(fmsbTab),paste0("./corStat_MSI/",dirnam,"_","MSI_fmsbInput.txt"),sep="\t",col.names=F,quote=F)
}
}
dd=future_lapply(files,corStatMSI) #多核并行运行
后续可以给大家讲解future_lapply多核并行的操作
下面为批量生成的结果MSI 文件txt，以其中一个单基因示范查看内容，分别展示了相关性和p值。
#批量绘制雷达图，读取前面处理好的数据
首先，使用list.files()函数从指定文件夹中获取符合条件的文件名列表。file1存储了符合".fmsbInput.txt"模式的文件名列表，file2存储了符合".corStat.txt"模式的文件名列表。
file1 <- list.files(paste0("./corStat_MSI"), pattern="*.fmsbInput.txt")
file2 <- list.files(paste0("./corStat_MSI"), pattern="*.corStat.txt")
for (i in 1:length(file1) ){
data=read.table(paste0("./corStat_MSI/",file1[i]),header=T,sep="\t",row.names=1,check.na mes=F) #read
a=file1[i]
a=gsub("|_MSI_fmsbInput\\.txt","",a)
接下来，通过循环遍历file1列表中的每个文件名，读取相应的文件数据，并对数据进行一些处理操作。data使用read.table()函数从文件中读取数据，并将数据存储在一个数据框中。
if (length(data)==0) {next}
maxValue=ceiling(max(abs(data))*10)/10
data=rbind(rep(maxValue,ncol(data)),rep(-maxValue,ncol(data)),data)
然后，通过对数据进行处理，计算最大值maxValue，并在数据的顶部和底部插入相应的最大值数据，以扩展雷达图的范围。
##color
colors="blue"
#name +sig ***
corStat=read.table(paste0("./corStat_MSI/",a,"_MSI_corStat.txt"),header=T,sep="\t",row.names=1,check.names=F)
if (length(corStat)==0) {next}
colnames(data)=paste0(colnames(data),corStat$sig)
接下来，根据需求设置雷达图的颜色colors。然后，使用read.table()函数读取与当前文件对应的".corStat.txt"文件的数据，并存储在corStat数据框中。然后，根据需求对数据框中列的名称进行处理，将其与corStat数据框中的sig列的值进行拼接。接下来，根据文件名设置dirname和dirnam变量。
dirname=file1[i]
tmp<- strsplit(dirname,split="_",fixed=TRUE)
dirnam<-unlist(lapply(tmp,head,1))
pdf(paste0("./RAD/",dirnam,"_MSI_radar.pdf"),height=7,width=7)
radarchart( data, axistype=1 ,
pcol=colors, #color
plwd=2 , #lwd
plty=1, #dotted or solid line
cglcol="grey", #Background line color
cglty=1,
caxislabels=seq(-maxValue,maxValue,maxValue/2), # scale
cglwd=1.2, #Background line thickness
axislabcol="green", #color
vlcex=0.8 #size
)
dev.off()
}
来给大家展示展示结果图吧，使用firefox命令。
那么这就是单基因在泛癌中的MSI微卫星不稳定性结果图。
MSI微卫星不稳定性的分析对于癌症的诊断和预后评估具有重要意义。MSI在一些肿瘤中与预后较好、对某些化疗药物敏感等特点相关。因此，对于结直肠癌等肿瘤，进行MSI微卫星不稳定性分析已成为临床的常规检测之一。
下期将为你带来更多R语言的骚操作技巧，以下推荐的是一个多功能的生信平台。
云生信平台链接：http://www.biocloudservice.com/home.html。
云生信平台链接：http://www.biocloudservice.com/home.html。
