2024-06-16

干货批量处理泛癌相关微卫星不稳定MSI分析

小师妹生信果 2023-09-26 19:00:18

转自公众号：R语言学徒
http://mp.weixin.qq.com/s?__biz=Mzg5MDk3Mzg4OA==&mid=2247486580&idx=1&sn=233c555aac94fe738e655ff68363400d

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]##namestmp<- 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 normalgroup=sapply(strsplit(row.names(exp1),"\-"),"[",4)group=sapply(strsplit(group,""),"[",1)group=gsub("2","1",group)exp1=exp1[group==0,]#intersectexp1=na.omit(exp1)sameSample=intersect(row.names(MSI),row.names(exp1))MSI=MSI[sameSample,]exp1=exp1[sameSample,]#COR TESToutTab=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$estimatepValue=corT$p.valuesig=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.names=F)   #reada=file1[i]a=gsub("|_MSI_fmsbInput\.txt","",a)

接下来，通过循环遍历file1列表中的每个文件名，读取相应的文件数据，并对数据进行一些处理操作。data使用read.table()函数从文件中读取数据，并将数据存储在一个数据框中。

if (length(data)==0) {next}maxValue=ceiling(max(abs(data))*10)/10data=rbind(rep(maxValue,ncol(data)),rep(-maxValue,ncol(data)),data)

然后，通过对数据进行处理，计算最大值maxValue，并在数据的顶部和底部插入相应的最大值数据，以扩展雷达图的范围。

##colorcolors="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))#printpdf(paste0("./RAD/",dirnam,"_MSI_radar.pdf"),height=7,width=7)radarchart( data, axistype=1 ,pcol=colors, #colorplwd=2 , #lwdplty=1, #dotted or solid linecglcol="grey", #Background line colorcglty=1,caxislabels=seq(-maxValue,maxValue,maxValue/2), # scalecglwd=1.2, #Background line thicknessaxislabcol="green", #colorvlcex=0.8 #size)dev.off()}